-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathultrapollution_PM25.R
98 lines (77 loc) · 5.67 KB
/
ultrapollution_PM25.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
# This script performs the econometric analyses for the ultrapollution project
# Loading libraries
library(stargazer)
library(sandwich)
library(lmtest)
# Specifying directories for data + results
direc <- paste('C:/Users/', username, '/Documents/Data/ultrapollution/ultra_data/', sep = '')
direc2 <- paste('C:/Users/', username, '/Documents/Data/ultrapollution/results/', sep = '')
# Reading in the data set
data <- read.csv(paste(direc, 'ultradata.csv', sep = ''))
# See which event types have at least 100 observations
event_types <- c()
for (d in unique(data$RACE_Distance)) {if (dim(data[which(data$RACE_Distance == d),][1])>100) {event_types <- c(event_types,d)}}
# Run OLS regressions
mod50k <- lm(log(Seconds+.001) ~ PM2.5 + factor(Gender) + Temperature*Humidity + Precipitation
+ WindSpeed + factor(FIPS_Race) + factor(RACE_Month)*factor(RACE_Year)
+ Total_Races + RACE_Finisher_Count + In_State
+ Travel_Distance + Ability + Population + Poverty_Rate + Unemployment_Rate + Income
+ Education_High_School + Education_Some_College + Education_Associates
+ Education_Bachelors + Education_Graduate + Altitude_Home + PM2.5_Home,
data = data[which(data$RACE_Distance == '50 KM'),])
mod50m <- lm(log(Seconds+.001) ~ PM2.5 + factor(Gender) + Temperature*Humidity + Precipitation
+ WindSpeed + factor(FIPS_Race) + factor(RACE_Month)*factor(RACE_Year)
+ Total_Races + RACE_Finisher_Count + In_State
+ Travel_Distance + Ability + Population + Poverty_Rate + Unemployment_Rate + Income
+ Education_High_School + Education_Some_College + Education_Associates
+ Education_Bachelors + Education_Graduate + Altitude_Home + PM2.5_Home,
data = data[which(data$RACE_Distance == '50 Miles'),])
mod100k <- lm(log(Seconds+.001) ~ PM2.5 + factor(Gender) + Temperature*Humidity + Precipitation
+ WindSpeed + factor(FIPS_Race) + factor(RACE_Month)*factor(RACE_Year)
+ Total_Races + RACE_Finisher_Count + In_State
+ Travel_Distance + Ability + Population + Poverty_Rate + Unemployment_Rate + Income
+ Education_High_School + Education_Some_College + Education_Associates
+ Education_Bachelors + Education_Graduate + Altitude_Home + PM2.5_Home,
data = data[which(data$RACE_Distance == '100 KM'),])
mod100m <- lm(log(Seconds+.001) ~ PM2.5 + factor(Gender) + Temperature*Humidity + Precipitation
+ WindSpeed + factor(FIPS_Race) + factor(RACE_Month)*factor(RACE_Year)
+ Total_Races + RACE_Finisher_Count + In_State
+ Travel_Distance + Ability + Population + Poverty_Rate + Unemployment_Rate + Income
+ Education_High_School + Education_Some_College + Education_Associates
+ Education_Bachelors + Education_Graduate + Altitude_Home + PM2.5_Home,
data = data[which(data$RACE_Distance == '100 Miles'),])
mod6h <- lm(log(Distance+.001) ~ PM2.5 + factor(Gender) + Temperature*Humidity + Precipitation
+ WindSpeed + factor(FIPS_Race) + factor(RACE_Month)*factor(RACE_Year)
+ Total_Races + RACE_Finisher_Count + In_State
+ Travel_Distance + Ability + Population + Poverty_Rate + Unemployment_Rate + Income
+ Education_High_School + Education_Some_College + Education_Associates
+ Education_Bachelors + Education_Graduate + Altitude_Home + PM2.5_Home,
data = data[which(data$RACE_Distance == '6 Hours'),])
mod12h <- lm(log(Distance+.001) ~ PM2.5 + factor(Gender) + Temperature*Humidity + Precipitation
+ WindSpeed + factor(FIPS_Race) + factor(RACE_Month)*factor(RACE_Year)
+ Total_Races + RACE_Finisher_Count + In_State
+ Travel_Distance + Ability + Population + Poverty_Rate + Unemployment_Rate + Income
+ Education_High_School + Education_Some_College + Education_Associates
+ Education_Bachelors + Education_Graduate + Altitude_Home + PM2.5_Home,
data = data[which(data$RACE_Distance == '12 Hours'),])
mod24h <- lm(log(Distance+.001) ~ PM2.5 + factor(Gender) + Temperature*Humidity + Precipitation
+ WindSpeed + factor(FIPS_Race) + factor(RACE_Month)*factor(RACE_Year)
+ Total_Races + RACE_Finisher_Count + In_State
+ Travel_Distance + Ability + Population + Poverty_Rate + Unemployment_Rate + Income
+ Education_High_School + Education_Some_College + Education_Associates
+ Education_Bachelors + Education_Graduate + Altitude_Home + PM2.5_Home,
data = data[which(data$RACE_Distance == '24 Hours'),])
xmod50k <- coeftest(mod50k, vcov = vcovCL, cluster = ~RACE_Name)
xmod50m <- coeftest(mod50m, vcov = vcovCL, cluster = ~RACE_Name)
xmod100k <- coeftest(mod100k, vcov = vcovCL, cluster = ~RACE_Name)
xmod100m <- coeftest(mod100m, vcov = vcovCL, cluster = ~RACE_Name)
xmod6h <- coeftest(mod6h, vcov = vcovCL, cluster = ~RACE_Name)
xmod12h <- coeftest(mod12h, vcov = vcovCL, cluster = ~RACE_Name)
xmod24h <- coeftest(mod24h, vcov = vcovCL, cluster = ~RACE_Name)
# Viewing results
write.csv(stargazer(xmod50k, xmod50m, xmod100k, xmod100m, xmod6h, xmod12h, xmod24h,
omit = c('FIPS_Race', 'RACE_Name', 'RACE_Month', 'RACE_Year')),
paste(direc2, 'PM_results_tex.txt', sep = ''), row.names = FALSE)
write.csv(stargazer(xmod50k, xmod50m, xmod100k, xmod100m, xmod6h, xmod12h, xmod24h, type = 'text',
omit = c('FIPS_Race', 'RACE_Name', 'RACE_Month', 'RACE_Year')),
paste(direc2, 'PM_results.txt', sep = ''), row.names = FALSE)