Illustration of Matching for birthwt data (a subset of variables, as noted)
Define the matrix of interest in terms of the call for bwt.lrr (class of Mar 17)
For that, Call: glm(formula = smoke ~ age + lwt + factor(race) + pt.1 + ui, family = "binomial", data = birthwt) …. But w/ pt.1 defined as binary
names(bwt.sub7) "age" "lwt" "race" "smoke" "ui" "bwt" "pt.1" (n=189)
bwt.match7=Matchby(Y=bwt/(28.35*16),Tr=smoke,X=bwt.sub7[,c(1,2,5,7)],by=list(race)) # note the object name here… and see summary below.
1 of 3 groups 2 of 3 groups
3 of 3 groups
> summary(bwt.match7)
Estimate... -1.0055
SE...... 0.24684
T-stat..... -4.0733
p.val...... 4.6344e-05
Original number of observations...... 189
Original number of treated obs...... 74
Matched number of observations...... 74
Matched number of observations (unweighted). 74
Ø table(smoke) shows that we seek 74 matches from 0’s for 1’s
Ø 0 1
115 74
bwt.match7 $est = -1.01
$se.standard = 0.247
$weights all ones… ignore
$orig.nobs = 189
$orig.wnobs = 189
$nobs 74 (order of these few summary lines not same as in original output… for when you try this)
$wnobs 74 #these two lines simply confirm the 1:1 matching for ‘full matching’
$estimand "ATT" #will discuss $version "Matchby"
$index.treated # i.e., the subscripts for the smokers (treated)
[1] 3 4 5 9 10 15 16 18 20 27 32 36 37 38 39 40 45 46 53 54 71 75 78 93 94 95 98 99 102 109 113 115 128 140 141 142 145 146 148 153 154 158 161 166 169 176 177 178 179 183 185 189 29 33 41 79 133 157 165 171 172 173 50 57 65 68 72 73 86 131 156 160 184 187
$index.control $ i.e., subscripts, matches for smokers, aligned as per indices above.
So, smoker (no. 3) is matched w/ non-smoker (no. 152), etc.
[1] 152 132 132 126 119 152 152 119 107 170 101 96 81 119 23 126 22 22 124 107 104 117 127 42 69 117 121 121 42 123 96 108 123 44 107 69 101 121 127 132 127 69 170 69 127 127 124 170 104 80 107 81 35 138 110 180 106 180 138 106 180 35 13 19 181 122 58 62 61 14 66 137 167 13
So, based on these results, we can use granova.ds to show results for the matched pairs: granova.ds((cbind(bwt[bwt.match7$index.treated],bwt[bwt.match7$index.control] )/(28.35*16))[,2:1],xlab="Birthwts in lbs. for babies born to non-smokers",ylab='Birthwts in lbs for babies born to smokers') # take some time to parse this line. Summary Stats
n 74.000mean(x) 7.116
mean(y) 6.111
mean(D=x-y) 1.005
SD(D) 2.138
ES(D) 0.470
r(x,y) 0.183 / r(x+y,d) 0.243
LL 95%CI 0.510
UL 95%CI 1.501 $compare CI w/
that for loess plot
t(D-bar) 4.046
df.t 73.000
pval.t 0.000
For comparison, recall the loess run: (where STRATIFICATION was used, not matching)
loess.bwt2 = loess.psa(bwt/(28.35*16), treatment=smoke,fitted(bwt.lrr),int=7)
for which the numerical results were: loess.bwt2
$ATE = -0.832
$se.wtd = 0.308
$CI95 -1.448 -0.216 #compare w/ result for Full Matching above.
(reverse signs for CI here to compare w/ those above)
$summary.strata
counts.0 counts.1 means.0 means.1 diff.means
1 24 3 6.42 6.44 0.0130
2 22 5 6.43 6.24 -0.1896
3 19 8 6.60 6.00 -0.6011
4 17 10 6.90 5.99 -0.9064
5 15 12 7.12 6.22 -0.8913
6 12 15 7.22 6.32 -0.9068
7 6 21 7.36 5.86 -1.4980
Next we focus on assessing balance (partially)… using a graphic.
cv.trans.bwt=cv.trans.psa(bwt.sub7[,-6],fcol=3)
since third variable is race, and that is a factor, we need to transform it to dummies, while leaving the rest unchanged. The command above gives us what we need, as:
cv.trans.bwt$covariates.transf
head(cv.trans.bwt$covariates.transf)
age lwt smoke ui pt.1 race_1 race_2 race_3
85 19 182 0 1 0 0 1 0
86 33 155 0 0 0 0 0 1
87 20 105 1 0 0 1 0 0
88 21 108 1 1 0 1 0 0
89 18 107 1 1 0 1 0 0
91 21 124 0 0 0 0 0 1
------
Use the cut command here:
cut.fit.bwt=cut2(fitted(bwt.lrr),quantile(fitted(bwt.lrr),prob=seq(.01,.99,len=7))) #parse, as usual...[NB: cut2 from Hmisc library!]
so we have table(cut.fit.bwt)
cut.fit.bwt
[0.0525,0.0732) [0.0732,0.1521) [0.1521,0.2909) [0.2909,0.4121) [0.4121,0.4936) [0.4936,0.5865) [0.5865,0.8489) [0.8489,0.8838]
2 31 31 30 31 31 31
Next, we assess COVARIATE BALANCE
cv.bal.bwt7=cv.bal.psa(cv.trans.bwt$covariates.transf,smoke,fitted(bwt.lrr), strata=cut.fit.bwt)
which yields a graphic and numerical summaries:
cv.bal.bwt7#=cv.bal.psa(cv.trans.bwt$covariates.transf,smoke,fitted(bwt.lrr),strata=cut.fit.bwt) $original.strata
[0.0525,0.0732) [0.0732,0.1521) [0.1521,0.2909) [0.2909,0.4121) [0.4121,0.4936) [0.4936,0.5865) [0.5865,0.8489) [0.8489,0.8838]
0 1 28 23 21 19 16 7 0
1 1 3 8 9 12 15 24 2
cv.bal.bwt7
$original.strata
[0.0525,0.0732) [0.0732,0.1521) [0.1521,0.2909) [0.2909,0.4121) [0.4121,0.4936) [0.4936,0.5865) [0.5865,0.8489) [0.8489,0.8838]
0 1 28 23 21 19 16 7 0
1 1 3 8 9 12 15 24 2
$strata.used
a b c d e f
0 28 23 21 19 16 7
1 3 8 9 12 15 24
Strata.Size 31 31 30 31 31 31
$mean.diff.strata.wtd
a b c d e f
race_3 0.00 0.05 0.02 0.03 0.00 0.00
race_1 0.00 0.01 0.00 0.02 0.01 0.04
pt.1 0.00 0.01 0.02 0.03 0.01 0.02
ui 0.01 0.01 0.02 0.00 0.01 0.02
age 0.29 0.35 0.03 0.22 0.18 0.02
lwt 0.88 2.30 0.90 0.73 1.45 3.43
race_2 0.00 0.07 0.02 0.04 0.01 0.04
$mean.diff.unadj
race_3 race_1 pt.1 ui age lwt race_2
0.32 0.32 0.14 0.05 0.48 2.76 0.00
$effect.sizes
stES_unadj stES_adj a b c d e f sd.ESs
race_3 0.72 0.10 0.00 0.74 0.25 0.38 0.00 0.00 0.30
race_1 0.67 0.07 0.00 0.18 0.00 0.20 0.14 0.51 0.19
pt.1 0.37 0.01 0.00 0.12 0.30 0.44 0.18 0.30 0.16
ui 0.13 0.03 0.10 0.09 0.35 0.06 0.17 0.41 0.15
age 0.09 0.12 0.33 0.40 0.04 0.25 0.20 0.02 0.15
lwt 0.09 0.07 0.17 0.44 0.18 0.14 0.28 0.65 0.20
race_2 0.01 0.04 0.00 1.19 0.32 0.76 0.19 0.71 0.44
$treatment.levels
[1] 0 1
$effects.strata.treatment ---- omitted here.
Note that the adjusted (blue) and unadjusted (red) ESs are given in the graphic below, summarized for each covariate. That the blue adjusted ESs never exceed .15 s.d. units. It means that, taken one-by-one the covariates are relatively well balanced using seven strata after propensity score adjustment. This is for the loess run; similar assessments of balance would be desirable for the matched data. See the functions in Matching package