CausalInferenceWithR
MiLin
4/7/2022
1 bcf
这个包是干什么的： Causal Inference for a Binary Treatment and Continuous Outcome using Bayesian Causal Forests
y :response variable z：Treatment variable x_control ：Design matrix for the “prognostic” function mu(x) x_moderate：Design matrix for the covariate-dependent treatment effects tau(x)
library(bcf)## Warning: package 'bcf' was built under R version 4.1.2## Loading required package: Rcpp# data generating process
p = 3 #two control variables and one moderator
n = 250
#
set.seed(1)
x = matrix(rnorm(n*p), nrow=n)
# create targeted selection
q = -1*(x[,1]>(x[,2])) + 1*(x[,1]<(x[,2]))
# generate treatment variable
pi = pnorm(q)
z = rbinom(n,1,pi)
# tau is the true (homogeneous) treatment effect
tau = (0.5*(x[,3] > -3/4) + 0.25*(x[,3] > 0) + 0.25*(x[,3]>3/4))
# generate the response using q, tau and z
mu = (q + tau*z)
# set the noise level relative to the expected mean function of Y
sigma = diff(range(q + tau*pi))/8
# draw the response variable with additive error
y = mu + sigma*rnorm(n)
# If you didn't know pi, you would estimate it here
pihat = pnorm(q)
bcf_fit = bcf(y, z, x, x, pihat, nburn=2000, nsim=2000)## Using 4 control variables.
## Using 3 potential effect moderators.
##
## Beginning MCMC:
## iteration: 0 sigma: 1
## iteration: 100 sigma: 0.311768
## iteration: 200 sigma: 0.306581
## iteration: 300 sigma: 0.273944
## iteration: 400 sigma: 0.252862
## iteration: 500 sigma: 0.28624
## iteration: 600 sigma: 0.284792
## iteration: 700 sigma: 0.270534
## iteration: 800 sigma: 0.281106
## iteration: 900 sigma: 0.266863
## iteration: 1000 sigma: 0.284123
## iteration: 1100 sigma: 0.267538
## iteration: 1200 sigma: 0.257672
## iteration: 1300 sigma: 0.279696
## iteration: 1400 sigma: 0.280306
## iteration: 1500 sigma: 0.293692
## iteration: 1600 sigma: 0.284637
## iteration: 1700 sigma: 0.267644
## iteration: 1800 sigma: 0.269394
## iteration: 1900 sigma: 0.268582
## iteration: 2000 sigma: 0.297118
## iteration: 2100 sigma: 0.293801
## iteration: 2200 sigma: 0.254708
## iteration: 2300 sigma: 0.284539
## iteration: 2400 sigma: 0.291441
## iteration: 2500 sigma: 0.25212
## iteration: 2600 sigma: 0.266714
## iteration: 2700 sigma: 0.289271
## iteration: 2800 sigma: 0.287449
## iteration: 2900 sigma: 0.264065
## iteration: 3000 sigma: 0.294549
## iteration: 3100 sigma: 0.277946
## iteration: 3200 sigma: 0.27893
## iteration: 3300 sigma: 0.27674
## iteration: 3400 sigma: 0.266013
## iteration: 3500 sigma: 0.251627
## iteration: 3600 sigma: 0.261951
## iteration: 3700 sigma: 0.268848
## iteration: 3800 sigma: 0.296868
## iteration: 3900 sigma: 0.248751
## time for loop: 7# Get posterior of treatment effects
tau_post = bcf_fit$tau
tauhat = colMeans(tau_post)
plot(tau, tauhat); abline(0,1)
2 causalweight
y : Dependent variable d : Treatment ,must be binary x : Confounders of the treatment and outcome, must not contain missings.
library(causalweight)## Loading required package: rangern = 10000
set.seed(100)
x = rnorm(n)
set.seed(101)
d = (0.25 * x + rnorm(n) > 0) * 1
set.seed(102)
y = 0.5 * d + 0.25 * x + rnorm(n)
output <-
treatweight(
y = y,
d = d,
x = x,
trim = 0.05,
ATET = FALSE,
logit = TRUE,
boot = 19
)
cat("ATE:",
round(c(output$effect), 3),
", standard error: ",
round(c(output$se), 3),
", p-value: ",
round(c(output$pval), 3))## ATE: 0.488 , standard error: 0.022 , p-value: 03 CausalImpact
Inferring Causal Effects using Bayesian Structural Time-Series Models
library(CausalImpact)## Loading required package: bsts## Loading required package: BoomSpikeSlab## Loading required package: Boom## Loading required package: MASS##
## Attaching package: 'Boom'## The following object is masked from 'package:stats':
##
## rWishart##
## Attaching package: 'BoomSpikeSlab'## The following object is masked from 'package:stats':
##
## knots## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: xts##
## Attaching package: 'bsts'## The following object is masked from 'package:BoomSpikeSlab':
##
## SuggestBurnset.seed(1)
x1 <- 100 + arima.sim(model = list(ar = 0.999), n = 52)
y <- 1.2 * x1 + rnorm(52)
y[41:52] <- y[41:52] + 10
data <- cbind(y, x1)
pre.period <- c(1, 40)
post.period <- c(41, 52)
impact <- CausalImpact(data, pre.period, post.period)
# Print and plot results
summary(impact)## Posterior inference {CausalImpact}
##
## Average Cumulative
## Actual 112 1340
## Prediction (s.d.) 102 (0.4) 1223 (4.8)
## 95% CI [101, 103] [1214, 1233]
##
## Absolute effect (s.d.) 9.8 (0.4) 117.1 (4.8)
## 95% CI [8.9, 11] [107.4, 126]
##
## Relative effect (s.d.) 9.6% (0.39%) 9.6% (0.39%)
## 95% CI [8.8%, 10%] [8.8%, 10%]
##
## Posterior tail-area probability p: 0.001
## Posterior prob. of a causal effect: 99.9%
##
## For more details, type: summary(impact, "report")summary(impact, "report")## Analysis report {CausalImpact}
##
##
## During the post-intervention period, the response variable had an average value of approx. 111.70. By contrast, in the absence of an intervention, we would have expected an average response of 101.94. The 95% interval of this counterfactual prediction is [101.16, 102.75]. Subtracting this prediction from the observed response yields an estimate of the causal effect the intervention had on the response variable. This effect is 9.76 with a 95% interval of [8.95, 10.53]. For a discussion of the significance of this effect, see below.
##
## Summing up the individual data points during the post-intervention period (which can only sometimes be meaningfully interpreted), the response variable had an overall value of 1.34K. By contrast, had the intervention not taken place, we would have expected a sum of 1.22K. The 95% interval of this prediction is [1.21K, 1.23K].
##
## The above results are given in terms of absolute numbers. In relative terms, the response variable showed an increase of +10%. The 95% interval of this percentage is [+9%, +10%].
##
## This means that the positive effect observed during the intervention period is statistically significant and unlikely to be due to random fluctuations. It should be noted, however, that the question of whether this increase also bears substantive significance can only be answered by comparing the absolute effect (9.76) to the original goal of the underlying intervention.
##
## The probability of obtaining this effect by chance is very small (Bayesian one-sided tail-area probability p = 0.001). This means the causal effect can be considered statistically significant.plot(impact)
plot(impact, "original")
plot(impact$model$bsts.model, "coefficients")
# For further output, type:
names(impact)## [1] "series" "summary" "report" "model"cit
Causal Inference Test
Causal Inference Test for a Binary Outcome
library(cit)
ss = 100
e1 = matrix(rnorm(ss),ncol=1)
e2 = matrix(rnorm(ss),ncol=1)
# Simulate genotypes, gene expression, covariates and a clinical trait
L = matrix(rbinom(ss*3,2,.5),ncol=3)
G = matrix( apply(.4*L, 1, sum) + e1,ncol=1)
T = matrix(.3*G + e2,ncol=1)
T = ifelse( T > median(T), 1, 0 )
C = matrix(matrix(rnorm(ss*2),ncol=1),ncol=2)
n.perm = 5
perm.index = matrix(NA, nrow=ss, ncol=n.perm )
for( j in 1:ncol(perm.index) ) perm.index[, j] = sample( 1:ss )
results = cit.bp(L, G, T, perm.index=perm.index, n.perm=n.perm)
results## perm p_cit p_TassocL p_TassocGgvnL p_GassocLgvnT p_LindTgvnG
## 1 0 0.8320000 0.18811034 0.1451558 0.004160814 0.832
## 2 1 0.7638609 0.76386086 0.7623847 0.276732180 0.604
## 3 2 0.9539845 0.95398448 0.1214126 0.436909584 0.929
## 4 3 0.9365019 0.93650193 0.8357338 0.499334650 0.328
## 5 4 0.9197164 0.09785651 0.4101963 0.919716440 0.463
## 6 5 0.8723721 0.08831553 0.4321350 0.872372064 0.546Causal Inference Test for a Continuous Outcome
# Sample Size
ss = 100
# Errors
e1 = matrix(rnorm(ss),ncol=1)
e2 = matrix(rnorm(ss),ncol=1)
# Simulate genotypes, gene expression, covariates, and clinical trait matrices
L = matrix(rbinom(ss*3,2,.5),ncol=3)
G = matrix( apply(.3*L, 1, sum) + e1,ncol=1)
T = matrix(.3*G + e2,ncol=1)
C = matrix(matrix(rnorm(ss*2),ncol=1),ncol=2)
n.perm = 5
perm.index = matrix(NA, nrow=ss, ncol=n.perm )
for( j in 1:ncol(perm.index) ) perm.index[, j] = sample( 1:ss )
results = cit.cp(L, G, T)
results## p_cit p_TassocL p_TassocGgvnL p_GassocLgvnT p_LindTgvnG
## 0.1156974441 0.1156974441 0.0343861941 0.0006577775 0.0312074251CMatching
Matching Algorithms for Causal Inference with Clustered Data
library(CMatching)## Loading required package: Matching## ##
## ## Matching (Version 4.9-11, Build Date: 2021-10-18)
## ## See http://sekhon.berkeley.edu/matching for additional documentation.
## ## Please cite software as:
## ## Jasjeet S. Sekhon. 2011. ``Multivariate and Propensity Score Matching
## ## Software with Automated Balance Optimization: The Matching package for R.''
## ## Journal of Statistical Software, 42(7): 1-52.
## ##data(schools)
# Kreft and De Leeuw, Introducing Multilevel Modeling, Sage (1988).
# The data set is the subsample of NELS-88 data consisting of 10 handpicked schools
# from the 1003 schools in the full data set.
# Suppose that the effect of homeworks on math score is unconfounded conditional on X
# and unobserved school features (we assume this only for illustrative purpouse).
# Let us consider the following variables:
X<-schools$ses # or X<-as.matrix(schools[,c("ses","white","public")])
Y<-schools$math
Tr<-ifelse(schools$homework>1,1,0)
Group<-schools$schid
# When Group is missing or there is only one Group CMatch returns
# the output of the Match function with a warning.
# Let us assume that the effect of homeworks (Tr) on math score (Y)
# is unconfounded conditional on X and other unobserved school features.
# Several strategies to handle unobserved group characteristics
# are described in Arpino & Cannas, 2016 (see References).
# Multivariate Matching on covariates in X
# default parameters: one-to-one matching on X with replacement with a caliper of 0.25
### Matching within schools
mw<-CMatch(type="within",Y=Y, Tr=Tr, X=X, Group=Group, caliper=0.1)
# compare balance before and after matching
bmw <- CMatchBalance(Tr~X,data=schools,match.out=mw)##
## ***** (V1) X *****
## Before Matching After Matching
## mean treatment........ 0.23211 0.25282
## mean control.......... -0.36947 0.25541
## std mean diff......... 61.315 -0.26998
##
## mean raw eQQ diff..... 0.61617 0.030444
## med raw eQQ diff..... 0.59 0.03
## max raw eQQ diff..... 1.12 0.09
##
## mean eCDF diff........ 0.17801 0.012232
## med eCDF diff........ 0.18324 0.011111
## max eCDF diff........ 0.35227 0.055556
##
## var ratio (Tr/Co)..... 1.2849 1.0154
## T-test p-value........ 3.4478e-07 0.61904
## KS Bootstrap p-value.. < 2.22e-16 0.996
## KS Naive p-value...... 1.979e-07 0.99907
## KS Statistic.......... 0.35227 0.055556 # calculate proportion of matched observations
(mw$orig.treated.nobs-mw$ndrops)/mw$orig.treated.nobs ## [1] 0.6640625 # check number of drops by school
mw$orig.dropped.nobs.by.group## Group
## 7472 7829 7930 24725 25456 25642 62821 68448 68493 72292
## 3 7 8 0 1 2 12 5 3 2 # examine output
mw # complete list of results ## $index.control
## [1] 4 3 8 35 27 35 30 27 30 30 67 65 52 46 72 70 84 74 69
## [20] 81 88 101 100 102 106 122 131 178 158 154 141 161 139 139 156 147 170 147
## [39] 183 170 139 139 164 139 164 170 187 154 187 164 147 187 164 147 164 170 149
## [58] 149 139 147 147 170 147 149 164 187 143 158 183 147 187 139 210 204 205 214
## [77] 202 217 214 204 225 224 220 249 243 254 259 247 254 258
##
## $index.treated
## [1] 5 7 16 28 31 34 38 41 42 43 44 45 48 49 76 78 83 85 86
## [20] 87 89 91 97 103 111 120 128 133 134 135 136 136 137 138 140 142 144 146
## [39] 150 151 153 155 157 159 160 160 163 166 167 168 171 171 172 173 174 176 179
## [58] 180 181 182 184 185 186 186 191 192 193 194 195 197 197 198 201 203 206 207
## [77] 211 215 216 218 226 231 238 241 242 244 248 251 253 256
##
## $index.dropped
## [1] 10 14 18 24 25 29 32 33 39 40 47 51 54 58 59 60 62 64 96
## [20] 118 126 132 145 148 152 162 165 169 175 188 189 190 196 199 208 209 213 219
## [39] 222 223 234 246 252
##
## $weights
## [1] 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
## [20] 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.5 0.5 1.0 1.0 1.0 1.0 1.0 1.0
## [39] 1.0 1.0 1.0 1.0 1.0 1.0 0.5 0.5 1.0 1.0 1.0 1.0 0.5 0.5 1.0 1.0 1.0 1.0 1.0
## [58] 1.0 1.0 1.0 1.0 1.0 0.5 0.5 1.0 1.0 1.0 1.0 1.0 0.5 0.5 1.0 1.0 1.0 1.0 1.0
## [77] 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
##
## $est
## [1] 4.341176
##
## $se
## [1] 1.832278
##
## $mdata
## $mdata$Y
## [1] 43 33 36 35 42 41 37 41 39 36 65 62 43 68 36 62 40 64 54 65 48 33 49 32 43
## [26] 39 47 68 56 68 61 61 61 62 68 71 60 70 67 62 64 67 69 69 63 63 71 65 69 61
## [51] 61 61 60 65 64 70 63 64 56 63 54 55 67 67 70 69 68 52 63 69 69 68 46 65 54
## [76] 58 34 60 50 65 40 37 60 68 50 44 43 62 64 43 42 53 64 32 42 32 35 42 35 35
## [101] 41 54 53 44 33 48 42 40 39 37 45 58 56 46 47 51 44 63 54 57 63 43 60 60 57
## [126] 66 56 66 67 56 60 60 63 60 63 56 58 57 58 63 66 58 63 66 63 56 63 63 60 66
## [151] 66 56 66 63 63 58 61 54 67 66 58 60 40 52 40 41 44 35 41 52 56 43 48 36 43
## [176] 37 34 37 37 47
##
## $mdata$Tr
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [112] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [149] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##
## $mdata$X
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] -0.74 -0.83 -0.50 -1.63 -0.02 -1.67 -1.16 0.03 -1.09 -1.10 0.98 0.24
## [2,] -0.72 -0.80 -0.51 -1.64 -0.02 -1.64 -1.07 -0.02 -1.07 -1.07 0.95 0.29
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] -0.03 0.39 -1.96 -0.48 -1.28 -0.02 -0.60 -0.13 -0.26 -0.75 -0.43 -0.21
## [2,] -0.11 0.35 -2.04 -0.45 -1.31 -0.07 -0.52 -0.14 -0.19 -0.71 -0.49 -0.21
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## [1,] 0.14 -1.05 -0.33 0.48 0.69 0.60 0.77 0.77 1.19 1.17 1.85 1
## [2,] 0.22 -1.11 -0.30 0.40 0.68 0.64 0.75 0.79 1.12 1.12 1.82 1
## [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## [1,] 1.43 0.98 0.18 1.48 1.18 1.13 1.34 1.11 1.40 1.40 1.03 0.56
## [2,] 1.42 1.00 0.26 1.42 1.12 1.12 1.38 1.12 1.38 1.42 1.02 0.64
## [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60]
## [1,] 1.03 1.30 1.01 1.01 1.32 0.97 1.31 1.48 0.86 0.85 1.09 1
## [2,] 1.02 1.38 1.00 1.02 1.38 1.00 1.38 1.42 0.86 0.86 1.12 1
## [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72]
## [1,] 0.98 1.50 0.93 0.93 1.34 1.05 1.8 0.68 0.23 1.01 1.01 1.18
## [2,] 1.00 1.42 1.00 0.86 1.38 1.02 1.8 0.68 0.26 1.00 1.02 1.12
## [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80] [,81] [,82] [,83] [,84]
## [1,] -0.54 0.90 -0.13 -0.88 -0.02 -2.04 -0.99 1.05 -0.29 -0.56 0.82 -0.86
## [2,] -0.55 0.97 -0.14 -0.92 0.01 -1.96 -0.92 0.97 -0.32 -0.49 0.76 -0.84
## [,85] [,86] [,87] [,88] [,89] [,90]
## [1,] -0.41 -0.19 -1.28 -0.67 -0.19 -0.23
## [2,] -0.44 -0.16 -1.19 -0.72 -0.16 -0.29
##
## $mdata$orig.weighted.treated.nobs
## [1] 128
##
## $mdata$orig.weighted.control.nobs
## [1] 132
##
##
## $orig.nobs
## [1] 260
##
## $orig.wnobs
## [1] 260
##
## $orig.weighted.treated.nobs
## [1] 128
##
## $orig.weighted.control.nobs
## [1] 132
##
## $orig.treated.nobs
## [1] 128
##
## $orig.control.nobs
## [1] 132
##
## $wnobs
## [1] 85
##
## $caliper
## [1] 0.1
##
## $intcaliper
## [1] 0.2157358
##
## $exact
## [1] FALSE
##
## $ndrops
## [1] 43
##
## $ndrops.matches
## [1] 43
##
## $estimand
## [1] "ATT"
##
## $orig.treated.nobs.by.group
##
## 7472 7829 7930 24725 25456 25642 62821 68448 68493 72292
## 6 14 12 7 5 4 52 13 6 9
##
## $orig.control.nobs.by.group
##
## 7472 7829 7930 24725 25456 25642 62821 68448 68493 72292
## 17 6 12 15 17 16 15 8 15 11
##
## $orig.ndrops.by.group
##
## 7472 7829 7930 25456 25642 62821 68448 68493 72292
## 3 7 8 1 2 12 5 3 2
##
## $orig.dropped.nobs.by.group
## Group
## 7472 7829 7930 24725 25456 25642 62821 68448 68493 72292
## 3 7 8 0 1 2 12 5 3 2
##
## $version
## [1] "standard"
##
## attr(,"class")
## [1] "CMatch" "Match" summary(mw) # basic statistics##
## Estimate... 4.3412
## SE......... 1.8323
##
## Original number of observations.............. 260
## Original number of treated obs............... 128
## Original number of treated obs by group......
##
## 7472 7829 7930 24725 25456 25642 62821 68448 68493 72292
## 6 14 12 7 5 4 52 13 6 9
##
## Matched number of observations............... 85
## Matched number of observations (unweighted). 90
##
## Caliper (SDs).......................................... 0.1
## Number of obs dropped by 'exact' or 'caliper' ......... 43
## Number of obs dropped by 'exact' or 'caliper' by group
##
## 7472 7829 7930 24725 25456 25642 62821 68448 68493 72292
## 3 7 8 0 1 2 12 5 3 2### Match preferentially within school
# i.e. first match within schools
# then (try to) match remaining units between schools
mpw <- CMatch(type="pwithin",Y=schools$math, Tr=Tr, X=schools$ses,
Group=schools$schid, caliper=0.1)
# examine covariate balance
bmpw<- CMatchBalance(Tr~ses,data=schools,match.out=mpw)##
## ***** (V1) ses *****
## Before Matching After Matching
## mean treatment........ 0.23211 0.23211
## mean control.......... -0.36947 0.23148
## std mean diff......... 61.315 0.063703
##
## mean raw eQQ diff..... 0.61617 0.032297
## med raw eQQ diff..... 0.59 0.02
## max raw eQQ diff..... 1.12 0.19
##
## mean eCDF diff........ 0.17801 0.0095918
## med eCDF diff........ 0.18324 0.0067568
## max eCDF diff........ 0.35227 0.054054
##
## var ratio (Tr/Co)..... 1.2849 1.0055
## T-test p-value........ 3.4478e-07 0.89232
## KS Bootstrap p-value.. < 2.22e-16 0.968
## KS Naive p-value...... 1.979e-07 0.98207
## KS Statistic.......... 0.35227 0.054054 # equivalent to MatchBalance(...) with mpw coerced to class "Match"
# proportion of matched observations
(mpw$orig.treated.nobs-mpw$ndrops) / mpw$orig.treated.nobs ## [1] 1# check drops by school
mpw$orig.dropped.nobs.by.group.after.pref.within## Group
## 7472 7829 7930 24725 25456 25642 62821 68448 68493 72292
## 0 0 0 0 0 0 0 0 0 0# proportion of matched observations after match-within only
(mpw$orig.treated.nobs-sum(mpw$orig.dropped.nobs.by.group.after.within)) / mpw$orig.treated.nobs## [1] 0.6640625# see complete output
mpw## $index.control
## [1] 4 3 8 35 27 35 30 27 30 30 67 65 52 46 72 70 84 74
## [19] 69 81 88 101 100 102 106 122 131 178 158 154 141 161 139 139 156 147
## [37] 170 147 183 170 139 139 164 139 164 170 187 154 187 164 147 187 164 147
## [55] 164 170 149 149 139 147 147 170 147 149 164 187 143 158 183 147 187 139
## [73] 210 204 205 214 202 217 214 204 225 224 220 249 243 254 259 247 254 258
## [91] 106 80 214 114 80 2 20 240 73 73 258 105 125 154 3 4 247 217
## [109] 170 3 4 247 154 202 104 187 73 106 217 143 61 63 53 230 230 143
## [127] 170 230 53 230 73 170 230 143 53 161 170 149 106 177 65 36 37 127
## [145] 61 63 202 73
##
## $index.treated
## [1] 5 7 16 28 31 34 38 41 42 43 44 45 48 49 76 78 83 85
## [19] 86 87 89 91 97 103 111 120 128 133 134 135 136 136 137 138 140 142
## [37] 144 146 150 151 153 155 157 159 160 160 163 166 167 168 171 171 172 173
## [55] 174 176 179 180 181 182 184 185 186 186 191 192 193 194 195 197 197 198
## [73] 201 203 206 207 211 215 216 218 226 231 238 241 242 244 248 251 253 256
## [91] 10 14 18 24 25 29 29 29 32 33 39 40 40 47 51 51 51 54
## [109] 58 59 59 59 60 62 64 64 96 118 126 132 145 145 148 148 152 162
## [127] 162 165 169 169 175 188 189 190 196 199 208 209 213 219 222 223 223 223
## [145] 234 234 246 252
##
## $index.dropped
## numeric(0)
##
## $weights
## [1] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [8] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [15] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [22] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [29] 1.0000000 1.0000000 0.5000000 0.5000000 1.0000000 1.0000000 1.0000000
## [36] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [43] 1.0000000 1.0000000 0.5000000 0.5000000 1.0000000 1.0000000 1.0000000
## [50] 1.0000000 0.5000000 0.5000000 1.0000000 1.0000000 1.0000000 1.0000000
## [57] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 0.5000000
## [64] 0.5000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 0.5000000
## [71] 0.5000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [78] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [85] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [92] 1.0000000 1.0000000 1.0000000 1.0000000 0.3333333 0.3333333 0.3333333
## [99] 1.0000000 1.0000000 1.0000000 0.5000000 0.5000000 1.0000000 0.3333333
## [106] 0.3333333 0.3333333 1.0000000 1.0000000 0.3333333 0.3333333 0.3333333
## [113] 1.0000000 1.0000000 0.5000000 0.5000000 1.0000000 1.0000000 1.0000000
## [120] 1.0000000 0.5000000 0.5000000 0.5000000 0.5000000 1.0000000 0.5000000
## [127] 0.5000000 1.0000000 0.5000000 0.5000000 1.0000000 1.0000000 1.0000000
## [134] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [141] 1.0000000 0.3333333 0.3333333 0.3333333 0.5000000 0.5000000 1.0000000
## [148] 1.0000000
##
## $mdata
## $mdata$Y
## [1] 43 33 36 35 42 41 37 41 39 36 65 62 43 68 36 62 40 64 54 65 48 33 49 32 43
## [26] 39 47 68 56 68 61 61 61 62 68 71 60 70 67 62 64 67 69 69 63 63 71 65 69 61
## [51] 61 61 60 65 64 70 63 64 56 63 54 55 67 67 70 69 68 52 63 69 69 68 46 65 54
## [76] 58 34 60 50 65 40 37 60 68 50 44 43 62 64 43 56 35 41 58 43 38 38 38 41 31
## [101] 44 44 44 63 57 57 57 53 70 51 51 51 67 70 70 70 58 47 42 62 69 69 63 63 51
## [126] 60 60 56 70 70 65 57 63 67 55 69 53 62 46 42 48 62 62 62 71 71 67 64 42 53
## [151] 64 32 42 32 35 42 35 35 41 54 53 44 33 48 42 40 39 37 45 58 56 46 47 51 44
## [176] 63 54 57 63 43 60 60 57 66 56 66 67 56 60 60 63 60 63 56 58 57 58 63 66 58
## [201] 63 66 63 56 63 63 60 66 66 56 66 63 63 58 61 54 67 66 58 60 40 52 40 41 44
## [226] 35 41 52 56 43 48 36 43 37 34 37 37 47 47 41 41 43 41 48 49 44 47 47 47 48
## [251] 47 57 53 42 37 35 56 53 42 37 57 44 60 58 47 47 35 61 44 33 48 40 40 61 56
## [276] 40 48 40 47 56 40 61 48 43 56 63 47 67 54 65 50 47 44 33 44 47
##
## $mdata$Tr
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [186] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [223] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [260] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##
## $mdata$X
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] -0.74 -0.83 -0.50 -1.63 -0.02 -1.67 -1.16 0.03 -1.09 -1.10 0.98 0.24
## [2,] -0.72 -0.80 -0.51 -1.64 -0.02 -1.64 -1.07 -0.02 -1.07 -1.07 0.95 0.29
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] -0.03 0.39 -1.96 -0.48 -1.28 -0.02 -0.60 -0.13 -0.26 -0.75 -0.43 -0.21
## [2,] -0.11 0.35 -2.04 -0.45 -1.31 -0.07 -0.52 -0.14 -0.19 -0.71 -0.49 -0.21
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## [1,] 0.14 -1.05 -0.33 0.48 0.69 0.60 0.77 0.77 1.19 1.17 1.85 1
## [2,] 0.22 -1.11 -0.30 0.40 0.68 0.64 0.75 0.79 1.12 1.12 1.82 1
## [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## [1,] 1.43 0.98 0.18 1.48 1.18 1.13 1.34 1.11 1.40 1.40 1.03 0.56
## [2,] 1.42 1.00 0.26 1.42 1.12 1.12 1.38 1.12 1.38 1.42 1.02 0.64
## [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60]
## [1,] 1.03 1.30 1.01 1.01 1.32 0.97 1.31 1.48 0.86 0.85 1.09 1
## [2,] 1.02 1.38 1.00 1.02 1.38 1.00 1.38 1.42 0.86 0.86 1.12 1
## [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72]
## [1,] 0.98 1.50 0.93 0.93 1.34 1.05 1.8 0.68 0.23 1.01 1.01 1.18
## [2,] 1.00 1.42 1.00 0.86 1.38 1.02 1.8 0.68 0.26 1.00 1.02 1.12
## [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80] [,81] [,82] [,83] [,84]
## [1,] -0.54 0.90 -0.13 -0.88 -0.02 -2.04 -0.99 1.05 -0.29 -0.56 0.82 -0.86
## [2,] -0.55 0.97 -0.14 -0.92 0.01 -1.96 -0.92 0.97 -0.32 -0.49 0.76 -0.84
## [,85] [,86] [,87] [,88] [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96]
## [1,] -0.41 -0.19 -1.28 -0.67 -0.19 -0.23 0.21 -1.58 -0.94 -1.24 -1.54 -0.39
## [2,] -0.44 -0.16 -1.19 -0.72 -0.16 -0.29 0.22 -1.54 -0.92 -1.25 -1.54 -0.39
## [,97] [,98] [,99] [,100] [,101] [,102] [,103] [,104] [,105] [,106] [,107]
## [1,] -0.39 -0.39 0.11 0.13 -0.29 -1.29 -1.29 0.65 -0.76 -0.76 -0.76
## [2,] -0.39 -0.39 0.06 0.06 -0.29 -1.29 -1.29 0.64 -0.80 -0.72 -0.72
## [,108] [,109] [,110] [,111] [,112] [,113] [,114] [,115] [,116] [,117]
## [1,] -1.85 1.48 -0.76 -0.76 -0.76 0.63 0.03 1.05 1.05 0.05
## [2,] -1.96 1.42 -0.80 -0.72 -0.72 0.64 0.01 1.08 1.02 0.06
## [,118] [,119] [,120] [,121] [,122] [,123] [,124] [,125] [,126] [,127]
## [1,] 0.18 -1.92 1.63 -0.12 -0.12 1.23 1.23 1.25 1.61 1.61
## [2,] 0.22 -1.96 1.80 -0.12 -0.12 1.21 1.25 1.25 1.80 1.42
## [,128] [,129] [,130] [,131] [,132] [,133] [,134] [,135] [,136] [,137]
## [1,] 1.28 1.23 1.23 0.09 1.55 1.26 1.68 1.22 0.81 1.55
## [2,] 1.25 1.21 1.25 0.06 1.42 1.25 1.80 1.21 0.79 1.42
## [,138] [,139] [,140] [,141] [,142] [,143] [,144] [,145] [,146] [,147]
## [1,] 0.86 0.18 -0.35 0.29 -0.99 -0.99 -0.99 -0.12 -0.12 0.03
## [2,] 0.86 0.22 -0.35 0.29 -1.01 -1.01 -0.97 -0.12 -0.12 0.01
## [,148]
## [1,] 0.09
## [2,] 0.06
##
## $mdata$orig.weighted.treated.nobs
## [1] 128
##
## $mdata$orig.weighted.control.nobs
## [1] 132
##
##
## $est
## [1] 5.507812
##
## $se
## [1] 2.204627
##
## $orig.nobs
## [1] 260
##
## $orig.wnobs
## [1] 260
##
## $orig.weighted.treated.nobs
## [1] 128
##
## $orig.weighted.control.nobs
## [1] 132
##
## $orig.treated.nobs
## [1] 128
##
## $orig.control.nobs
## [1] 132
##
## $wnobs
## [1] 128
##
## $caliper
## [1] 0.1
##
## $exact
## [1] FALSE
##
## $ndrops
## [1] 0
##
## $ndrops.matches
## [1] 0
##
## $estimand
## [1] "ATT"
##
## $orig.treated.nobs.by.group
##
## 7472 7829 7930 24725 25456 25642 62821 68448 68493 72292
## 6 14 12 7 5 4 52 13 6 9
##
## $orig.dropped.nobs.by.group.after.within
## Group
## 7472 7829 7930 24725 25456 25642 62821 68448 68493 72292
## 3 7 8 0 1 2 12 5 3 2
##
## $orig.ndrops.by.group
## < table of extent 0 >
##
## $orig.dropped.nobs.by.group.after.pref.within
## Group
## 7472 7829 7930 24725 25456 25642 62821 68448 68493 72292
## 0 0 0 0 0 0 0 0 0 0
##
## $version
## [1] "standard"
##
## attr(,"class")
## [1] "CMatch" "Match"# or use summary method for main results
summary(mpw) ##
## Estimate... 5.5078
## SE......... 2.2046
##
## Original number of observations.............. 260
## Original number of treated obs............... 128
## Original number of treated obs by group......
##
## 7472 7829 7930 24725 25456 25642 62821 68448 68493 72292
## 6 14 12 7 5 4 52 13 6 9
##
## Matched number of observations............... 128
## Matched number of observations (unweighted). 148
##
## Caliper (SDs).......................................... 0.1
## Number of obs dropped by 'exact' or 'caliper' ......... 0
## Number of obs dropped by 'exact' or 'caliper' by group
##
## 7472 7829 7930 24725 25456 25642 62821 68448 68493 72292
## 0 0 0 0 0 0 0 0 0 0#### Propensity score matching
# estimate the ps model
mod <- glm(Tr~ses+parented+public+sex+race+urban,
family=binomial(link="logit"),data=schools)
eps <- fitted(mod)
# eg 1: within school propensity score matching
psmw <- CMatch(type="within",Y=schools$math, Tr=Tr, X=eps,
Group=schools$schid, caliper=0.1)
# equivalent to direct call at MatchW(Y=schools$math, Tr=Tr, X=eps,
# Group=schools$schid, caliper=0.1)
# eg 2: preferential within school propensity score matching
psmw <- CMatch(type="pwithin",Y=schools$math, Tr=Tr, X=eps, Group=schools$schid, caliper=0.1)
# Other strategies for controlling unobserved cluster covariates
# via different specifications of propensity score (see Arpino and Mealli):
# eg 3: propensity score matching using ps estimated from a logit model with dummies for hospitals
mod <- glm(Tr ~ ses + parented + public + sex + race + urban
+schid - 1,family=binomial(link="logit"),data=schools)
eps <- fitted(mod)
dpsm <- CMatch(type="within",Y=schools$math, Tr=Tr, X=eps, Group=NULL, caliper=0.1)## Warning in MatchW(Y, Tr, X, Group, estimand, M, exact, caliper, weights, : There
## is only one group or no group has been specified: same output of Match# this is equivalent to run Match with X=eps
# eg4: propensity score matching using ps estimated from multilevel logit model
# (random intercept at the hospital level)
require(lme4)## Loading required package: lme4## Loading required package: Matrixmod<-glmer(Tr ~ ses + parented + public + sex + race + urban + (1 | schid),
family=binomial(link="logit"), data=schools)
eps <- fitted(mod)
mpsm<-CMatch(type="within",Y=schools$math, Tr=Tr, X=eps, Group=NULL, caliper=0.1)## Warning in MatchW(Y, Tr, X, Group, estimand, M, exact, caliper, weights, : There
## is only one group or no group has been specified: same output of Match# this is equivalent to run Match with X=epsbartCause
Causal Inference using Bayesian Additive Regression Trees