McCandless_MN4111_Wksht1
Conor McCandless
2024-04-08
Load libraries
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(wooldridge)
library(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerlibrary(lmtest)## Warning: package 'lmtest' was built under R version 4.3.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.3.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(sandwich)## Warning: package 'sandwich' was built under R version 4.3.3library(purrr)
library(ggplot2)
library(zoo)
options(scipen=999)Question 1 (i) Motivation is the unobserved variable that is confounding because it is correlated with pros and cons and on the selection. Because less voted Marines are most likely to join in order to get the training making it a negative treatment bias.
- Because the original controls resulted in a negative bias, the addition of the control for the AMTB brought us closer its true value.
Question 2 (i)
data("bwght")
summary("bwght")## Length Class Mode
## 1 character characterwork<-bwght
hist(work$bwght)
hist(work$cigs)
(ii)
p1<-ggplot(data=work, aes(x=cigs, y=bwght))+geom_point()
p1
p2<-ggplot(data = work, aes(x = cigs, y = bwght)) +
stat_summary_bin(fun='mean', bins=20, size=4, geom='point') +
stat_smooth(method = "lm", se= FALSE)
p2## `geom_smooth()` using formula = 'y ~ x'
The more cigs you smoke, the lower the birth weight.
Family income goes up, birth weight goes up. Education, access to care, general health nutrition.
Family income goes up, cigerette use goes down
mod1<-lm(bwght~cigs, data=work)
summary(mod1)##
## Call:
## lm(formula = bwght ~ cigs, data = work)
##
## Residuals:
## Min 1Q Median 3Q Max
## -96.772 -11.772 0.297 13.228 151.228
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 119.77190 0.57234 209.267 < 0.0000000000000002 ***
## cigs -0.51377 0.09049 -5.678 0.0000000166 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.13 on 1386 degrees of freedom
## Multiple R-squared: 0.02273, Adjusted R-squared: 0.02202
## F-statistic: 32.24 on 1 and 1386 DF, p-value: 0.00000001662se_robust <- as.vector(coeftest(mod1, vcov. = vcovHC(mod1, type = "HC1"))[,"Std. Error"])
stargazer(mod1, type="text")##
## ===============================================
## Dependent variable:
## ---------------------------
## bwght
## -----------------------------------------------
## cigs -0.514***
## (0.090)
##
## Constant 119.772***
## (0.572)
##
## -----------------------------------------------
## Observations 1,388
## R2 0.023
## Adjusted R2 0.022
## Residual Std. Error 20.129 (df = 1386)
## F Statistic 32.235*** (df = 1; 1386)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01 se=list(se_robust)
mod2<-lm(bwght~cigs, data=work)
summary(mod2)##
## Call:
## lm(formula = bwght ~ cigs, data = work)
##
## Residuals:
## Min 1Q Median 3Q Max
## -96.772 -11.772 0.297 13.228 151.228
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 119.77190 0.57234 209.267 < 0.0000000000000002 ***
## cigs -0.51377 0.09049 -5.678 0.0000000166 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.13 on 1386 degrees of freedom
## Multiple R-squared: 0.02273, Adjusted R-squared: 0.02202
## F-statistic: 32.24 on 1 and 1386 DF, p-value: 0.00000001662se_robust2 <- as.vector(coeftest(mod2, vcov. = vcovHC(mod2, type = "HC1"))[,"Std. Error"])
stargazer(mod2, type="text")##
## ===============================================
## Dependent variable:
## ---------------------------
## bwght
## -----------------------------------------------
## cigs -0.514***
## (0.090)
##
## Constant 119.772***
## (0.572)
##
## -----------------------------------------------
## Observations 1,388
## R2 0.023
## Adjusted R2 0.022
## Residual Std. Error 20.129 (df = 1386)
## F Statistic 32.235*** (df = 1; 1386)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01se=list(se_robust)
stargazer(mod1,mod2, type= "text", se=list(se_robust, se_robust2))##
## ============================================================
## Dependent variable:
## ----------------------------
## bwght
## (1) (2)
## ------------------------------------------------------------
## cigs -0.514*** -0.514***
## (0.088) (0.088)
##
## Constant 119.772*** 119.772***
## (0.575) (0.575)
##
## ------------------------------------------------------------
## Observations 1,388 1,388
## R2 0.023 0.023
## Adjusted R2 0.022 0.022
## Residual Std. Error (df = 1386) 20.129 20.129
## F Statistic (df = 1; 1386) 32.235*** 32.235***
## ============================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Cigs in first model is -.514 and -.463 when controlled for family income
Number of classes becuase it is a downstream effect and correlated with being eligible.
#Estimate regressions using lm()
mod0 = lm(bwght ~ cigs , data=work)
mod1 = lm(bwght ~ cigs+lfaminc , data=work)
mod2 = lm(bwght ~ cigs+lfaminc+fatheduc+motheduc , data=work)
mod3 = lm(bwght ~ cigs+lfaminc+fatheduc+motheduc+male+parity+white, data=work)
# Save cluster and het. robust standard errors
se0 = as.vector(coeftest(mod0,vcov = vcovCL, type="HC1")[,"Std. Error"]) # Here we just pull the existing standard errors from the regression object
se1 = as.vector(coeftest(mod1,vcov = vcovCL, type="HC1")[,"Std. Error"]) # Cluster SEs at "foreign" level
se2 = as.vector(coeftest(mod2,vcov = vcovCL, type="HC1")[,"Std. Error"])
se3 = as.vector(coeftest(mod3,vcov = vcovCL, type="HC1")[,"Std. Error"])
## Calculate outcome means
# Using map/pluck
outcome_means = list(mod0, mod1, mod2, mod3) %>% map(pluck,"fitted.values") %>% map_dbl(mean)
# Using lapply
#outcome_means = lapply(list(mod0, mod1, mod2, mod3), function(reg) mean(reg$model[,1]))
#outcome_means = unlist(outcome_means)
# Make table
stargazer(mod0, mod1, mod2, mod3, type = "text",
title = "Effects of cigarettes on birthweight",
omit.stat = c("f", "ser", "adj.rsq"),
add.lines = list(c("Outcome mean", round(outcome_means, 2))),
dep.var.labels.include = FALSE,
dep.var.caption = "",
omit = c("Constant"),
covariate.labels = c("Cigs Smoked per day", "Ln(family income)", "Father Education",
"Mother Education", "Male infant", "Parity", "White"))##
## Effects of cigarettes on birthweight
## ===========================================================
## (1) (2) (3) (4)
## -----------------------------------------------------------
## Cigs Smoked per day -0.514*** -0.467*** -0.589*** -0.597***
## (0.090) (0.091) (0.111) (0.110)
##
## Ln(family income) 1.851*** 1.458 1.221
## (0.595) (0.901) (0.912)
##
## Father Education 0.500* 0.415
## (0.281) (0.279)
##
## Mother Education -0.445 -0.337
## (0.319) (0.318)
##
## Male infant 3.825***
## (1.143)
##
## Parity 1.918***
## (0.655)
##
## White 4.637***
## (1.612)
##
## -----------------------------------------------------------
## Outcome mean 118.7 118.7 119.53 119.53
## Observations 1,388 1,388 1,191 1,191
## R2 0.023 0.030 0.033 0.054
## ===========================================================
## Note: *p<0.1; **p<0.05; ***p<0.01 se = list(se0, se1, se2, se3)
notes = c("Standard errors in parentheses",
"Outcome in all specifications is birthweight in ounces")
notes.align = "l"