NYC Motor Vehicle Collisions – Street-Level Heat Map

StreetLevelMap In this post I will extend a previous analysis creating a borough-level heat map of NYC motor vehicle collisions. The data is from NYC Open Data. In particular, I will go from borough-level to street-level collisions. The processing of the code is very similar to the previous analysis, with a few more functions that map streets to colors. Below, I load the ggmap package, and the data, and only keep collisions with longitude and latitude information.

library(ggmap)

d=read.csv('.../NYPD_Motor_Vehicle_Collisions.csv')
d_clean=d[which(regexpr(',',d$LOCATION)!=-1),]

#### 1. Clean Data ####
# get long and lat coordinates from concatenated "location" var
comm=regexpr(',',d_clean$LOCATION)
d_clean$loc=as.character(d_clean$LOCATION)
d_clean$lat=as.numeric(substr(d_clean$loc,2,comm-1))
d_clean$long=as.numeric(substr(d_clean$loc,comm+1,nchar(d_clean$loc)-1))

# create year variable
d_clean$year=substr(d_clean$DATE,7,10)

Continue reading NYC Motor Vehicle Collisions – Street-Level Heat Map

Simulating Endogeneity

Introduction The topic in this post is endogeneity, which can severely bias regression estimates. I will specifically simulate endogeneity caused by an omitted variable. In future posts in this series, I’ll simulate other specification issues such as heteroskedasticity, multicollinearity, and collider bias.

The Data-Generating Process

Consider the data-generating process (DGP) of some outcome variable Y :

Y = a+\beta x+c z + \epsilon_1

\epsilon_1 \sim N(0,\sigma^{2})

For the simulation, I set parameter values for a , \beta , and c and simulate positively correlated independent variables, x and z (N=500).

# simulation parameters
set.seed(144);
ss=500; trials=5000;
a=50; b=.5; c=.01;
d=25; h=.9;

# generate two independent variables
x=rnorm(n=ss,mean=1000,sd=50);
z=d+h*x+rnorm(ss,0,10)

The Simulation

The simulation will estimate the two models below. The first model is correct in the sense that it includes all terms in the actual DGP. However, the second model omits a variable that is present in the DGP. Instead, the variable is obsorbed into the error term \epsilon_1 .

(1)\thinspace Y = a+\beta x+c z + \epsilon_1

(2) \thinspace Y = a+\beta x + \epsilon_1

This second model will yield a biased estimator of \beta . The variance will also be biased. This is because x is endogenous, which is a fancy way of saying it is correlated with the error term, \epsilon_1 . Since cor(x,z)>0 and \epsilon_1=\epsilon + cz , then cor(x,\epsilon_1)>0 . To show this, I run the simulation below with 5,000 iterations. For each iteration, I construct the outcome variable, Y using the DGP. I then run a regression to estimate \beta , first with model 1 and then with model 2.

sim=function(endog){
 # assume normal error with constant variance to start
 e=rnorm(n=ss,mean=0,sd=10)
 y=a+b*x+c*z+e
 # Select data generation process
 if(endog==TRUE){ fit=lm(y~x) }else{ fit=lm(y~x+z)}
 return(fit$coefficients)
}

# run simulation - with and wihtout endogeneity
sim_results=t(replicate(trials,sim(endog=FALSE)))
sim_results_endog=t(replicate(trials,sim(endog=TRUE)))

Simulation Results This simulation yields two different sampling distributions for \beta . Note that I have set the true value to \beta=.5 . When z is not omitted, the simulation yields the green sampling distribution, centered around the true value. The average value across all simulations is 0.4998. When z is omitted, the simulation yields the red sampling distribution, centered around 0.5895. It’s biased upward from the true value of .5 by .0895. Moreover, the variance of the biased sampling distribution is much smaller than the true variance around \beta . This compromises the ability to perform any meaningful inferences about the true parameter. End2Bias Analysis The bias in \beta can be derived analytically. Consider that in model 1 (presented above), x and z are related by:

(3)\thinspace z = d+hx+\epsilon_2

Substituting z in equation 1 with equation 3 and re-ordering:

 Y = a+\beta x+c (d+hx+\epsilon_2) + \epsilon_1

 (4)\thinspace Y = (a+cd)+(\beta+ch) x + (\epsilon_1+c\epsilon_2)

When omitting variable z , it is actually equation 4 that is estimated. It can be seen that \beta is biased by the quantity ch . In this case, since x and z are positively correlated by construction and their slope coefficients are positive, the bias will be positive. According to the parameters of the simulation, the “true” bias should be  ch=.09. Here is the distribution of the bias, it is centered around .0895, very close to the true bias value. bias2 The derivation above also lets us determine the direction of bias from knowing the correlation of x and z as well as the sign of c (the true partial effect of z on y ). If both are the same sign, then the estimate of \beta will be biased upward. If the signs differ, then the estimate of \beta will be biased downward. Conclusion The case above was pretty general, but has particular applications. For example, if we believe that an individual’s income is a function of years of education and year of work experience, then omitting one variable will bias the slope estimate of the other.

Motor Vehicle Collision Density in NYC

DensityByBorough

In a previous post, I visualized crime density in Boston using R’s ggmap package. In this post, I use ggmap to visualize vehicle accidents in New York City. R code and data are included in this post.

The data comes from NYC Open Data. My data cut ranges from 2012 to 2015. The data tracks the type of vehicle, the name of the street on which the accident occurs, as well as the longitude and latitude coordinates of the accident. Both coordinates are saved as a single character variable called “LOCATION”.

Below, I load the ggmap and gridExtra packages, load the data, drop all accidents without location coordinates, and parse the LOCATION variable to get the longitude and latitude coordinates. I also parse the date variable to create a year variable and use that variable to create two data sets: one with all vehicle accidents in 2013 and another with all vehicle accidents in 2014

 
library(ggmap)
library(gridExtra)

d=read.csv('.../NYPD_Motor_Vehicle_Collisions.csv')

d_clean=d[which(regexpr(',',d$LOCATION)!=-1),]

comm=regexpr(',',d_clean$LOCATION)
d_clean$loc=as.character(d_clean$LOCATION)
d_clean$lat=as.numeric(substr(d_clean$loc,2,comm-1))
d_clean$long=as.numeric(substr(d_clean$loc,comm+1,nchar(d_clean$loc)-1))
d_clean$year=substr(d_clean$DATE,7,10)

d_2013=d_clean[which(d_clean$year=='2013'),c('long','lat')]
d_2014=d_clean[which(d_clean$year=='2014'),c('long','lat')]

Continue reading Motor Vehicle Collision Density in NYC