CS112 Assignment 3, Fall 2020
Ha Tran Nguyen Phuong
11/01/2020
A few important notes
Option 1 for submitting your assignment: This method is actually preferred. This is an RMarkdown document. Did you know you can open this document in RStudio, edit it by adding your answers and code, and then knit it to a pdf? To submit your answers to this assignment, simply knit this file as a pdf and submit it as a pdf on Forum. All of your code must be included in the resulting pdf file, i.e., don’t set echo = FALSE in any of your code chunks. To learn more about RMarkdown, watch the videos from session 1 and session 2 of the CS112B optional class. This is also a cheat sheet for using Rmarkdown. If you have questions about RMarkdown, please post them on Perusall. Try knitting this document in your RStudio. You should be able to get a pdf file. At any step, you can try knitting the document and recreate a pdf. If you get error, you might have incomplete code.
Option 2 for submitting your assignment: If you are not comfortable with RMarkdown, you can also choose the Google Doc version of this assignment, make a copy of it and edit the Google doc (include your code, figures, results, and explanations) and at the end download your Google Doc as a pdf and submit the pdf file.
Note: Either way (if you use Rmd and knit as pdf OR if you use Google Doc and download as pdf) you should make sure you put your name on top of the document.
Note: The first time you run this document you may get error that some packages don’t exist. If you don’t have the packages listed on top of this document, install them first and you won’t get those errors.
Note: If you work with others or get help from others on the assignment, please provide the names of your partners at the top of this document. Working together is fine, but all must do and understand their own work and have to mention the collaboration on top of this document.
Note: Don’t change seed in the document. The function set.seed() has already been set at the beginning of this document to 928. Changing the see again to a different number will make your results not replicable.
QUESTION 1: Data Generating Example
The following code, creates a toy dataset with a treatment variable, \(D\), an outcome variable, \(Y\), and other variables \(V_1\) to \(V_4\).
n = 1000
## Generating a random data set here
#Syntax for the normal distribution here is rnorm(sample size, mean, SD)
V1 = rnorm(n, 45, 10)
#getting a binary variable
V2 = sample(c(1,0),
replace = TRUE,
size = n,
prob = c(.4,.6))
V3 = rnorm(n, V1/10, 1)
V4 = rnorm(n, 0, 1)
D = as.numeric(pnorm(rnorm(n, .01*V1 + .8*V2 + 0.3*V3 + V4, 1), .45 + .32 + .3*4.5, 1) > .5)
Y = rnorm(n, .8*D - 0.45*V2 - .4*V3 + 2, .2)
# combining everything in a data frame
df = data.frame(V1, V2, V3, V4, D, Y)STEP 1
From the variables \(V_1\), \(V_2\), \(V_3\), and \(V_4\), which one(s) are not confounding variable(s) (covariates that causes confounding)? Remember, a rule of thumb (although not a perfect rule) is that the variable is correlated with both the treatment variable and the outcome variable. Explain!
V2 and V3 are confounding variable because they affect both treatment variable D and outcome variable Y (D is affected by 0.8V2 and 0.3 V3, while Y is affected by -0.45V2 and 0.4V3)
STEP 2
Can you figure out the true treatment effect by looking at the data generating process above? Yes, it’s 0.8 i.e. treatment increases Y by 0.8 units. This is because the systematic component of Y includes 0.8D. #### STEP 3
Plot the outcome variable against the treatment variable. Make sure you label your axes. Do you see a trend?
# Your code here
plot(df$D, df$Y, xlab = "Treatment", ylab ="Outcome")
abline(lm(D ~ Y, data=df))
From the trendline, we can see a slight positive relationship i.e. there is a positive treatment effect. #### STEP 4
Are the variables \(V_1\), \(V_2\), \(V_3\), and \(V_4\) balanced across the treatment and control groups? You can use any R function from any package to check this (for instance, you can check the cobalt package). Make sure you check all variables.
Note: This is optional but you can use the gridExtra package and its grid.arrange() function to put all the 4 graphs in one 2 x 2 graph. Read more about the package and how to use it here: https://cran.r-project.org/web/packages/egg/vignettes/Ecosystem.html. Set nrow = 2.
# Your code here
w.out <- weightit(D~V1 + V2 + V3 + V4, data = df)
plot1 = bal.plot(w.out, "V1", which = "unadjusted")
plot2 = bal.plot(w.out, "V2", which = "unadjusted")
plot3 = bal.plot(w.out, "V3", which = "unadjusted")
plot4 = bal.plot(w.out, "V4", which = "unadjusted")
grid.arrange(plot1, plot2, plot3, plot4, nrow = 2)
MatchBalance(D~V1 + V2 + V3 + V4, data = df)##
## ***** (V1) V1 *****
## before matching:
## mean treatment........ 47.145
## mean control.......... 42.457
## std mean diff......... 46.267
##
## mean raw eQQ diff..... 4.7386
## med raw eQQ diff..... 4.6719
## max raw eQQ diff..... 8.2566
##
## mean eCDF diff........ 0.13197
## med eCDF diff........ 0.14894
## max eCDF diff........ 0.22444
##
## var ratio (Tr/Co)..... 1.1115
## T-test p-value........ 1.4033e-13
## KS Bootstrap p-value.. < 2.22e-16
## KS Naive p-value...... 2.3252e-11
## KS Statistic.......... 0.22444
##
##
## ***** (V2) V2 *****
## before matching:
## mean treatment........ 0.52749
## mean control.......... 0.28684
## std mean diff......... 48.155
##
## mean raw eQQ diff..... 0.24236
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.12033
## med eCDF diff........ 0.12033
## max eCDF diff........ 0.24066
##
## var ratio (Tr/Co)..... 1.2185
## T-test p-value........ 4.4409e-15
##
##
## ***** (V3) V3 *****
## before matching:
## mean treatment........ 4.919
## mean control.......... 4.0853
## std mean diff......... 61.744
##
## mean raw eQQ diff..... 0.84184
## med raw eQQ diff..... 0.82372
## max raw eQQ diff..... 1.7613
##
## mean eCDF diff........ 0.16484
## med eCDF diff........ 0.17683
## max eCDF diff........ 0.26201
##
## var ratio (Tr/Co)..... 0.954
## T-test p-value........ < 2.22e-16
## KS Bootstrap p-value.. < 2.22e-16
## KS Naive p-value...... 2.5535e-15
## KS Statistic.......... 0.26201
##
##
## ***** (V4) V4 *****
## before matching:
## mean treatment........ 0.56565
## mean control.......... -0.49923
## std mean diff......... 127.27
##
## mean raw eQQ diff..... 1.0694
## med raw eQQ diff..... 1.0629
## max raw eQQ diff..... 1.8961
##
## mean eCDF diff........ 0.31869
## med eCDF diff........ 0.35388
## max eCDF diff........ 0.50036
##
## var ratio (Tr/Co)..... 1.0242
## T-test p-value........ < 2.22e-16
## KS Bootstrap p-value.. < 2.22e-16
## KS Naive p-value...... 0
## KS Statistic.......... 0.50036
##
##
## Before Matching Minimum p.value: < 2.22e-16
## Variable Name(s): V1 V3 V4 Number(s): 1 3 4From bal.plot, we can see the distribution is not very balanced for all 4 values - they have roughly the same distribution, but different means. Using MatchBalance, we can also see minimum p-value before matching. We see the minimum p-value is less than 2.22e^-16 for V1, V3 and V4, which would mean the difference between control and treatment group (both in mean and distribution) is highly statistically significant. Thus the values are definitely not balanced. #### STEP 5
Write code that would simply calculate the Prima Facie treatment effect in the data above. What’s the Prima Facie treatment effect? Note that the outcome variable is binary.
# Your code here
mean(df[df$D==1,]$Y) - mean(df[df$D==0,]$Y)## [1] 0.3423726The Prima Facie Treatment effect is 0.3423. This is the effect when we have the naive assumption that the actual average treatment effect is equal to the difference between average outcome forthe treatment group and the control group. However, as we have the formula and the true treatment effect, we know this is not the true treatment effect because we have confounders (V2 and V3) and our data is not well-matched.
STEP 6
Explain why the Prima Facie effect is not the true average causal effect from Step 2.
We have 2 confounders, V2 and V3, which affect both the treatment variable and the outcome variable. Thus, if you have a particular value of V2, you are more likely to be in treatment and also more likely to have a better outcome Y, for instance. So it is actually not the treatment variable D that makes you have a higher Y, but the variable V2 that causes this increase. Thus, the Prima Facie effect is different from the true average causal effect.
STEP 7
We can use matching to create a better balance of the covariates. Use a propensity score model that includes all the variables \(V_1\), \(V_2\), \(V_3\), and \(V_4\).
# Your code here
#creating prospensity score
ps <- glm(D ~ V1+V2+V3+V4, family = 'binomial', data = df)
#specifying X and Y variable
X <- ps$fitted
Y <- df$Y
Tr <- df$D
#finding the treatment effect using Matching
mout2 <- Match(Y = Y, Tr = Tr, X = X)
summary(mout2)##
## Estimate... 0.81326
## AI SE...... 0.16221
## T-stat..... 5.0137
## p.val...... 5.3395e-07
##
## Original number of observations.............. 1000
## Original number of treated obs............... 491
## Matched number of observations............... 491
## Matched number of observations (unweighted). 624STEP 8
Check the balance of the covariates. Do you see any improvements after matching?
# Your code here
#checking if they are balanced
mout <- Match(Tr = Tr, X = X)
mb <- MatchBalance(D ~ V1 + V2 + V3 + V4, data = df, match.out = mout)##
## ***** (V1) V1 *****
## Before Matching After Matching
## mean treatment........ 47.145 47.145
## mean control.......... 42.457 46.073
## std mean diff......... 46.267 10.572
##
## mean raw eQQ diff..... 4.7386 2.4364
## med raw eQQ diff..... 4.6719 1.7626
## max raw eQQ diff..... 8.2566 9.4025
##
## mean eCDF diff........ 0.13197 0.065311
## med eCDF diff........ 0.14894 0.054487
## max eCDF diff........ 0.22444 0.17147
##
## var ratio (Tr/Co)..... 1.1115 1.6438
## T-test p-value........ 1.4033e-13 0.060663
## KS Bootstrap p-value.. < 2.22e-16 < 2.22e-16
## KS Naive p-value...... 2.3252e-11 2.1513e-08
## KS Statistic.......... 0.22444 0.17147
##
##
## ***** (V2) V2 *****
## Before Matching After Matching
## mean treatment........ 0.52749 0.52749
## mean control.......... 0.28684 0.58847
## std mean diff......... 48.155 -12.201
##
## mean raw eQQ diff..... 0.24236 0.049679
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.12033 0.02484
## med eCDF diff........ 0.12033 0.02484
## max eCDF diff........ 0.24066 0.049679
##
## var ratio (Tr/Co)..... 1.2185 1.0292
## T-test p-value........ 4.4409e-15 0.036802
##
##
## ***** (V3) V3 *****
## Before Matching After Matching
## mean treatment........ 4.919 4.919
## mean control.......... 4.0853 4.8504
## std mean diff......... 61.744 5.0783
##
## mean raw eQQ diff..... 0.84184 0.29996
## med raw eQQ diff..... 0.82372 0.27618
## max raw eQQ diff..... 1.7613 1.27
##
## mean eCDF diff........ 0.16484 0.058721
## med eCDF diff........ 0.17683 0.044872
## max eCDF diff........ 0.26201 0.14103
##
## var ratio (Tr/Co)..... 0.954 1.0788
## T-test p-value........ < 2.22e-16 0.36665
## KS Bootstrap p-value.. < 2.22e-16 < 2.22e-16
## KS Naive p-value...... 2.5535e-15 8.1531e-06
## KS Statistic.......... 0.26201 0.14103
##
##
## ***** (V4) V4 *****
## Before Matching After Matching
## mean treatment........ 0.56565 0.56565
## mean control.......... -0.49923 0.38581
## std mean diff......... 127.27 21.493
##
## mean raw eQQ diff..... 1.0694 0.15808
## med raw eQQ diff..... 1.0629 0.06365
## max raw eQQ diff..... 1.8961 1.8961
##
## mean eCDF diff........ 0.31869 0.04643
## med eCDF diff........ 0.35388 0.019231
## max eCDF diff........ 0.50036 0.21795
##
## var ratio (Tr/Co)..... 1.0242 1.4111
## T-test p-value........ < 2.22e-16 6.192e-06
## KS Bootstrap p-value.. < 2.22e-16 < 2.22e-16
## KS Naive p-value...... < 2.22e-16 2.6801e-13
## KS Statistic.......... 0.50036 0.21795
##
##
## Before Matching Minimum p.value: < 2.22e-16
## Variable Name(s): V1 V3 V4 Number(s): 1 3 4
##
## After Matching Minimum p.value: < 2.22e-16
## Variable Name(s): V1 V3 V4 Number(s): 1 3 4While after matching the minimum p.value stays the same for all the 3 variables, they are dependent on KS bootstrap p-value. If we look at the standard mean difference, there is reduction in all variables i.e. the value gets closer to 0. The T-test value also increase for all 4 variables, which indicate that they are better matched.
STEP 9
What is the treatment effect after matching? Is this surprising given your answer to Step 2. Is the treatment effect found in this Step closer to the treatment effect in Step 2 than the treatment effect before matching?
summary(mout2)##
## Estimate... 0.81326
## AI SE...... 0.16221
## T-stat..... 5.0137
## p.val...... 5.3395e-07
##
## Original number of observations.............. 1000
## Original number of treated obs............... 491
## Matched number of observations............... 491
## Matched number of observations (unweighted). 624The treatment effect after matching is roughly 0.81, which is very close to the actual treatment effect and a lot better than the treatment effect before matching (which was only 0.31). Thus, even though the matching isn’t the best (because after matching minimum p-value for the 3 variables were still very high), it is still significantly better than without matching.
QUESTION 2: Daughters
Read Section 5 (which is only about 1 page long!) of Iacus, King, Porro (2011), Multivariate Matching Methods That Are Monotonic Imbalance Bounding, JASA, V 106, N. 493, available here: https://gking.harvard.edu/files/gking/files/cem_jasa.pdf. Don’t worry about the “CEM” elements. Focus on the “daughters” case study.
Data for this case study is available in “doughters” below.
daughters = read.csv("http://bit.ly/daughters_data") %>%
clean_names()STEP 1
Before doing any matching, run a regression, with the outcome variable being nowtot, the treatment variable being hasgirls, and the independent vars mentioned below: - dems, - repubs, - christian, - age, - srvlng, - demvote
Show the regression specification. Use the regression to estimate a treatment effect and confidence interval. Check the balance of this not-matched data set using any method of your choice (balance tables, balance plots, love plots, etc).
# Your code here
#creating linear regression model
lm1 <- lm(nowtot ~ hasgirls + dems + repubs + christian + age + srvlng + demvote, data = daughters)
summary(lm1)##
## Call:
## lm(formula = nowtot ~ hasgirls + dems + repubs + christian +
## age + srvlng + demvote, data = daughters)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56.028 -10.322 -1.517 11.208 69.642
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.6991 18.6306 2.077 0.038390 *
## hasgirls -0.4523 1.9036 -0.238 0.812322
## dems -8.1022 17.5861 -0.461 0.645238
## repubs -55.1069 17.6340 -3.125 0.001901 **
## christian -13.3961 3.7218 -3.599 0.000357 ***
## age 0.1260 0.1117 1.128 0.259938
## srvlng -0.2251 0.1355 -1.662 0.097349 .
## demvote 87.5501 8.4847 10.319 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.19 on 422 degrees of freedom
## Multiple R-squared: 0.7821, Adjusted R-squared: 0.7784
## F-statistic: 216.3 on 7 and 422 DF, p-value: < 2.2e-16#treatment effect estimate
lm1$coefficients['hasgirls']## hasgirls
## -0.4522678#finding confidence interval
confint(lm1,'hasgirls',level=0.95)## 2.5 % 97.5 %
## hasgirls -4.19406 3.289525#checking balance
w.out_daughter <- weightit(hasgirls ~ dems + repubs + christian + age + srvlng + demvote, data = daughters)
plot1 = bal.plot(w.out_daughter, "dems", which = "unadjusted")
plot2 = bal.plot(w.out_daughter, "repubs", which = "unadjusted")
plot3 = bal.plot(w.out_daughter, "christian", which = "unadjusted")
plot4 = bal.plot(w.out_daughter, "age", which = "unadjusted")
plot5 = bal.plot(w.out_daughter, "srvlng", which = "unadjusted")
plot6 = bal.plot(w.out_daughter, "demvote", which = "unadjusted")
grid.arrange(plot1, plot2, plot3, plot4, plot5, plot6,nrow = 2)
MatchBalance(hasgirls ~ dems + repubs + christian + age + srvlng + demvote, data = daughters)##
## ***** (V1) dems *****
## before matching:
## mean treatment........ 0.45833
## mean control.......... 0.50847
## std mean diff......... -10.047
##
## mean raw eQQ diff..... 0.050847
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.025071
## med eCDF diff........ 0.025071
## max eCDF diff........ 0.050141
##
## var ratio (Tr/Co)..... 0.98809
## T-test p-value........ 0.35571
##
##
## ***** (V2) repubs *****
## before matching:
## mean treatment........ 0.53846
## mean control.......... 0.49153
## std mean diff......... 9.4
##
## mean raw eQQ diff..... 0.042373
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.023468
## med eCDF diff........ 0.023468
## max eCDF diff........ 0.046936
##
## var ratio (Tr/Co)..... 0.98911
## T-test p-value........ 0.3873
##
##
## ***** (V3) christian *****
## before matching:
## mean treatment........ 0.9391
## mean control.......... 0.94915
## std mean diff......... -4.1958
##
## mean raw eQQ diff..... 0.016949
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.005025
## med eCDF diff........ 0.005025
## max eCDF diff........ 0.01005
##
## var ratio (Tr/Co)..... 1.1787
## T-test p-value........ 0.68107
##
##
## ***** (V4) age *****
## before matching:
## mean treatment........ 52.628
## mean control.......... 49.178
## std mean diff......... 38.385
##
## mean raw eQQ diff..... 3.661
## med raw eQQ diff..... 4
## max raw eQQ diff..... 7
##
## mean eCDF diff........ 0.075348
## med eCDF diff........ 0.075538
## max eCDF diff........ 0.17807
##
## var ratio (Tr/Co)..... 0.71552
## T-test p-value........ 0.0020402
## KS Bootstrap p-value.. 0.004
## KS Naive p-value...... 0.0087659
## KS Statistic.......... 0.17807
##
##
## ***** (V5) srvlng *****
## before matching:
## mean treatment........ 8.5865
## mean control.......... 8.7458
## std mean diff......... -2.1085
##
## mean raw eQQ diff..... 0.66949
## med raw eQQ diff..... 0
## max raw eQQ diff..... 5
##
## mean eCDF diff........ 0.017181
## med eCDF diff........ 0.01445
## max eCDF diff........ 0.051608
##
## var ratio (Tr/Co)..... 0.77347
## T-test p-value........ 0.85956
## KS Bootstrap p-value.. 0.792
## KS Naive p-value...... 0.97653
## KS Statistic.......... 0.051608
##
##
## ***** (V6) demvote *****
## before matching:
## mean treatment........ 0.49929
## mean control.......... 0.50602
## std mean diff......... -5.2747
##
## mean raw eQQ diff..... 0.011441
## med raw eQQ diff..... 0.01
## max raw eQQ diff..... 0.08
##
## mean eCDF diff........ 0.015928
## med eCDF diff........ 0.010811
## max eCDF diff........ 0.048512
##
## var ratio (Tr/Co)..... 1.1269
## T-test p-value........ 0.61103
## KS Bootstrap p-value.. 0.948
## KS Naive p-value...... 0.98776
## KS Statistic.......... 0.048512
##
##
## Before Matching Minimum p.value: 0.0020402
## Variable Name(s): age Number(s): 4The treatment effect is: -0.4522, and the confidence interval is (-4.19406, 3.289525). Since this confidence interval includes 0, we can conclude that the treatment effect is not statistically significant.
However, we can look at the bal.plot and the MatchBalance to see the difference in distribution of the covariates. For instance, we can see a clear difference in distribution for the variable age before matching. This is confirmed when we use the MatchBalance to find the minimum p-value, which gives a statically significant value of 0.0020402 for the variable ‘age’.
STEP 2
Then, do genetic matching. Use the same variables as in the regression above. Make sure to set estimand = "ATT". What’s your treatment effect?
Note: For replicability purposes, we need to choose a see for the GenMatch() function. However, setting seed for GenMatch() is different. The usual practice of typing set.seed(some number) before the GenMatch line doesn’t ensure stochastic stability. To set seeds for GenMatch, you have to run GenMatch() including instructions for genoud within GenMatch(), e.g.: GenMatch(Tr, X, unif.seed = 123, int.seed = 92485...). You can find info on these .seed elements in the documentation for genoud(). The special .seed elements should only be present in the GenMatch() line, not elsewhere (because nothing besides GenMatch() runs genoud.
Note: When you use the GenMatch() function, wrap everything inside the following function invisible(capture.output()). This will reduce the unnecessary output produced from the GenMatch() function. For instance, you can say: invisible(capture.output(genout_daughters <- GenMatch(...)))
# Your code here
#specifying our outcome variable Y, treatment variable Tr and covariates X
attach(daughters)
Y = nowtot
Tr = hasgirls
X = cbind(dems, repubs, christian, age, srvlng, demvote)
detach(daughters)
head(X)## dems repubs christian age srvlng demvote
## [1,] 1 0 1 58 7 0.57
## [2,] 1 0 0 54 15 0.60
## [3,] 0 1 1 31 1 0.43
## [4,] 1 0 1 51 1 0.52
## [5,] 1 0 1 39 7 0.59
## [6,] 0 1 1 68 27 0.28#using genetic matching
invisible(capture.output(genout_daughters <- GenMatch(Tr = Tr, X = X, estimand = "ATT", unif.seed = 123, int.seed = 123)))
mout_daughters <- Match(Y = Y, Tr = Tr, X = X, estimand = "ATT", Weight.matrix = genout_daughters)
summary(mout_daughters)##
## Estimate... 0.59295
## AI SE...... 2.1834
## T-stat..... 0.27157
## p.val...... 0.78596
##
## Original number of observations.............. 430
## Original number of treated obs............... 312
## Matched number of observations............... 312
## Matched number of observations (unweighted). 312After matching, the estimated treatment effect is 0.592 i.e. Having a girl increases the legislator’s support for National Organization of Women by 0.59 unit. However, given that the standard error is 2.18, and the p-value is 0.786, this is not a statistically significant result.
STEP 3
Summarize (in 5-15 sentences) the genetic matching procedure and results, including what you matched on, what you balanced on, and what your balance results were. Provide output for MatchBalance() in the body of your submission.
GenMatch() matches on the control versus treatment group for hasgirls, in this case, control group is people who do not have a girl, and treatment group is people who do. It attempts to balance them based on the 6 covariates that were also used in the linear regression model, namely dems, repubs, christian, age, srvlng and demvote. It uses the genetic algorithm to find the appropriate scales such that the minimum p-value of the difference between the treatment and control group is the highest that it can be. This is because a low p-value suggests that the difference is statistically significant, and make the covariate a potential confounder for both outcome and treatment, thus we want our p-value to be as high as possible to ensure there is no significant difference between the treatment and control group. GenMatch helps us achieve that, because before matching, the lowest p-value is 0.0020402 for the covariate age, which is smaller than the standard 0.05 level of significance and indicates there is statistically significant different between the two groups in terms of age. After matching, the lowest p-value is 0.31731 for the covariate ‘dem’, which means there’s no longer a statistically significant difference in any of the covariates (since the lowest p-value is still larger than 0.05). Since as a general rule of thumb, matching is achieved when all p-values are larger than 0.15, this means the genetic matching procedure produces a good match, which we can then use to estimate the treatment effect for the treated.
#checking the matching balance
mb_daughters <- MatchBalance(hasgirls ~ dems + repubs + christian + age + srvlng + demvote, data = daughters, match.out = mout_daughters, nboots =1000)##
## ***** (V1) dems *****
## Before Matching After Matching
## mean treatment........ 0.45833 0.45833
## mean control.......... 0.50847 0.46154
## std mean diff......... -10.047 -0.64223
##
## mean raw eQQ diff..... 0.050847 0.0032051
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.025071 0.0016026
## med eCDF diff........ 0.025071 0.0016026
## max eCDF diff........ 0.050141 0.0032051
##
## var ratio (Tr/Co)..... 0.98809 0.99897
## T-test p-value........ 0.35571 0.31731
##
##
## ***** (V2) repubs *****
## Before Matching After Matching
## mean treatment........ 0.53846 0.53846
## mean control.......... 0.49153 0.53846
## std mean diff......... 9.4 0
##
## mean raw eQQ diff..... 0.042373 0
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 0
##
## mean eCDF diff........ 0.023468 0
## med eCDF diff........ 0.023468 0
## max eCDF diff........ 0.046936 0
##
## var ratio (Tr/Co)..... 0.98911 1
## T-test p-value........ 0.3873 1
##
##
## ***** (V3) christian *****
## Before Matching After Matching
## mean treatment........ 0.9391 0.9391
## mean control.......... 0.94915 0.9391
## std mean diff......... -4.1958 0
##
## mean raw eQQ diff..... 0.016949 0
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 0
##
## mean eCDF diff........ 0.005025 0
## med eCDF diff........ 0.005025 0
## max eCDF diff........ 0.01005 0
##
## var ratio (Tr/Co)..... 1.1787 1
## T-test p-value........ 0.68107 1
##
##
## ***** (V4) age *****
## Before Matching After Matching
## mean treatment........ 52.628 52.628
## mean control.......... 49.178 52.519
## std mean diff......... 38.385 1.2124
##
## mean raw eQQ diff..... 3.661 0.5641
## med raw eQQ diff..... 4 1
## max raw eQQ diff..... 7 4
##
## mean eCDF diff........ 0.075348 0.012251
## med eCDF diff........ 0.075538 0.0064103
## max eCDF diff........ 0.17807 0.038462
##
## var ratio (Tr/Co)..... 0.71552 0.99112
## T-test p-value........ 0.0020402 0.47987
## KS Bootstrap p-value.. 0.007 0.9
## KS Naive p-value...... 0.0087659 0.97513
## KS Statistic.......... 0.17807 0.038462
##
##
## ***** (V5) srvlng *****
## Before Matching After Matching
## mean treatment........ 8.5865 8.5865
## mean control.......... 8.7458 8.6955
## std mean diff......... -2.1085 -1.443
##
## mean raw eQQ diff..... 0.66949 0.5
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 5 9
##
## mean eCDF diff........ 0.017181 0.013067
## med eCDF diff........ 0.01445 0.011218
## max eCDF diff........ 0.051608 0.051282
##
## var ratio (Tr/Co)..... 0.77347 0.94914
## T-test p-value........ 0.85956 0.55124
## KS Bootstrap p-value.. 0.8 0.512
## KS Naive p-value...... 0.97653 0.80655
## KS Statistic.......... 0.051608 0.051282
##
##
## ***** (V6) demvote *****
## Before Matching After Matching
## mean treatment........ 0.49929 0.49929
## mean control.......... 0.50602 0.49936
## std mean diff......... -5.2747 -0.0503
##
## mean raw eQQ diff..... 0.011441 0.0091026
## med raw eQQ diff..... 0.01 0.01
## max raw eQQ diff..... 0.08 0.08
##
## mean eCDF diff........ 0.015928 0.013787
## med eCDF diff........ 0.010811 0.0096154
## max eCDF diff........ 0.048512 0.044872
##
## var ratio (Tr/Co)..... 1.1269 1.1284
## T-test p-value........ 0.61103 0.97778
## KS Bootstrap p-value.. 0.92 0.803
## KS Naive p-value...... 0.98776 0.91194
## KS Statistic.......... 0.048512 0.044872
##
##
## Before Matching Minimum p.value: 0.0020402
## Variable Name(s): age Number(s): 4
##
## After Matching Minimum p.value: 0.31731
## Variable Name(s): dems Number(s): 1STEP 4
Is your treatment effect different from the one reported before matching? By how much? If your numbers are different, provide some explanation as to why the two numbers are different. If they’re not, provide an explanation why they’re not different.
Yes, they differ. GenMatch gives an estimate of 0.78526 , whereas before matching using linear regression gets an estimate of -0.4522678, so there is both a difference in values (over 1.2 point difference) and in terms of effect (before matching, having a girl will have a negative impact on support for the National Organization of Women, after matching, it has a positive impact). This is because before matching, there are a lot of confounding variables. Looking at the result of the linear model, we find that the covariates christian, repubs and demvotes all have significant p-values, which mean they also influence the outcome Y. Without matching, they also influence the treatment X (because the mean difference is large, even though the p-value is not statistically significant). By matching on these covariates, we remove the confounders and give a more accurate estimation of the actual treatment effect on the treated i.e. having a girl on support for National Organization of Women.
STEP 5
Change the parameters in your genetic matching algorithm to improve the balance of the covariates. Consider rerunning with M = 2 or 3 or more. Consider increasing other parameters in the GenMatch() function such as the number of generations and population size, caliper, etc. Try 10 different ways but don’t report the results of the genetic matching weights or the balance in this document. Only report the treatment effect of your top 3 matches. For instance, run the Match() function three times for your top 3 genout objects. Make sure the summary reports the treatment effect estimate, the standard error, and the confidence interval. Do you see a large variation in your treatment effect between your top 3 models?
Note: For replicability purposes, we need to choose a see for the GenMatch() function. However, setting seed for GenMatch() is different. The usual practice of typing set.seed(123) before the GenMatch line doesn’t ensure stochastic stability. To set seeds for GenMatch, you have to run GenMatch() including instructions for genoud within GenMatch(), e.g.: GenMatch(Tr, X, unif.seed = 123, int.seed = 92485...). You can find info on these .seed elements in the documentation for genoud(). The special .seed elements should only be present in the GenMatch() line, not elsewhere (because nothing besides GenMatch() runs genoud.
Note: When you use the GenMatch() function, wrap everything inside the following function invisible(capture.output()). This will reduce the unnecessary output produced with the GenMatch() function. For instance, you can say: invisible(capture.output(genout_daughters <- GenMatch(...)))
Note: In the matching assignment, you may find that the Genetic Matching step takes a while, e.g., hours. If you have to reduce pop.size to e.g., 10 or 16 to ensure it stops after only an hour or two, that’s fine. Running your computer for an hour or two is a good thing. Running it for a full day or more is unnecessary overkill (and if this is your situation, change hyperparameters like pop.size to reduce run-time). For example, we suggest you modify the pop.size (e.g., you can set it to 5!), max.generations (set it to 2!), and wait.generations (set it to 1!) and that should expedite things.
Note: Can you set a caliper for one confounding variable, and not others (e.g., only set a caliper for “age”)? No and yes. No, strictly speaking you can’t. But yes, practically speaking you can, if you set other calipers (for the other confounders) that are so wide as to not induce any constraints. E.g., in GenMatch, and Match, you could set caliper = c(1e16, 1e16, 0.5, 1e16) and this would induce a certain meaningful caliper for the third confounder in X, without constraining the other confounders (because 1e16 implies a caliper that is so enormously wide that it does not, in practical terms, serve as a caliper at all).
# Your code here
result_summary = function(match_output){
summary(match_output)
ci_upper <- match_output$est + 1.96*match_output$se
ci_lower <- match_output$est - 1.96*match_output$se
sprintf("Estimate treatment effect: %.2f. Confidence interval: (%.2f, %.2f)", match_output$est, ci_upper, ci_lower)
}
print("5th trial results: ")## [1] "5th trial results: "invisible(capture.output(genout_daughters5 <- GenMatch(Tr = Tr, X = X, caliper = 0.2, estimand = "ATT", unif.seed = 123, int.seed = 92485)))
mout_daughters5 <- Match(Y = Y, Tr = Tr, X = X, estimand = "ATT", caliper = 0.2, Weight.matrix = genout_daughters5)
result_summary(mout_daughters5)##
## Estimate... -2.6667
## AI SE...... 0.39094
## T-stat..... -6.8212
## p.val...... 9.0268e-12
##
## Original number of observations.............. 430
## Original number of treated obs............... 312
## Matched number of observations............... 30
## Matched number of observations (unweighted). 30
##
## Caliper (SDs)........................................ 0.2 0.2 0.2 0.2 0.2 0.2
## Number of obs dropped by 'exact' or 'caliper' 282## [1] "Estimate treatment effect: -2.67. Confidence interval: (-1.90, -3.43)"print("6th trial results: ")## [1] "6th trial results: "invisible(capture.output(genout_daughters6 <- GenMatch(Tr = Tr, X = X, caliper = 0.2, estimand = "ATT", unif.seed = 123, int.seed = 92485, wait.generations = 5)))
mout_daughters6 <- Match(Y = Y, Tr = Tr, X = X, estimand = "ATT", caliper = 0.2, Weight.matrix = genout_daughters6)
result_summary(mout_daughters6)##
## Estimate... -2.6667
## AI SE...... 0.39094
## T-stat..... -6.8212
## p.val...... 9.0268e-12
##
## Original number of observations.............. 430
## Original number of treated obs............... 312
## Matched number of observations............... 30
## Matched number of observations (unweighted). 30
##
## Caliper (SDs)........................................ 0.2 0.2 0.2 0.2 0.2 0.2
## Number of obs dropped by 'exact' or 'caliper' 282## [1] "Estimate treatment effect: -2.67. Confidence interval: (-1.90, -3.43)"print("7th trial results: ")## [1] "7th trial results: "invisible(capture.output(genout_daughters7 <- GenMatch(Tr = Tr, X = X, caliper = 0.2, estimand = "ATT", unif.seed = 123, int.seed = 92485, wait.generations = 5, pop.size = 200)))
mout_daughters7 <- Match(Y = Y, Tr = Tr, X = X, estimand = "ATT", caliper = 0.2, Weight.matrix = genout_daughters7)
result_summary(mout_daughters7)##
## Estimate... -2.6667
## AI SE...... 0.39094
## T-stat..... -6.8212
## p.val...... 9.0268e-12
##
## Original number of observations.............. 430
## Original number of treated obs............... 312
## Matched number of observations............... 30
## Matched number of observations (unweighted). 30
##
## Caliper (SDs)........................................ 0.2 0.2 0.2 0.2 0.2 0.2
## Number of obs dropped by 'exact' or 'caliper' 282## [1] "Estimate treatment effect: -2.67. Confidence interval: (-1.90, -3.43)"There is no variation in the 3 best models: they all give the estimated treatment effect to be -2.67, and the same confidence interval (-1.90, -3.43). This is probably because they all end up with the same dataset (which only has 30 matched observations). However, there minimum p-value is really really high (0.77 for all 3) which indicates this dataset is very balanced.
STEP 6
Repeat everything you’ve done for Steps 1-2, including the regression, genetic algorithm, code and estimating the treatment effect EXCEPT this time change the definition of treatment to cover 2 girls or more, and change the definition of control to cover 2 boys or more. Exclude all observations that don’t meet these requirements. Be sure to explain (in a sentence or two) what you’re doing with your new treatment and control definitions. Do your new definitions change anything?
Note: Definition of the new treatment variable is as follows: Individuals in the treatment group should be having 2 or more girls and no boys, and individuals in the control group should be having 2 or more boys and no girls. What I had in mind was that such a situation increased the “dosage” of treatment vs. the “dosage” of control (and Rosenbaum alluded to this kind of observational design logic in one of the recently assigned articles). Note that you can’t have the same units in the treatment group AND the control group – we should all know by now that such a situation would be wrong.
# Your code here
#preprocessing
#finding the data in treatment group
X_treatment <- daughters[daughters$ngirls >= 2 & daughters$nboys == 0, ]
#finding the data in control group
X_control <- daughters[daughters$nboys >= 2 & daughters$ngirls ==0, ]
#creating a new dataset with only the treatment and control group
daughters_mod <- rbind(X_treatment, X_control)
#creating a new variable defining whether each person is in treatment or control
daughters_mod$treat <- ifelse(daughters_mod$ngirls >=2, 1, 0)
head(daughters_mod)#linear regression
lm_mod <- lm(nowtot ~ treat + dems + repubs + christian + age + srvlng + demvote, data = daughters_mod)
summary(lm_mod)##
## Call:
## lm(formula = nowtot ~ treat + dems + repubs + christian + age +
## srvlng + demvote, data = daughters_mod)
##
## Residuals:
## Min 1Q Median 3Q Max
## -51.636 -6.608 0.928 8.992 46.235
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -18.4332 16.3697 -1.126 0.26335
## treat 12.2925 3.5008 3.511 0.00072 ***
## dems 49.9284 4.4002 11.347 < 2e-16 ***
## repubs NA NA NA NA
## christian -3.4303 7.7211 -0.444 0.65798
## age -0.2558 0.2180 -1.173 0.24395
## srvlng 0.3908 0.2654 1.473 0.14461
## demvote 86.6631 17.9813 4.820 6.32e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.82 on 84 degrees of freedom
## Multiple R-squared: 0.839, Adjusted R-squared: 0.8275
## F-statistic: 72.97 on 6 and 84 DF, p-value: < 2.2e-16#treatment effect estimate
lm_mod$coefficients['treat']## treat
## 12.29254#finding confidence interval
confint(lm_mod,'treat',level=0.95)## 2.5 % 97.5 %
## treat 5.330876 19.25421attach(daughters_mod)
Y = nowtot
Tr = treat
X = cbind(dems, repubs, christian, age, srvlng, demvote)
detach(daughters_mod)
head(X)## dems repubs christian age srvlng demvote
## [1,] 1 0 1 51 1 0.52
## [2,] 1 0 1 39 7 0.59
## [3,] 0 1 1 70 11 0.31
## [4,] 1 0 1 38 5 0.71
## [5,] 1 0 1 37 3 0.51
## [6,] 1 0 0 55 15 0.65#using genetic matching
invisible(capture.output(genout_daughters_mod <- GenMatch(Tr = Tr, X = X, estimand = "ATT", unif.seed = 123, int.seed = 92485)))
mout_daughters_mod <- Match(Y = Y, Tr = Tr, X = X, estimand = "ATT", Weight.matrix = genout_daughters_mod)
summary(mout_daughters_mod)##
## Estimate... 13.298
## AI SE...... 4.2203
## T-stat..... 3.1509
## p.val...... 0.0016274
##
## Original number of observations.............. 91
## Original number of treated obs............... 47
## Matched number of observations............... 47
## Matched number of observations (unweighted). 47mb_daughters_mod <- MatchBalance(treat ~ dems + repubs + christian + age + srvlng + demvote, data = daughters_mod, match.out = mout_daughters_mod)##
## ***** (V1) dems *****
## Before Matching After Matching
## mean treatment........ 0.61702 0.61702
## mean control.......... 0.40909 0.61702
## std mean diff......... 42.317 0
##
## mean raw eQQ diff..... 0.20455 0
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 0
##
## mean eCDF diff........ 0.10397 0
## med eCDF diff........ 0.10397 0
## max eCDF diff........ 0.20793 0
##
## var ratio (Tr/Co)..... 0.97609 1
## T-test p-value........ 0.04806 1
##
##
## ***** (V2) repubs *****
## Before Matching After Matching
## mean treatment........ 0.38298 0.38298
## mean control.......... 0.59091 0.38298
## std mean diff......... -42.317 0
##
## mean raw eQQ diff..... 0.22727 0
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 0
##
## mean eCDF diff........ 0.10397 0
## med eCDF diff........ 0.10397 0
## max eCDF diff........ 0.20793 0
##
## var ratio (Tr/Co)..... 0.97609 1
## T-test p-value........ 0.04806 1
##
##
## ***** (V3) christian *****
## Before Matching After Matching
## mean treatment........ 0.91489 0.91489
## mean control.......... 0.97727 0.91489
## std mean diff......... -22.116 0
##
## mean raw eQQ diff..... 0.068182 0
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 0
##
## mean eCDF diff........ 0.03119 0
## med eCDF diff........ 0.03119 0
## max eCDF diff........ 0.062379 0
##
## var ratio (Tr/Co)..... 3.5005 1
## T-test p-value........ 0.1887 1
##
##
## ***** (V4) age *****
## Before Matching After Matching
## mean treatment........ 51.213 51.213
## mean control.......... 51.977 51.66
## std mean diff......... -8.2171 -4.8024
##
## mean raw eQQ diff..... 1.5682 1.1277
## med raw eQQ diff..... 1 1
## max raw eQQ diff..... 6 4
##
## mean eCDF diff........ 0.0288 0.03125
## med eCDF diff........ 0.025629 0.021277
## max eCDF diff........ 0.07882 0.10638
##
## var ratio (Tr/Co)..... 0.79756 1.0521
## T-test p-value........ 0.71354 0.75411
## KS Bootstrap p-value.. 0.966 0.886
## KS Naive p-value...... 0.99893 0.953
## KS Statistic.......... 0.07882 0.10638
##
##
## ***** (V5) srvlng *****
## Before Matching After Matching
## mean treatment........ 7.9574 7.9574
## mean control.......... 10.568 7.8936
## std mean diff......... -38.919 0.95154
##
## mean raw eQQ diff..... 2.8636 0.70213
## med raw eQQ diff..... 2 0
## max raw eQQ diff..... 14 2
##
## mean eCDF diff........ 0.079008 0.024113
## med eCDF diff........ 0.091876 0.021277
## max eCDF diff........ 0.13926 0.06383
##
## var ratio (Tr/Co)..... 0.48803 0.89279
## T-test p-value........ 0.13925 0.83488
## KS Bootstrap p-value.. 0.532 0.99
## KS Naive p-value...... 0.77019 0.99998
## KS Statistic.......... 0.13926 0.06383
##
##
## ***** (V6) demvote *****
## Before Matching After Matching
## mean treatment........ 0.51745 0.51745
## mean control.......... 0.48727 0.51191
## std mean diff......... 24.982 4.58
##
## mean raw eQQ diff..... 0.045909 0.025106
## med raw eQQ diff..... 0.045 0.02
## max raw eQQ diff..... 0.11 0.07
##
## mean eCDF diff........ 0.10233 0.050236
## med eCDF diff........ 0.091393 0.042553
## max eCDF diff........ 0.25193 0.12766
##
## var ratio (Tr/Co)..... 1.0419 1.2252
## T-test p-value........ 0.23199 0.56203
## KS Bootstrap p-value.. 0.062 0.7
## KS Naive p-value...... 0.11171 0.83838
## KS Statistic.......... 0.25193 0.12766
##
##
## Before Matching Minimum p.value: 0.04806
## Variable Name(s): dems repubs Number(s): 1 2
##
## After Matching Minimum p.value: 0.56203
## Variable Name(s): demvote Number(s): 6After matching, the estimated treatment effect is 13.298 and the p-value is 0.001, which means there is a statistically significant different between the control group (having 2+ boys and no girls) and treatment group (having 2+ girls and no boys) i.e. being “treated” will increase the average support for the National Organization of Women by 13.298 units. This is not too different from the linear regression model estimate (12.29) but we have reasons to believe in the Matching values more, because after matching the minimum p-value is 0.562, which is larger than the standard requirement that all p-values be larger than 0.15.
STEP 7
It is NOT wise to match or balance on “totchi”. What is the reason? Hint: You will have to look at what variables mean in the data set to be able to answer this question.
Totchi is the total children. This is a variable that is influenced by total number of girls and total number of boys in the dataset. However, we are using the number of girls and number of boys as our treatment units. So matching “totchi” would mean matching the treatment variable, which does not make sense because treatment variable is not supposed to be balanced.
QUESTION 3: COPD
Most causal studies on the health effects of smoking are observational studies (well, for very obvious reasons). In this exercise, we are specifically after answer the following question: Does smoking increase the risk of chronic obstructive pulmonary disease (COPD)? To learn more about the disease, read here: https://www.cdc.gov/copd/index.html
We’ll use a sub-sample of the 2015 BRFSS survey (pronounced bur-fiss), which stands for Behavioral Risk Factor Surveillance System. The data is collected through a phone survey across American citizens regarding their health-related risk behaviors and chronic health conditions. Although, the entire survey has over 400,000 records and over 300 variables, we only sample 5,000 observations and 7 variables.
Let’s load the data first and take a look at the first few rows:
brfss = read.csv("http://bit.ly/BRFSS_data") %>%
clean_names()
head(brfss)A summary of the variables is as follows:
- copd: Ever told you have chronic obstructive pulmonary disease (COPD)?
- smoke: Adults who are current smokers (0 = no, 1 = yes)
- race: Race group
- age: age group
- sex: gender
- wtlbs: weight in pounds (lbs)
- avedrnk2: During the past 30 days, when you drank, how many drinks did you drink on average?
STEP 1
Check the balance of covariates before matching using any method of your choice. You can look at balance tables, balance plots, or love plots from any package of your choice. Do you see a balance across the covariates?
Note: This is optional but you can use the gridExtra package and its grid.arrange() function to put all the 4 graphs in one 2 x 2 graph. Read more about the package and how to use it here: https://cran.r-project.org/web/packages/egg/vignettes/Ecosystem.html. Set nrow = 2.
# Your code here
w.out2 <- weightit(smoke~race + age + sex + wtlbs + avedrnk2, data = brfss, estimand = "ATT", m.threshold = 4)
plot1a = bal.plot(w.out2, "race", which = "unadjusted")
plot2a = bal.plot(w.out2, "age", which = "unadjusted")
plot3a = bal.plot(w.out2, "sex", which = "unadjusted")
plot4a = bal.plot(w.out2, "wtlbs", which = "unadjusted")
plot5a = bal.plot(w.out2, "avedrnk2", which = "unadjusted")
grid.arrange(plot1a, plot2a, plot3a, plot4a, plot5a, nrow = 3)
MatchBalance(smoke~race + age + sex + wtlbs + avedrnk2, data = brfss)##
## ***** (V1) raceHispanic *****
## before matching:
## mean treatment........ 0.070866
## mean control.......... 0.067249
## std mean diff......... 1.4088
##
## mean raw eQQ diff..... 0.0026247
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.0018087
## med eCDF diff........ 0.0018087
## max eCDF diff........ 0.0036174
##
## var ratio (Tr/Co)..... 1.0508
## T-test p-value........ 0.71939
##
##
## ***** (V2) raceOther *****
## before matching:
## mean treatment........ 0.068241
## mean control.......... 0.047664
## std mean diff......... 8.1551
##
## mean raw eQQ diff..... 0.019685
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.010289
## med eCDF diff........ 0.010289
## max eCDF diff........ 0.020577
##
## var ratio (Tr/Co)..... 1.4023
## T-test p-value........ 0.034312
##
##
## ***** (V3) raceWhite *****
## before matching:
## mean treatment........ 0.74147
## mean control.......... 0.83388
## std mean diff......... -21.094
##
## mean raw eQQ diff..... 0.091864
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.046207
## med eCDF diff........ 0.046207
## max eCDF diff........ 0.092414
##
## var ratio (Tr/Co)..... 1.3853
## T-test p-value........ 5.4811e-08
##
##
## ***** (V4) age25-34 *****
## before matching:
## mean treatment........ 0.18241
## mean control.......... 0.10854
## std mean diff......... 19.116
##
## mean raw eQQ diff..... 0.073491
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.036936
## med eCDF diff........ 0.036936
## max eCDF diff........ 0.073873
##
## var ratio (Tr/Co)..... 1.543
## T-test p-value........ 7.0411e-07
##
##
## ***** (V5) age35-44 *****
## before matching:
## mean treatment........ 0.1811
## mean control.......... 0.12482
## std mean diff......... 14.605
##
## mean raw eQQ diff..... 0.05643
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.02814
## med eCDF diff........ 0.02814
## max eCDF diff........ 0.056279
##
## var ratio (Tr/Co)..... 1.359
## T-test p-value........ 0.00016078
##
##
## ***** (V6) age45-54 *****
## before matching:
## mean treatment........ 0.18373
## mean control.......... 0.18263
## std mean diff......... 0.28224
##
## mean raw eQQ diff..... 0.0013123
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.00054686
## med eCDF diff........ 0.00054686
## max eCDF diff........ 0.0010937
##
## var ratio (Tr/Co)..... 1.0057
## T-test p-value........ 0.94281
##
##
## ***** (V7) age55-64 *****
## before matching:
## mean treatment........ 0.21916
## mean control.......... 0.22534
## std mean diff......... -1.4934
##
## mean raw eQQ diff..... 0.0065617
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.003091
## med eCDF diff........ 0.003091
## max eCDF diff........ 0.006182
##
## var ratio (Tr/Co)..... 0.98138
## T-test p-value........ 0.70477
##
##
## ***** (V8) age65+ *****
## before matching:
## mean treatment........ 0.16404
## mean control.......... 0.31029
## std mean diff......... -39.467
##
## mean raw eQQ diff..... 0.14698
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.073123
## med eCDF diff........ 0.073123
## max eCDF diff........ 0.14625
##
## var ratio (Tr/Co)..... 0.64147
## T-test p-value........ < 2.22e-16
##
##
## ***** (V9) sexMale *****
## before matching:
## mean treatment........ 0.54593
## mean control.......... 0.48679
## std mean diff......... 11.872
##
## mean raw eQQ diff..... 0.059055
## med raw eQQ diff..... 0
## max raw eQQ diff..... 1
##
## mean eCDF diff........ 0.029573
## med eCDF diff........ 0.029573
## max eCDF diff........ 0.059146
##
## var ratio (Tr/Co)..... 0.99332
## T-test p-value........ 0.002627
##
##
## ***** (V10) wtlbs *****
## before matching:
## mean treatment........ 177.76
## mean control.......... 177.99
## std mean diff......... -0.49852
##
## mean raw eQQ diff..... 2.1432
## med raw eQQ diff..... 0
## max raw eQQ diff..... 50
##
## mean eCDF diff........ 0.0077522
## med eCDF diff........ 0.0056792
## max eCDF diff........ 0.025879
##
## var ratio (Tr/Co)..... 1.1223
## T-test p-value........ 0.89836
## KS Bootstrap p-value.. 0.692
## KS Naive p-value...... 0.78001
## KS Statistic.......... 0.025879
##
##
## ***** (V11) avedrnk2 *****
## before matching:
## mean treatment........ 2.9843
## mean control.......... 1.9672
## std mean diff......... 40.633
##
## mean raw eQQ diff..... 1
## med raw eQQ diff..... 1
## max raw eQQ diff..... 4
##
## mean eCDF diff........ 0.052876
## med eCDF diff........ 0.01035
## max eCDF diff........ 0.26493
##
## var ratio (Tr/Co)..... 1.9119
## T-test p-value........ < 2.22e-16
## KS Bootstrap p-value.. < 2.22e-16
## KS Naive p-value...... 0
## KS Statistic.......... 0.26493
##
##
## Before Matching Minimum p.value: < 2.22e-16
## Variable Name(s): age65+ avedrnk2 Number(s): 8 11(My distribution for sex variable is weird) From the plots, we can see that the avedrnk2 variable and the age variable are not well-distributed. The age variable has significantly more 65+ people in the contorl than treatment group, and the avedrnk for control is a lot lower than treatment i.e. people in treatment are likely to have more drinks. This is also confirmed by the minimum p-value (< 2.22e^-16) for these 2 values.
STEP 2
Now, let’s do Mahalanobis distance matching. Note that you can use the same old Match() function. Use all covariates in the data set to match on, however, make sure you convert all categorical variables into factor variables (Google to see how). We are going to specify estimand = "ATT" in the Match() function. What’s the treatment effect after matching?
# Your code here
#redefining the age variable as the average value in that age group
brfss[brfss$age == "18-24",]$age = 21
brfss[brfss$age == "25-34",]$age = 30
brfss[brfss$age == "35-44",]$age = 40
brfss[brfss$age == "45-54",]$age = 50
brfss[brfss$age == "55-64",]$age = 60
brfss[brfss$age == "65+",]$age = 70
head(brfss)#defining outcome, treatment and covariates
attach(brfss)
Y <- ifelse(copd == "No", 0, 1)
Tr <- smoke
X <- cbind(as.factor(race), as.integer(age), as.factor(sex), wtlbs, avedrnk2)
detach(brfss)
#use Matching
mout_brfss <- Match(Y = Y, Tr = Tr, X = X, estimand = "ATT")
summary(mout_brfss)##
## Estimate... 0.086697
## AI SE...... 0.015663
## T-stat..... 5.535
## p.val...... 3.1118e-08
##
## Original number of observations.............. 5000
## Original number of treated obs............... 762
## Matched number of observations............... 762
## Matched number of observations (unweighted). 2248mb_brfss <- MatchBalance(smoke ~ race + as.integer(age) + sex + wtlbs + avedrnk2, data = brfss, match.out = mout_brfss)##
## ***** (V1) raceHispanic *****
## Before Matching After Matching
## mean treatment........ 0.070866 0.070866
## mean control.......... 0.067249 0.073491
## std mean diff......... 1.4088 -1.0222
##
## mean raw eQQ diff..... 0.0026247 0.00088968
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.0018087 0.00044484
## med eCDF diff........ 0.0018087 0.00044484
## max eCDF diff........ 0.0036174 0.00088968
##
## var ratio (Tr/Co)..... 1.0508 0.96702
## T-test p-value........ 0.71939 0.15716
##
##
## ***** (V2) raceOther *****
## Before Matching After Matching
## mean treatment........ 0.068241 0.068241
## mean control.......... 0.047664 0.062992
## std mean diff......... 8.1551 2.0804
##
## mean raw eQQ diff..... 0.019685 0.0017794
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.010289 0.00088968
## med eCDF diff........ 0.010289 0.00088968
## max eCDF diff........ 0.020577 0.0017794
##
## var ratio (Tr/Co)..... 1.4023 1.0773
## T-test p-value........ 0.034312 0.045288
##
##
## ***** (V3) raceWhite *****
## Before Matching After Matching
## mean treatment........ 0.74147 0.74147
## mean control.......... 0.83388 0.74541
## std mean diff......... -21.094 -0.89863
##
## mean raw eQQ diff..... 0.091864 0.0013345
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.046207 0.00066726
## med eCDF diff........ 0.046207 0.00066726
## max eCDF diff........ 0.092414 0.0013345
##
## var ratio (Tr/Co)..... 1.3853 1.0101
## T-test p-value........ 5.4811e-08 0.083062
##
##
## ***** (V4) as.integer(age) *****
## Before Matching After Matching
## mean treatment........ 47.996 47.996
## mean control.......... 53.637 48.214
## std mean diff......... -36.926 -1.4288
##
## mean raw eQQ diff..... 5.6352 0.12989
## med raw eQQ diff..... 10 0
## max raw eQQ diff..... 10 10
##
## mean eCDF diff........ 0.094374 0.0022242
## med eCDF diff........ 0.12065 0.0020018
## max eCDF diff........ 0.15243 0.0044484
##
## var ratio (Tr/Co)..... 1.0145 1.0185
## T-test p-value........ < 2.22e-16 0.022559
## KS Bootstrap p-value.. < 2.22e-16 0.996
## KS Naive p-value...... 1.8474e-13 1
## KS Statistic.......... 0.15243 0.0044484
##
##
## ***** (V5) sexMale *****
## Before Matching After Matching
## mean treatment........ 0.54593 0.54593
## mean control.......... 0.48679 0.54724
## std mean diff......... 11.872 -0.26341
##
## mean raw eQQ diff..... 0.059055 0.00044484
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.029573 0.00022242
## med eCDF diff........ 0.029573 0.00022242
## max eCDF diff........ 0.059146 0.00044484
##
## var ratio (Tr/Co)..... 0.99332 1.0005
## T-test p-value........ 0.002627 0.31731
##
##
## ***** (V6) wtlbs *****
## Before Matching After Matching
## mean treatment........ 177.76 177.76
## mean control.......... 177.99 178.06
## std mean diff......... -0.49852 -0.65757
##
## mean raw eQQ diff..... 2.1432 0.72645
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 50 35
##
## mean eCDF diff........ 0.0077522 0.0033568
## med eCDF diff........ 0.0056792 0.002669
## max eCDF diff........ 0.025879 0.01379
##
## var ratio (Tr/Co)..... 1.1223 1.1091
## T-test p-value........ 0.89836 0.41466
## KS Bootstrap p-value.. 0.664 0.886
## KS Naive p-value...... 0.78001 0.98312
## KS Statistic.......... 0.025879 0.01379
##
##
## ***** (V7) avedrnk2 *****
## Before Matching After Matching
## mean treatment........ 2.9843 2.9843
## mean control.......... 1.9672 2.8828
## std mean diff......... 40.633 4.0546
##
## mean raw eQQ diff..... 1 0.039591
## med raw eQQ diff..... 1 0
## max raw eQQ diff..... 4 1
##
## mean eCDF diff........ 0.052876 0.0026394
## med eCDF diff........ 0.01035 0.0013345
## max eCDF diff........ 0.26493 0.012456
##
## var ratio (Tr/Co)..... 1.9119 1.0428
## T-test p-value........ < 2.22e-16 8.7839e-11
## KS Bootstrap p-value.. < 2.22e-16 0.612
## KS Naive p-value...... < 2.22e-16 0.99492
## KS Statistic.......... 0.26493 0.012456
##
##
## Before Matching Minimum p.value: < 2.22e-16
## Variable Name(s): as.integer(age) avedrnk2 Number(s): 4 7
##
## After Matching Minimum p.value: 8.7839e-11
## Variable Name(s): avedrnk2 Number(s): 7The average treatment effect estimate after matching is 0.086697. This has a statistically significant p-value of 3.1118e^-8, and a small standard error of 0.015663.
STEP 3
Provide a few sentences talking about the number of treated units dropped, and a few more sentences talking about the balance obtained.
No treated units were dropped because there are 762 treated observations, and all 762 observations were matched. This is also because we use the estimand as “ATT” which is the average treatment effect for the treated, so we find the matches for all the treated units, and drop all the control units that are not needed. The balance obtained is not great, because the minimum p-value is still 8.78e^11 for the variable avedrnk2, which is statistically significant.
STEP 4
Now, let’s do another Mahalanobis distance matching. Use all covariates in the data set in the propensity score estimation. However, this time make sure you specify estimand = "ATE" in the Match() function. What’s the treatment effect after matching?
# Your code here
mout_brfss <- Match(Y = Y, Tr = Tr, X = X, estimand = "ATE")
summary(mout_brfss)##
## Estimate... 0.1323
## AI SE...... 0.017839
## T-stat..... 7.4163
## p.val...... 1.2057e-13
##
## Original number of observations.............. 5000
## Original number of treated obs............... 762
## Matched number of observations............... 5000
## Matched number of observations (unweighted). 7759mb_brfss <- MatchBalance(smoke ~ race + age + sex + wtlbs + avedrnk2, data = brfss, match.out = mout_brfss)##
## ***** (V1) raceHispanic *****
## Before Matching After Matching
## mean treatment........ 0.070866 0.0656
## mean control.......... 0.067249 0.0682
## std mean diff......... 1.4088 -1.0501
##
## mean raw eQQ diff..... 0.0026247 0.0016755
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.0018087 0.00083774
## med eCDF diff........ 0.0018087 0.00083774
## max eCDF diff........ 0.0036174 0.0016755
##
## var ratio (Tr/Co)..... 1.0508 0.96456
## T-test p-value........ 0.71939 0.006697
##
##
## ***** (V2) raceOther *****
## Before Matching After Matching
## mean treatment........ 0.068241 0.0494
## mean control.......... 0.047664 0.05
## std mean diff......... 8.1551 -0.27685
##
## mean raw eQQ diff..... 0.019685 0.00051553
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.010289 0.00025777
## med eCDF diff........ 0.010289 0.00025777
## max eCDF diff........ 0.020577 0.00051553
##
## var ratio (Tr/Co)..... 1.4023 0.98862
## T-test p-value........ 0.034312 0.57748
##
##
## ***** (V3) raceWhite *****
## Before Matching After Matching
## mean treatment........ 0.74147 0.8232
## mean control.......... 0.83388 0.8204
## std mean diff......... -21.094 0.73387
##
## mean raw eQQ diff..... 0.091864 0.0019332
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.046207 0.00096662
## med eCDF diff........ 0.046207 0.00096662
## max eCDF diff........ 0.092414 0.0019332
##
## var ratio (Tr/Co)..... 1.3853 0.98777
## T-test p-value........ 5.4811e-08 0.004253
##
##
## ***** (V4) age30 *****
## Before Matching After Matching
## mean treatment........ 0.18241 0.1248
## mean control.......... 0.10854 0.12147
## std mean diff......... 19.116 1.0085
##
## mean raw eQQ diff..... 0.073491 0.0023199
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.036936 0.0011599
## med eCDF diff........ 0.036936 0.0011599
## max eCDF diff........ 0.073873 0.0023199
##
## var ratio (Tr/Co)..... 1.543 1.0235
## T-test p-value........ 7.0411e-07 0.16387
##
##
## ***** (V5) age40 *****
## Before Matching After Matching
## mean treatment........ 0.1811 0.1463
## mean control.......... 0.12482 0.1314
## std mean diff......... 14.605 4.2157
##
## mean raw eQQ diff..... 0.05643 0.010182
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.02814 0.0050909
## med eCDF diff........ 0.02814 0.0050909
## max eCDF diff........ 0.056279 0.010182
##
## var ratio (Tr/Co)..... 1.359 1.0943
## T-test p-value........ 0.00016078 4.0312e-10
##
##
## ***** (V6) age50 *****
## Before Matching After Matching
## mean treatment........ 0.18373 0.1748
## mean control.......... 0.18263 0.18333
## std mean diff......... 0.28224 -2.2466
##
## mean raw eQQ diff..... 0.0013123 0.0057997
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.00054686 0.0028999
## med eCDF diff........ 0.00054686 0.0028999
## max eCDF diff........ 0.0010937 0.0057997
##
## var ratio (Tr/Co)..... 1.0057 0.96342
## T-test p-value........ 0.94281 4.1995e-05
##
##
## ***** (V7) age60 *****
## Before Matching After Matching
## mean treatment........ 0.21916 0.2281
## mean control.......... 0.22534 0.2262
## std mean diff......... -1.4934 0.45276
##
## mean raw eQQ diff..... 0.0065617 0.0015466
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.003091 0.0007733
## med eCDF diff........ 0.003091 0.0007733
## max eCDF diff........ 0.006182 0.0015466
##
## var ratio (Tr/Co)..... 0.98138 1.0059
## T-test p-value........ 0.70477 0.35041
##
##
## ***** (V8) age70 *****
## Before Matching After Matching
## mean treatment........ 0.16404 0.2812
## mean control.......... 0.31029 0.2874
## std mean diff......... -39.467 -1.3789
##
## mean raw eQQ diff..... 0.14698 0.0042531
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.073123 0.0021266
## med eCDF diff........ 0.073123 0.0021266
## max eCDF diff........ 0.14625 0.0042531
##
## var ratio (Tr/Co)..... 0.64147 0.98694
## T-test p-value........ < 2.22e-16 9.2901e-05
##
##
## ***** (V9) sexMale *****
## Before Matching After Matching
## mean treatment........ 0.54593 0.4954
## mean control.......... 0.48679 0.496
## std mean diff......... 11.872 -0.11999
##
## mean raw eQQ diff..... 0.059055 0.00038665
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.029573 0.00019332
## med eCDF diff........ 0.029573 0.00019332
## max eCDF diff........ 0.059146 0.00038665
##
## var ratio (Tr/Co)..... 0.99332 0.99998
## T-test p-value........ 0.002627 0.1797
##
##
## ***** (V10) wtlbs *****
## Before Matching After Matching
## mean treatment........ 177.76 177.32
## mean control.......... 177.99 178
## std mean diff......... -0.49852 -1.5865
##
## mean raw eQQ diff..... 2.1432 1.3949
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 50 142
##
## mean eCDF diff........ 0.0077522 0.0058235
## med eCDF diff........ 0.0056792 0.0044464
## max eCDF diff........ 0.025879 0.035572
##
## var ratio (Tr/Co)..... 1.1223 0.94923
## T-test p-value........ 0.89836 7.5639e-06
## KS Bootstrap p-value.. 0.656 < 2.22e-16
## KS Naive p-value...... 0.78001 0.00010895
## KS Statistic.......... 0.025879 0.035572
##
##
## ***** (V11) avedrnk2 *****
## Before Matching After Matching
## mean treatment........ 2.9843 2.1307
## mean control.......... 1.9672 2.1067
## std mean diff......... 40.633 1.2505
##
## mean raw eQQ diff..... 1 0.033896
## med raw eQQ diff..... 1 0
## max raw eQQ diff..... 4 6
##
## mean eCDF diff........ 0.052876 0.0017094
## med eCDF diff........ 0.01035 0.00038665
## max eCDF diff........ 0.26493 0.020492
##
## var ratio (Tr/Co)..... 1.9119 0.96929
## T-test p-value........ < 2.22e-16 6.5542e-05
## KS Bootstrap p-value.. < 2.22e-16 0.02
## KS Naive p-value...... < 2.22e-16 0.076905
## KS Statistic.......... 0.26493 0.020492
##
##
## Before Matching Minimum p.value: < 2.22e-16
## Variable Name(s): age70 avedrnk2 Number(s): 8 11
##
## After Matching Minimum p.value: < 2.22e-16
## Variable Name(s): wtlbs Number(s): 10The average treatment effect estimate after matching is 0.1323, which is statistically significant with a p-value of 1.2057e-13. #### STEP 5
Are your answers in Steps 2 and 4 different? Why? What does the matching process do differently in each step? Which answer do you trust more? Yes, these answers are different, because the first is matching only for the treated (so only 752 values), but the second is matching for both the treated and the control group (5000 values in total). I think the ATT is more accurate both empirically (because if we look at the MatchBalance, after matching ATE still has a minimum p-value <2.22e^-16 while ATT has a minimum p-value of 8.7839e-11, so ATT is better matched, and would give more accurate results) and in general (because we need to match all the observations, it would probably be a worse match).
BONUS QUESTION: Sensitivity Analysis
STEP 1
Use the BRFSS data set from the previous question. Now, identify the critical value of Gamma as we discussed in the class. Do it using rbounds: https://nbviewer.jupyter.org/gist/viniciusmss/a156c3f22081fb5c690cdd58658f61fa
# Your code here
library(rbounds)
mout_brfss$est## [,1]
## [1,] 0.132296?psens
psens(mout_brfss, Gamma=5, GammaInc=.1)$boundsSTEP 2
Then, write a paragraph explaining what you found. Your paragraph should include numbers obtained from your analysis and explain what those numbers mean as simply as you can.
After matching, we may still have covariates that we did not collect data on, observe and match, but nonetheless affect the treatment and outcome variable. Sensitivity analysis attempts to address how likely is this possibility. Since only at the gamma of 3.7 the upper bound is 0.0594 > 0.05, this is the critical value. Gamma represents the ratio in odds of treatment for 2 subjects with the same observed covariates but a different unobserved covariate, and sensitivity analysis looks at how large gamma can be before the conclusion of the study changes i.e. p-value is larger than 0.05. So in this case, one subject needs to be 3.7 times more likely as another to receive the treatment due to an unobserved covariate before the study conclusion becomes non-significant. Since this is a large number of gamma (it is unlikely to find any unobserved covariate that affects odds of treatment this much) the study is likely robust to unobserved treatment.
End of Assignment
Final Steps
Before finalizing your project you’ll want be sure there are comments in your code chunks and text outside of your code chunks to explain what you’re doing in each code chunk. These explanations are incredibly helpful for someone who doesn’t code or someone unfamiliar to your project.
You have two options for submission:
- You can complete this .rmd file, knit it to pdf and submit the resulting .pdf file on Forum.
- You can complete the Google Doc version of this assignment, include your code, graphs, results, and your explanations wherever necessary and download the Google Doc as a pdf file and submit the pdf file on Forum. If you choose this method, you need to make sure you will provide a link to an .R script file where you code can be found (you can host your code on Github or Google Drive). Note that links to Google Docs are not accepted as your final submission.
Knitting your R Markdown Document
Last but not least, you’ll want to Knit your .Rmd document into a pdf document. If you get an error, take a look at what the error says and edit your .Rmd document. Then, try to Knit again! Troubleshooting these error messages will teach you a lot about coding in R. If you get any error that doesn’t make sense to you, post it on Perusall.
Good Luck! The CS112 Teaching Team