Lab 4 Stat209 2/27/15
Rogosa Lab 4 R-session (base: secs 1-3)
Preliminaries, install MatchIt and get lab data
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
# I Started from a relatively clean install, get MatchIt
> install.packages("MatchIt")
Installing package into ‘C:/Users/rag/Documents/R/win-library/3.1’
(as ‘lib’ is unspecified)
--- Please select a CRAN mirror for use in this session ---
trying URL 'http://cran.cnr.Berkeley.edu/bin/windows/contrib/3.1/MatchIt_2.4-21.zip'
Content type 'application/zip' length 78183 bytes (76 Kb)
opened URL
downloaded 76 Kb
package ‘MatchIt’ successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\rag\AppData\Local\Temp\RtmpeUK0DY\downloaded_packages
#You could then force an install of optmatch (which MatchIt may do for you
#automatically when you first call an optmatch function depending on the state of the license)
> install.packages("optmatch")
Installing package into ‘C:/Users/rag/Documents/R/win-library/3.1’
(as ‘lib’ is unspecified)
also installing the dependency ‘digest’
trying URL 'http://cran.cnr.Berkeley.edu/bin/windows/contrib/3.1/digest_0.6.8.zip'
Content type 'application/zip' length 151920 bytes (148 Kb)
opened URL
downloaded 148 Kb
trying URL 'http://cran.cnr.Berkeley.edu/bin/windows/contrib/3.1/optmatch_0.9-3.zip'
Content type 'application/zip' length 596350 bytes (582 Kb)
opened URL
downloaded 582 Kb
package ‘digest’ successfully unpacked and MD5 sums checked
package ‘optmatch’ successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\rag\AppData\Local\Temp\RtmpeUK0DY\downloaded_packages
> library(MatchIt)
Loading required package: MASS
> library(optmatch)
Loading required package: digest
You're loading optmatch, by B. Hansen and M. Fredrickson.
The optmatch package makes essential use of D. P. Bertsekas
and P. Tseng's RELAX-IV algorithm and code, as well as
Bertsekas' AUCTION algorithm and code. Using the software
to 'satisfy in any part commercial delivery requirements to
government or industry' requires a special agreement with
Dr. Bertsekas. For more information, enter
relaxinfo() at the command line.
> data(nuclearplants) # we don't need this, just for illustration
Now to get going on our Lab 4 tasks
> data(lalonde) # get lalonde data from MatchIt, there are different versions in existence under this name (see M-B text)
# or simply
help(lalonde) #produces
-----------------------------------
lalonde package:MatchIt R Documentation
Data from National Supported Work Demonstration and PSID, as analyzed by Dehejia and Wahba (1999).
Description:
This is a subsample of the data from the treated group in the
National Supported Work Demonstration (NSW) and the comparison
sample from the Current Population Survey (CPS). This data was
previously analyzed extensively by Lalonde (1986) and Dehejia and
Wahba (1999). The full dataset is available at . [note: broken link still in current documentation]
Usage:
data(lalonde)
Format:
A data frame with 313 [sic, 614] observations (185 treated, 429 control).
There are 10 variables measured for each individual. "treat" is
the treatment assignment (1=treated, 0=control). "age" is age in
years. "educ" is education in number of years of schooling.
"black" is an indicator for African-American (1=African-American,
0=not). "hispan" is an indicator for being of Hispanic origin
(1=Hispanic, 0=not). "married" is an indicator for married
(1=married, 0=not married). "nodegree" is an indicator for
whether the individual has a high school degree (1=no degree,
0=degree). "re74" is income in 1974, in U.S. dollars. "re75" is
income in 1975, in U.S. dollars. "re78" is income in 1978, in
U.S. dollars.
References:
Lalonde, R. (1986). Evaluating the econometric evaluations of
training programs with experimental data. American Economic
Review 76: 604-620.
Dehejia, R.H. and Wahba, S. (1999). Causal Effects in
Nonexperimental Studies: Re-Evaluating the Evaluation of
Training Programs. Journal of the American Statistical
Association 94: 1053-1062.
-----------------------------------------------
> dim(lalonde) # we have 614 cases, 10 vars; see lab script part 2
[1] 614 10
> names(lalonde)
[1] "treat" "age" "educ" "black" "hispan" "married"
[7] "nodegree" "re74" "re75" "re78"
> attach(lalonde)
> table(treat) # so these summaries synch with data description
treat
0 1
429 185
> lalonde[1:10,] # here's the first 10 observations, all from the treatment
treat age educ black hispan married nodegree re74 re75 re78
NSW1 1 37 11 1 0 1 1 0 0 9930.0
NSW2 1 22 9 0 1 0 1 0 0 3595.9
NSW3 1 30 12 1 0 0 0 0 0 24909.5
NSW4 1 27 11 1 0 0 1 0 0 7506.1
NSW5 1 33 8 1 0 0 1 0 0 289.8
NSW6 1 22 9 1 0 0 1 0 0 4056.5
NSW7 1 23 12 1 0 0 0 0 0 0.0
NSW8 1 32 11 1 0 0 1 0 0 8472.2
NSW9 1 22 16 1 0 0 0 0 0 2164.0
NSW10 1 33 12 0 0 1 0 0 0 12418.1
> # brute-force matching-- Lab4 script item 3.1
> numObs = dim(lalonde)[1]
> numObs
[1] 614
> # that's the 429 + 185, not going to distinguish between T/C here
#Do it yourself matching function, Lab text part 3#######
> indiv1 = lalonde[1,]
> indiv1 # the target for this matching exercise
treat age educ black hispan married nodegree re74 re75 re78
NSW1 1 37 11 1 0 1 1 0 0 9930
> # aside indiv1 is in the treatment group, the first 185 cases are all in treatment
> # for example, examine lalonde[185,] lalonde[186,]
> matches = matrix( nrow = numObs, ncol = 1, data = FALSE )
> # set up a dummy matrix and then loop through the data with specified calipers
> for ( i in 1:numObs) {
+ matches[i,1] <- all(
+ c(lalonde[i,4:7] == indiv1[4:7],
+ (lalonde[i,3] < indiv1[3] + 3) && (lalonde[i,3] > indiv1[3] - 3),
+ (lalonde[i,2] < indiv1[2] + 3) && (lalonde[i,2] > indiv1[2] - 3)
+ ) )
+ }
# Karen's function picks out exact matches on the four dichotomous black hisp married nodegree
and caliper matching on age and educ
> sum(as.numeric(matches) )
[1] 4
# this picks up a match to itself (the target indiv1)
> # we only got 3 matches to others cases 183, 191 212, the latter 2 being control group members
> matches # you can print out all 614
[,1]
[1,] TRUE
[2,] FALSE
[182,] FALSE
[183,] TRUE
[190,] FALSE
[191,] TRUE
[211,] FALSE
[212,] TRUE
[613,] FALSE
[614,] FALSE
> lalonde[191,]
treat age educ black hispan married nodegree re74 re75 re78
PSID6 0 37 9 1 0 1 1 13685 12756 17833
> lalonde[212,]
treat age educ black hispan married nodegree re74 re75 re78
PSID27 0 36 9 1 0 1 1 13256 8457 0
> lalonde[183,]
treat age educ black hispan married nodegree re74 re75 re78
NSW183 1 35 9 1 0 1 1 13602 13831 12804
> indiv1
treat age educ black hispan married nodegree re74 re75 re78
NSW1 1 37 11 1 0 1 1 0 0 9930
> # if we took re74 re75 into account we wouldn't have any control matches for indiv1
# aside, opmatch demo
> # optmatch (Ben Hansen) could also find matches to person 1
>
# I just used age and education for a demo; if I had further specified the dichotomous variables
I would have an alternative to the results above
> m1 = matchit(treat ~ age + educ, data = lalonde, method = "optimal", ratio = 2)
Loading required package: optmatch
m1.data = match.data(m1)
# so I created 185 subclasses each with 3 members, (2 control, 1 treatment), fun
> names(m1.data)
[1] "treat" "age" "educ" "black" "hispan" "married" "nodegree" "re74" "re75" "re78"
[11] "distance" "weights" "subclass"
# pick off subclass 1; two control members (not the same as above as i didn't use the other vars)
> s1 = subset(m1.data, subset = subclass < 2, select = c(age, educ, distance, subclass))
> dim(s1)
[1] 3 4
> s1
age educ distance subclass
NSW1 37 11 0.2537 1
PSID9 38 9 0.2467 1
PSID400 37 8 0.2499 1
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
> # Lab script 3.2
# now do the logistic regression that computes propensity scores (matching packages will do this for you)
> glm.lalonde = glm( treat ~ age + educ + black + hispan + married + nodegree + re74 + re75, data = lalonde, family = binomial)
> summary(glm.lalonde)
Call:
glm(formula = treat ~ age + educ + black + hispan + married +
nodegree + re74 + re75, family = binomial, data = lalonde)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.764 -0.474 -0.286 0.751 2.717
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.73e+00 1.02e+00 -4.65 3.3e-06 ***
age 1.58e-02 1.36e-02 1.16 0.2452
educ 1.61e-01 6.51e-02 2.48 0.0133 *
black 3.07e+00 2.87e-01 10.70 < 2e-16 ***
hispan 9.84e-01 4.26e-01 2.31 0.0208 *
married -8.32e-01 2.90e-01 -2.87 0.0042 **
nodegree 7.07e-01 3.38e-01 2.09 0.0362 *
re74 -7.18e-05 2.87e-05 -2.50 0.0125 *
re75 5.34e-05 4.63e-05 1.15 0.2488
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 751.49 on 613 degrees of freedom
Residual deviance: 487.84 on 605 degrees of freedom
AIC: 505.8
Number of Fisher Scoring iterations: 5
# best case is counter-intuitive; you really want not-wonderful predictions (analog to Rsq) so that
# you have both treat=1 and treat=0 cases at levels of propen
> propen = fitted(glm.lalonde) # now we have the propensity scores, Lab script calls these propScores
> tapply(propen, treat, quantile) # look at overlap via 5-number summary (or side-by-side boxplots)
not good overlap, as noted in class handout snippet
$`0`
0% 25% 50% 75% 100%
0.00908 0.03888 0.07585 0.19514 0.78917
$`1`
0% 25% 50% 75% 100%
0.02495 0.52646 0.65368 0.72660 0.85315
> # we are fitting prob(treat = 1) fits for those in treatment group
> # will be larger, we need some overlap for matching purposes
> boxplot(propen ~ treat) #gives side-by-side boxplots, you can add labels, not pretty compared to class exs
> tapply(propen, treat, stem)
The decimal point is 1 digit(s) to the left of the | # treat = 0
0 | 11111111111112222222222222222222222222222222222222222222233333333333+57
0 | 55555555555555555555555555555556666666666666666666666666777777777777+39
1 | 000000000000000011111111111111112222222222223334444
1 | 56677888888889
2 | 0000011122244
2 | 55667788
3 | 0
3 | 6778899
4 | 1233444
4 | 568999
5 | 1144
5 | 5679
6 | 11112234444
6 | 555556667778888899999999
7 | 111112222444
7 | 55555555689
The decimal point is 1 digit(s) to the left of the | # treat = 1
0 | 244
0 | 6788899
1 | 000122334444
1 | 55
2 | 012
2 | 579
3 |
3 | 8899
4 | 0134
4 | 5559
5 | 1112334
5 | 6666677778899999
6 | 00111122222334444
6 | 55555555555666666778888888899999
7 | 000000001111222222222223333344444
7 | 555555566666777777777778888888999
8 | 1122
8 | 5
$`0`
NULL
$`1`
NULL
> quantile(propen) # overall distrib
0% 25% 50% 75% 100%
0.00908 0.04854 0.12068 0.63872 0.85315
# or you can have MatchIt compute propensity scores for you and do a simple match (as in the Aspirin Study (Love))
######subclassification (by hand) #################
> # a common use of the propensity scores (backed by theory, class handout week8))
> # is to stratify by quintiles
> # the simple-minded way I do it is to use "cut", Lab script is fancier programming
> ?cut # this is a simple function to create bins
> k = 1:4
> quantile(propen, k/5) # defines 5 bins (equal sample size) for the 614 propensity scores
20% 40% 60% 80%
0.04015 0.08721 0.26978 0.67085
> propbin = cut(propen, c(0, .04015,.08721,.26978,.67085,1))
> propen[1:10] # show propen values that get put into bins for first 10 cases
NSW1 NSW2 NSW3 NSW4 NSW5 NSW6 NSW7 NSW8 NSW9 NSW10
0.63877 0.22463 0.67824 0.77632 0.70164 0.69907 0.65368 0.78972 0.77984 0.04292
> propbin[1:10] # here are the bins for those values
[1] (0.27,0.671] (0.0872,0.27] (0.671,1] (0.671,1]
[5] (0.671,1] (0.671,1] (0.27,0.671] (0.671,1]
[9] (0.671,1] (0.0401,0.0872]
Levels: (0,0.0401] (0.0401,0.0872] (0.0872,0.27] (0.27,0.671] (0.671,1]
> # levels are the 5 quintiles
> table(propbin, treat) # either way you display it, we do not have good overlap in the bottom
two quintiles, lower estimated probability for being in treatment
for treatment cases (see your side-by-side boxplot)
treat
propbin 0 1
(0,0.0401] 122 1
(0.0401,0.0872] 116 7
(0.0872,0.27] 101 21
(0.27,0.671] 53 71
(0.671,1] 37 85
> tapply(propbin, treat,table) # show same thing sideways
$`0`
(0,0.0401] (0.0401,0.0872] (0.0872,0.27] (0.27,0.671] (0.671,1]
122 116 101 53 37
$`1`
(0,0.0401] (0.0401,0.0872] (0.0872,0.27] (0.27,0.671] (0.671,1]
1 7 21 71 85
> table(treat)
treat
0 1
429 185
> #we have (real) overlap only in the upper 3 quintiles of the propensity score
##############################
## an important diagnostic step we skipped over is to check balance (as we did in class)
## on each of the confounding vars (i.e. covariates) within each subclassification, c.f Rosebaum-Rubin ex
# a quick look in not overwhelmingly impressive, but better than nothing
# for the two groups
> tapply(propen, treat, mean) # of course these should differ
0 1
0.1822248 0.5774355
> tapply(age, treat, mean)
0 1
28.03030 25.81622
> tapply(educ, treat, mean)
0 1
10.23543 10.34595
> tapply(married, treat, mean) # big diff here
0 1
0.5128205 0.1891892
> tapply(nodegree, treat, mean)
0 1
0.5967366 0.7081081
> tapply(re74, treat, mean) # big diff
0 1
5619.237 2095.574
> tapply(re75, treat, mean)
0 1
2466.484 1532.055
# after this initial quintile subclassification on propensity score, look at balance within quintile
# compare to overall above
# keep in mind the top 3 quintiles is where we have data for treat-control comparison
# first look at propen score
> tapply(propen, list(propbin, treat), mean)
0 1
(0,0.0401] 0.02603935 0.02495179
(0.0401,0.0872] 0.06106601 0.06589413
(0.0872,0.27] 0.13761477 0.14830101
(0.27,0.671] 0.51792739 0.57075500
(0.671,1] 0.71796548 0.73766399
> tapply(propen, list(propbin, treat), median)
0 1
(0,0.0401] 0.02592168 0.02495179
(0.0401,0.0872] 0.05951047 0.07382733
(0.0872,0.27] 0.11515593 0.13557377
(0.27,0.671] 0.53948683 0.59404299
(0.671,1] 0.71467967 0.73802599
> names(lalonde)
[1] "treat" "age" "educ" "black" "hispan" "married" "nodegree" "re74" "re75"
[10] "re78"
> tapply(age, list(propbin, treat), median)
0 1
(0,0.0401] 29 27
(0.0401,0.0872] 26 23
(0.0872,0.27] 20 23
(0.27,0.671] 24 25
(0.671,1] 19 25
> tapply(age, list(propbin, treat), mean) # not better than overall
0 1
(0,0.0401] 32.09016 27.00000
(0.0401,0.0872] 28.93103 24.00000
(0.0872,0.27] 23.92079 24.95238
(0.27,0.671] 28.52830 26.18310
(0.671,1] 22.32432 25.85882
> tapply(educ, list(propbin, treat), mean) # not bad, but not better than overall
0 1
(0,0.0401] 9.877049 13.00000
(0.0401,0.0872] 10.086207 10.71429
(0.0872,0.27] 10.910891 10.52381
(0.27,0.671] 9.886792 10.05634
(0.671,1] 10.540541 10.48235
> tapply(married, list(propbin, treat), mean) # alot better than overall esp in top 3
0 1
(0,0.0401] 0.9754098 1.00000000
(0.0401,0.0872] 0.5603448 0.28571429
(0.0872,0.27] 0.1188119 0.19047619
(0.27,0.671] 0.4528302 0.38028169
(0.671,1] 0.0000000 0.01176471
> tapply(nodegree, list(propbin, treat), mean) # better than overall esp in top 3
0 1
(0,0.0401] 0.5245902 0.0000000
(0.0401,0.0872] 0.5775862 0.2857143
(0.0872,0.27] 0.6534653 0.7142857
(0.27,0.671] 0.5660377 0.5915493
(0.671,1] 0.7837838 0.8470588
> tapply(re74, list(propbin, treat), mean) # better than overall
0 1
(0,0.0401] 11879.0344 9381.5660
(0.0401,0.0872] 4716.3099 1490.1237
(0.0872,0.27] 1550.6141 3177.4622
(0.27,0.671] 4319.4811 3540.6004
(0.671,1] 777.6701 585.4043
> tapply(re75, list(propbin, treat), mean) # better than overall
0 1
(0,0.0401] 3762.2162 853.7225
(0.0401,0.0872] 2617.8838 1500.3862
(0.0872,0.27] 1254.2185 2009.7170
(0.27,0.671] 2504.0532 1858.2195
(0.671,1] 974.7579 1152.1901
##############################
> # lab text does a more sophisticated quintile function (see below), more exact as we see
> # IMPLEMENT LAB TEXT 3.3
> quintilePropScores = quantile( propen, probs = seq( from = 0, to = 1, by = .2) )
> quintilePropScores
0% 20% 40% 60% 80% 100%
0.00908 0.04015 0.08721 0.26978 0.67085 0.85315
# has bin breaks to greater signif digits, that makes a small difference it turns out
> lalonde$PropScores = propen # make var names consistent
#################################
> table(propbin, treat)
treat
propbin 0 1
(0,0.0401] 122 1
(0.0401,0.0872] 116 7
(0.0872,0.27] 101 21
(0.27,0.671] 53 71
(0.671,1] 37 85
######### now start looking at outcome
> tapply(re78, list(propbin, treat),mean) # here are the mean diffs in re78 the outcome
stratified by propensity quintile
# direction of mean diffs favors treatment, job training
0 1
(0,0.0401] 10467 0
(0.0401,0.0872] 5797 7919
(0.0872,0.27] 6043 9211
(0.27,0.671] 4977 5819
(0.671,1] 4666 6030
# contrast that with the comparison ignoring any attempt to balance, in the other direction, but not significant
> tapply(re78, treat, mean)
0 1
6984.170 6349.144
> t.test(re78 ~ treat)
Welch Two Sample t-test
data: re78 by treat
t = 0.9377, df = 326.412, p-value = 0.3491
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-697.192 1967.244
sample estimates:
mean in group 0 mean in group 1
6984.170 6349.144
> lm(re78 ~ treat, data = lalonde) # another way to look
Call:
lm(formula = re78 ~ treat, data = lalonde)
Coefficients:
(Intercept) treat
6984 -635
##### do t-tests by subclassification (strata)
> tapply(re78, list(propbin, treat),t.test) # can't do a full set of t-tests
Error in t.test.default(X[[6L]], ...) : not enough 'x' observations
> # so for the 3 upper quintiles is the mean difference significant?
attributes for case 1, bin and propensity score
> propbin[1]
[1] (0.27,0.671]
Levels: (0,0.0401] (0.0401,0.0872] (0.0872,0.27] (0.27,0.671] (0.671,1]
> propen[1]
NSW1
0.6388
> bins = c("(0,0.0401]", "(0.0401,0.0872]", "(0.0872,0.27]", "(0.27,0.671]", "(0.671,1]")
> # i.e. make a string out of the bin designations, then pull-off entries
#####actually you can just do t.test(re78[propbin == propbin[5]] ~ treat[propbin == propbin[5]]) # t-test for quintile 5
> t.test(re78[propbin == bins[5]] ~ treat[propbin == bins[5]]) # t-test for quintile 5
Welch Two Sample t-test
data: re78[propbin == bins[5]] by treat[propbin == bins[5]]
t = -0.9844, df = 100.7, p-value = 0.3273
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-4113 1385
sample estimates:
mean in group 0 mean in group 1
4666 6030
# since we are doing 3 of these best practice would be to Bonferroni
and do each at Type 1 .017 instead of .05. Also Karen did these one-sided which you might justify,
won't make a difference here; see Lab text version below
> t.test(re78[propbin == bins[4]] ~ treat[propbin == bins[4]])
Welch Two Sample t-test
data: re78[propbin == bins[4]] by treat[propbin == bins[4]]
t = -0.7516, df = 108.7, p-value = 0.4539
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3061 1378
sample estimates:
mean in group 0 mean in group 1
4977 5819
> t.test(re78[propbin == bins[3]] ~ treat[propbin == bins[3]])
Welch Two Sample t-test
data: re78[propbin == bins[3]] by treat[propbin == bins[3]]
t = -1.536, df = 23.57, p-value = 0.1379
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-7428 1093
sample estimates:
mean in group 0 mean in group 1
6043 9211
> # even one-sided versions of these tests are not going to be signif; should also Bonferroni the Type I error
> t.test(re78[propbin == bins[2]] ~ treat[propbin == bins[2]]) # can also try bin 2 but really don't have enough data
Welch Two Sample t-test
data: re78[propbin == bins[2]] by treat[propbin == bins[2]]
t = -1.083, df = 7.38, p-value = 0.3128
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-6708 2462
sample estimates:
mean in group 0 mean in group 1
5797 7919
## so in summary we can't find any evidence for the effectiveness of job training looking at each of the subclasses
##### lmer, a better way to do the t-tests ########################
# Refer to week 4, lab2, TH1 for our background on lmer
# What we are going to do here is carry out the t-test within each of the 5 subclassifications (level-1 model)
# and then allow lmer to 'average' (weighting by the precision basically) those 5 comparisons. There are no
# conditioning variables in the level 2 model so the fixed effects are just the 'averages' of the level 1 effects
# being careful I added propen and the bin classification to the data set
> detach(lalonde)
> lalonde$propen = propen
> lalonde$bins = propbin
> lalonde[1:10,]
treat age educ black hispan married nodegree re74 re75 re78 propen bins
NSW1 1 37 11 1 0 1 1 0 0 9930.0460 0.63876993 (0.27,0.671]
NSW2 1 22 9 0 1 0 1 0 0 3595.8940 0.22463424 (0.0872,0.27]
NSW3 1 30 12 1 0 0 0 0 0 24909.4500 0.67824388 (0.671,1]
NSW4 1 27 11 1 0 0 1 0 0 7506.1460 0.77632408 (0.671,1]
NSW5 1 33 8 1 0 0 1 0 0 289.7899 0.70163874 (0.671,1]
NSW6 1 22 9 1 0 0 1 0 0 4056.4940 0.69906990 (0.671,1]
NSW7 1 23 12 1 0 0 0 0 0 0.0000 0.65368426 (0.27,0.671]
NSW8 1 32 11 1 0 0 1 0 0 8472.1580 0.78972311 (0.671,1]
NSW9 1 22 16 1 0 0 0 0 0 2164.0220 0.77983825 (0.671,1]
NSW10 1 33 12 0 0 1 0 0 0 12418.0700 0.04292461 (0.0401,0.0872]
> str(lalonde)
'data.frame': 614 obs. of 12 variables:
$ treat : int 1 1 1 1 1 1 1 1 1 1 ...
$ age : int 37 22 30 27 33 22 23 32 22 33 ...
$ educ : int 11 9 12 11 8 9 12 11 16 12 ...
$ black : int 1 0 1 1 1 1 1 1 1 0 ...
$ hispan : int 0 1 0 0 0 0 0 0 0 0 ...
$ married : int 1 0 0 0 0 0 0 0 0 1 ...
$ nodegree: int 1 1 0 1 1 1 0 1 0 0 ...
$ re74 : num 0 0 0 0 0 0 0 0 0 0 ...
$ re75 : num 0 0 0 0 0 0 0 0 0 0 ...
$ re78 : num 9930 3596 24909 7506 290 ...
$ propen : num 0.639 0.225 0.678 0.776 0.702 ...
$ bins : Factor w/ 5 levels "(0,0.0401]","(0.0401,0.0872]",..: 4 3 5 5 5 5 4 5 5 2 ...
> propen.lmer = lmer(re78 ~ treat + (1 + treat|bins), data = lalonde)
> summary(propen.lmer)
Linear mixed model fit by REML
Formula: re78 ~ treat + (1 + treat | bins)
Data: lalonde
AIC BIC logLik deviance REMLdev
12649 12676 -6319 12668 12637
Random effects:
Groups Name Variance Std.Dev. Corr
bins (Intercept) 5208931 2282.3
treat 2069968 1438.7 -1.000
Residual 52597983 7252.4
Number of obs: 614, groups: bins, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) 6434.3 1090.2 5.902
treat 383.9 951.5 0.403
Correlation of Fixed Effects:
(Intr)
treat -0.796
# so here we have an overall estimate of the effect of the treat on re78 of positive $384, but
# far from significant. Much smaller point estimate than in some of the individual strata
################################
##Ancova with propensity score alternative ##############################
#One not terrible, but not great, use of propensity scores which you likely will encounter is to use
# the propensity score as a covariate.
# Not great, but not the worst analysis possible. An analogy could be made to covariance vs blocking
# where here we use regression smoothing rather than blocks (strat). Or to nonrandom assignment
# (week 5) where at each level of propensity, use a (biased) coin flip to assign to treat or control.
# One reason the ancova is not great is that the OLS regression will give most influence to extremes
# of propen, where we have no overlap.
> propen_ancova = lm(re78 ~ treat + propen)
> summary(propen_ancova)
Call:
lm(formula = re78 ~ treat + propen)
Residuals:
Min 1Q Median 3Q Max
-8925 -5632 -2148 3993 54953
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7833.8 430.7 18.188 < 2e-16 ***
treat 1207.6 834.1 1.448 0.148200
propen -4662.4 1319.4 -3.534 0.000441 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7402 on 611 degrees of freedom
Multiple R-squared: 0.02152, Adjusted R-squared: 0.01832
F-statistic: 6.72 on 2 and 611 DF, p-value: 0.001298
# so bigger point estimate (in the same direction) than lmer on strata, but still not significant.
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
### Karen's original functions (for analysis the lmer full is better; Karen followed the t-test in Gary King's tutorial)
> # lab text does a more sophisticated quintile function
> # IMPLEMENT LAB TEXT 3.3
> quintilePropScores = quantile( propen, probs = seq( from = 0, to = 1, by = .2) )
> quintilePropScores[1:10]
0% 20% 40% 60% 80% 100%
0.00908 0.04015 0.08721 0.26978 0.67085 0.85315 NA NA NA NA
> quintilePropScores
0% 20% 40% 60% 80% 100%
0.00908 0.04015 0.08721 0.26978 0.67085 0.85315
> lalonde$PropScores = propen
> lalonde$PropScores[1:10]
[1] 0.63877 0.22463 0.67824 0.77632 0.70164 0.69907 0.65368 0.78972 0.77984 0.04292
> # here's a pretty fancy do-loop version from Lab script
> treated = list(5)
> controls = list(5)
> for ( i in 1:5 ) {
+ treated[[i]] = subset( x = lalonde,
+ subset = PropScores >= quintilePropScores[i] &
+ PropScores < quintilePropScores[i+1] &
+ treat == 1)
+ controls[[i]] = subset( x = lalonde,
+ subset = PropScores >= quintilePropScores[i] &
+ PropScores < quintilePropScores[i+1] &
+ treat == 0)
+ }
> for ( i in 1:5 ) {
+ t.test( treated[[i]]$re78,
+ controls[[i]]$re78,
+ alternative = "greater" )
+ }
Error in t.test.default(treated[[i]]$re78, controls[[i]]$re78, alternative = "greater") :
not enough 'x' observations
# so I change the loop in Lab text to exclude testing bin 1
> for ( i in 2:5 ) {
+ t.test( treated[[i]]$re78,
+ controls[[i]]$re78,
+ alternative = "greater" )
+ }
> ptests = for ( i in 2:5 ) {
+ t.test( treated[[i]]$re78,
+ controls[[i]]$re78,
+ alternative = "greater" )
+ }
> ptests
Welch Two Sample t-test
data: treated[[i]]$re78 and controls[[i]]$re78
t = 1.120, df = 108.3, p-value = 0.1326
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
-735.8 Inf
sample estimates:
mean of x mean of y
6027 4499
> t.test( treated[[5]]$re78,
+ controls[[5]]$re78,
+ alternative = "greater" )
Welch Two Sample t-test
data: treated[[5]]$re78 and controls[[5]]$re78
t = 1.120, df = 108.3, p-value = 0.1326
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
-735.8 Inf
sample estimates:
mean of x mean of y
6027 4499
> t.test( treated[[4]]$re78,
+ controls[[4]]$re78,
+ alternative = "greater" )
Welch Two Sample t-test
data: treated[[4]]$re78 and controls[[4]]$re78
t = 0.6149, df = 103.4, p-value = 0.27
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
-1193 Inf
sample estimates:
mean of x mean of y
5819 5117
> t.test( treated[[3]]$re78,
+ controls[[3]]$re78,
+ alternative = "greater" )
Welch Two Sample t-test
data: treated[[3]]$re78 and controls[[3]]$re78
t = 1.536, df = 23.57, p-value = 0.06895
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
-363.7 Inf
sample estimates:
mean of x mean of y
9211 6043
> length(treated[[3]]$re78)
[1] 21
> mean(treated[[3]]$re78)
[1] 9211
> length(controls[[3]]$re78)
[1] 101
> mean(controls[[3]]$re78)
[1] 6043
> # those match my bins
> length(treated[[4]]$re78)
[1] 71
> mean(treated[[4]]$re78)
[1] 5819
> length(controls[[4]]$re78)
[1] 51
> mean(controls[[4]]$re78)
[1] 5117
> # our different bins procedure are a little discrepant and that matters, 2 controls
##########################################################
####Doing this by hand vs letting MatchIt do all
##
##if you type at the R-prompt 'demo(subclass)' you will be taken
##on a tour of MatchIt output. This demo uses a subset of our matching
##variables, hardwired to 6 subclasses, but these fail to be of equal
##size (even close) like they are supposed to be. Shows you the balance statistics
##and some charming plots, (e.g. absolute standardized diffs)if you let it run.
##For pedagogy I chose to do all this step by step.
# The matchit JSS paper, linked in week 8, provides better explanation of using the subclass option
# than the main matchit docs
#####################*****end Lab script part 3 ***********###################################