Extra Credit Appendix

library(sf)
library(tmap)
library(tigris)
options(tigris_class = "sf")
ny_county <- counties("NY", cb = TRUE)

Using FIPS code '36' for state 'NY'

NYmap <- geo_join(ny_county, NY_Climate2, by = "NAME")

st_crs<- : replacing crs does not reproject data; use st_transform for that

tm_shape(NYmap) + tm_polygons("Pop2", palette = "Greens") + tm_text("NAME", size = "AREA")

library(plyr)
library(dplyr)
library(knitr)
library(DT)
library(xtable)
nyef5 <- NY_Climate2 %>% 
  arrange(desc(Pop2)) %>%
  top_n(n = 10)

Selecting by Pop2

kable(nyef5)

|NAME        |Pop2 |
|:-----------|:----|
|Bronx       |9    |
|Erie        |9    |
|Kings       |9    |
|Nassau      |9    |
|Queens      |9    |
|Suffolk     |9    |
|Westchester |9    |
|Monroe      |8    |
|New York    |5    |
|Onondaga    |5    |
|Richmond    |5    |

library(texreg)
nyef1 <- lm(Income ~ Pop, data = NY_Efficiency)
htmlreg(nyef1)

Statistical models

	Model 1
(Intercept)	23.32***
	(3.47)
Pop2	0.05
	(6.18)
Pop3	3.50
	(6.27)
Pop4	7.66
	(5.15)
Pop5	6.73
	(4.89)
Pop8	3.68
	(6.74)
Pop9	29.36***
	(3.56)
R2	0.07
Adj. R2	0.07
Num. obs.	2299
RMSE	34.68
p < 0.001, p < 0.01, p < 0.05

with(NY_Efficiency, plot(Pop, Income))
abline(nyef1)

only using the first two of 7 regression coefficients

library(merTools)
library(lme4)
library(ggplot2)
library(Zelig)
NY_Efficiency$Pop <- as.factor(NY_Efficiency$Pop)
z.nye2 <- zelig(Unemploy ~ Pop, model = "poisson", data = NY_Efficiency, cite = F)
z.nye2$setx(Pop = "1")
z.nye2$setx1(Pop = "9")
z.nye2$sim()
z.nye2$graph()

fd2 <- z.nye2$get_qi(xvalue="x1", qi="fd")
summary(fd2)

       V1        
 Min.   :-3.121  
 1st Qu.:-2.372  
 Median :-2.192  
 Mean   :-2.200  
 3rd Qu.:-2.029  
 Max.   :-1.408  

library(merTools)
library(lme4)
library(ggplot2)
library(Zelig)
library(texreg)
NY_Efficiency$Pop <- as.factor(NY_Efficiency$Pop)
z.nye3 <- zelig(Edu ~ Pop, model = "poisson", data = NY_Efficiency, cite = F)
z.nye3$setx(Pop = "1")
z.nye3$setx1(Pop = "9")
z.nye3$sim()
htmlreg(z.nye3, doctype = FALSE)

Statistical models

	Model 1
(Intercept)	3.11***
	(0.02)
Pop2	0.12**
	(0.04)
Pop3	0.21***
	(0.04)
Pop4	0.53***
	(0.03)
Pop5	0.35***
	(0.03)
Pop8	0.53***
	(0.03)
Pop9	0.57***
	(0.02)
AIC	23037.31
BIC	23077.51
Log Likelihood	-11511.65
Deviance	10558.68
Num. obs.	2307
p < 0.001, p < 0.01, p < 0.05

fd3 <- z.nye3$get_qi(xvalue="x1", qi="fd")
summary(fd3)

   V1       

Min. :15.73
1st Qu.:16.80
Median :17.14
Mean :17.12
3rd Qu.:17.45
Max. :18.56

library(plyr)
library(knitr)
library(DT)
library(xtable)
nyef4 <- NY_Efficiency %>% 
  dplyr::group_by(County)%>%
  dplyr::summarize(p_county = n()) %>%
  arrange(desc(p_county)) %>%
  top_n(n = 10)

Selecting by p_county

kable(nyef4)

|County      | p_county|
|:-----------|--------:|
|NEW YORK    |      502|
|KINGS       |      418|
|BRONX       |      304|
|QUEENS      |      271|
|WESTCHESTER |      192|
|SUFFOLK     |       87|
|RICHMOND    |       77|
|NASSAU      |       73|
|ALBANY      |       57|
|ERIE        |       42|

library(ggplot2)
library(texreg)
NY_Efficiency$TPCost <- as.numeric(gsub('[$,]', '', NY_Efficiency$TPCost))
nyef9 <- lm(TPCost ~ Pop, data = NY_Efficiency)
htmlreg(nyef9)

Statistical models

	Model 1
(Intercept)	472356.56
	(331769.37)
Pop2	478437.35
	(607040.58)
Pop3	599028.19
	(616640.06)
Pop4	1029039.31*
	(503284.92)
Pop5	513437.04
	(477253.08)
Pop8	1640769.71*
	(663538.73)
Pop9	715484.87*
	(341121.69)
R2	0.00
Adj. R2	0.00
Num. obs.	2307
RMSE	3447848.38
p < 0.001, p < 0.01, p < 0.05

ggplot(NY_Efficiency, aes(x = Pop, y = TPCost)) + 
  geom_point()

library(Zelig)
NY_Energy_Efficiency$FYSavings <- as.numeric(gsub('[$,]', '', NY_Energy_Efficiency$FYSavings))
NY_Energy_Efficiency$TPCost <- as.numeric(gsub('[$,]', '', NY_Energy_Efficiency$TPCost))
z.nyef6 <- zelig(FYSavings ~ TPCost, model = "poisson", data = NY_Energy_Efficiency, cite = F) 
summary(z.nyef6)

Model: 

Call:
z5$zelig(formula = FYSavings ~ TPCost, data = NY_Energy_Efficiency)

Deviance Residuals: 
    Min       1Q   Median       3Q  
-3167.1   -318.6   -196.8    -15.0  
    Max  
 5258.4  

Coefficients:
             Estimate Std. Error z value
(Intercept) 1.124e+01  7.392e-05  152076
TPCost      7.464e-08  3.917e-12   19056
            Pr(>|z|)
(Intercept)   <2e-16
TPCost        <2e-16

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 528356365  on 2306  degrees of freedom
Residual deviance: 377712188  on 2305  degrees of freedom
  (6 observations deleted due to missingness)
AIC: 377737792

Number of Fisher Scoring iterations: 7

Next step: Use 'setx' method

z.nyef6$setx()
z.nyef6$sim()
plot(z.nyef6)

library(texreg)
NY_Energy_Efficiency$TPCost <- as.numeric(gsub('[$,]', '', NY_Energy_Efficiency$TPCost))
nyef8 <- lm(EEReduction ~ TPCost, data = NY_Energy_Efficiency)
htmlreg(nyef8)

Statistical models

	Model 1
(Intercept)	495711.92***
	(33059.32)
TPCost	0.13***
	(0.01)
R2	0.08
Adj. R2	0.08
Num. obs.	2307
RMSE	1504240.08
p < 0.001, p < 0.01, p < 0.05

with(NY_Energy_Efficiency, plot(TPCost, EEReduction))
abline(nyef8)