Causal analysis of Sesame Street experiment. See Chapters 18 and 21 in Regression and Other Stories.
library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("rstanarm")
library("brms")Set random seed for reproducability
SEED <- 1234sesame <- read.csv(root("Sesame/data","sesame.csv"))
head(sesame) rownames id site sex age viewcat setting viewenc prebody prelet preform
1 1 1 1 1 66 1 2 1 16 23 12
2 2 2 1 2 67 3 2 1 30 26 9
3 3 3 1 1 56 3 2 2 22 14 9
4 4 4 1 1 49 1 2 2 23 11 10
5 5 5 1 1 69 4 2 2 32 47 15
6 6 6 1 2 54 3 2 2 29 26 10
prenumb prerelat preclasf postbody postlet postform postnumb postrelat
1 40 14 20 18 30 14 44 14
2 39 16 22 30 37 17 39 14
3 9 9 8 21 46 15 40 9
4 14 9 13 21 14 13 19 8
5 51 17 22 32 63 18 54 14
6 33 14 14 27 36 14 39 16
postclasf peabody agecat encour X_Isite_2 X_Isite_3 X_Isite_4 X_Isite_5
1 23 62 1 1 0 0 0 0
2 22 8 1 1 0 0 0 0
3 19 32 1 0 0 0 0 0
4 15 27 0 0 0 0 0 0
5 21 71 1 0 0 0 0 0
6 24 32 1 0 0 0 0 0
regular watched encouraged y pretest
1 0 0 1 30 23
2 1 1 1 37 26
3 1 1 0 46 14
4 0 0 0 14 11
5 1 1 0 63 47
6 1 1 0 36 26(sesame_tab <- table(sesame[,c('watched','encouraged')])) encouraged
watched 0 1
0 40 14
1 48 138round(prop.table(sesame_tab, margin=2), digits=2) encouraged
watched 0 1
0 0.45 0.09
1 0.55 0.91Estimate the intent-to-treat (ITT) effect of the instrument (encouragement) on the treatment (regular watching), that is, percentage of children actually induced to watch Sesame Street by the intervention
itt_zt <- stan_glm(watched ~ encouraged, data=sesame, seed=SEED, refresh=0)
print(itt_zt, digits=2)stan_glm
family: gaussian [identity]
formula: watched ~ encouraged
observations: 240
predictors: 2
------
Median MAD_SD
(Intercept) 0.55 0.04
encouraged 0.36 0.05
Auxiliary parameter(s):
Median MAD_SD
sigma 0.38 0.02
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregEstimate the intent-to-treat (ITT) estimate on the outcome (post-treatment letter identification)
itt_zy <- stan_glm(postlet ~ encouraged, data=sesame, refresh=0)
print(itt_zy, digits=1)stan_glm
family: gaussian [identity]
formula: postlet ~ encouraged
observations: 240
predictors: 2
------
Median MAD_SD
(Intercept) 24.9 1.4
encouraged 2.8 1.8
Auxiliary parameter(s):
Median MAD_SD
sigma 13.3 0.6
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregCalculate Wald estimate, ie the ratio of the above two estimates
wald_est <- coef(itt_zy)["encouraged"] / coef(itt_zt)["encouraged"]
round(wald_est, digits=1)encouraged
7.9 The first step is to regress the "treatment" variable---an indicator for regular watching (watched)---on the randomized instrument, encouragement to watch (encouraged).
fit_2a <- stan_glm(watched ~ encouraged, data=sesame, seed=SEED, refresh=0)
print(fit_2a, digits=2)stan_glm
family: gaussian [identity]
formula: watched ~ encouraged
observations: 240
predictors: 2
------
Median MAD_SD
(Intercept) 0.55 0.04
encouraged 0.36 0.05
Auxiliary parameter(s):
Median MAD_SD
sigma 0.38 0.02
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregsummary(fit_2a$fitted, digits=2) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.55 0.55 0.91 0.77 0.91 0.91 sesame$watched_hat <- fit_2a$fittedThen we plug predicted values of watched into the equation predicting the letter recognition outcome.
fit_2b <- stan_glm(postlet ~ watched_hat, data=sesame, seed=SEED, refresh=0)
print(fit_2b, digits = 1)stan_glm
family: gaussian [identity]
formula: postlet ~ watched_hat
observations: 240
predictors: 2
------
Median MAD_SD
(Intercept) 20.6 3.9
watched_hat 7.9 4.9
Auxiliary parameter(s):
Median MAD_SD
sigma 13.3 0.6
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregTwo stage approach with adjusting for covariates in an instrumental variables framework.
The first step is to regress the "treatment" variable on the randomized instrument, encouragement to watch (encouraged) and pre-treatment variables.
fit_3a <- stan_glm(watched ~ encouraged + prelet + as.factor(site) + setting,
data=sesame, seed=SEED, refresh=0)
print(fit_3a, digits=2)stan_glm
family: gaussian [identity]
formula: watched ~ encouraged + prelet + as.factor(site) + setting
observations: 240
predictors: 8
------
Median MAD_SD
(Intercept) 0.66 0.11
encouraged 0.34 0.05
prelet 0.01 0.00
as.factor(site)2 0.03 0.07
as.factor(site)3 -0.11 0.07
as.factor(site)4 -0.34 0.07
as.factor(site)5 -0.29 0.10
setting -0.05 0.05
Auxiliary parameter(s):
Median MAD_SD
sigma 0.35 0.02
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregsummary(fit_3a$fitted, digits=2) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.21 0.64 0.81 0.77 0.97 1.10 Then we plug predicted values of watched into the equation predicting the letter recognition outcome.
watched_hat_3 <- fit_3a$fitted
fit_3b <- stan_glm(postlet ~ watched_hat_3 + prelet + as.factor(site) + setting,
data=sesame, seed=SEED, refresh=0)
print(fit_3b, digits = 1)stan_glm
family: gaussian [identity]
formula: postlet ~ watched_hat_3 + prelet + as.factor(site) + setting
observations: 240
predictors: 8
------
Median MAD_SD
(Intercept) 1.3 4.7
watched_hat_3 14.0 4.2
prelet 0.7 0.1
as.factor(site)2 8.4 1.8
as.factor(site)3 -4.0 1.8
as.factor(site)4 1.0 2.5
as.factor(site)5 2.7 2.9
setting 1.6 1.5
Auxiliary parameter(s):
Median MAD_SD
sigma 9.7 0.4
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregUse the predictor matrix from this second-stage regression.
X_adj <- X <- model.matrix(fit_3b)
X_adj[,"watched_hat_3"] <- sesame$watched
n <- nrow(X)
p <- ncol(X)RMSE1 <- sqrt(sum((sesame$postlet - X %*% coef(fit_3b))^2)/(n-p))
RMSE2 <- sqrt(sum((sesame$postlet - X_adj %*% coef(fit_3b))^2)/(n-p))
se_adj <- se(fit_3b)["watched_hat_3"] * sqrt(RMSE1 / RMSE2)
print(se_adj, digits=2)watched_hat_3
4.2 f1 <- bf(watched ~ encour)
f2 <- bf(postlet ~ watched)IV_brm_a <- brm(f1 + f2, data=sesame, seed=SEED)print(IV_brm_a, digits=1) Family: MV(gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: watched ~ encour
postlet ~ watched
Data: sesame (Number of observations: 240)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
watched_Intercept 0.5 0.0 0.5 0.6 1.0 4360 2975
postlet_Intercept 20.6 3.8 13.5 28.4 1.0 1976 2405
watched_encour 0.4 0.1 0.3 0.5 1.0 4335 2665
postlet_watched 7.9 4.8 -1.9 16.9 1.0 1958 2150
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_watched 0.4 0.0 0.4 0.4 1.0 4121 2791
sigma_postlet 12.7 0.7 11.4 14.2 1.0 2971 2170
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
rescor(watched,postlet) 0.2 0.2 -0.1 0.5 1.0 2108
Tail_ESS
rescor(watched,postlet) 2117
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).f1 <- bf(watched ~ encour + prelet + setting + factor(site))
f2 <- bf(postlet ~ watched + prelet + setting + factor(site))IV_brm_b <- brm(f1 + f2, data=sesame, seed=SEED)print(IV_brm_b, digits=1) Family: MV(gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: watched ~ encour + prelet + setting + factor(site)
postlet ~ watched + prelet + setting + factor(site)
Data: sesame (Number of observations: 240)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
watched_Intercept 0.7 0.1 0.5 0.9 1.0 4353 2782
postlet_Intercept 1.4 4.6 -8.1 10.2 1.0 1845 1696
watched_encour 0.3 0.1 0.2 0.4 1.0 4902 3139
watched_prelet 0.0 0.0 -0.0 0.0 1.0 6567 2378
watched_setting -0.1 0.1 -0.2 0.0 1.0 4974 2889
watched_factorsite2 0.0 0.1 -0.1 0.2 1.0 3803 2777
watched_factorsite3 -0.1 0.1 -0.2 0.0 1.0 3255 2857
watched_factorsite4 -0.3 0.1 -0.5 -0.2 1.0 3805 3189
watched_factorsite5 -0.3 0.1 -0.5 -0.1 1.0 3981 3031
postlet_watched 13.9 3.9 6.5 21.8 1.0 1772 1648
postlet_prelet 0.7 0.1 0.5 0.9 1.0 4456 2728
postlet_setting 1.5 1.4 -1.3 4.4 1.0 3386 2894
postlet_factorsite2 8.4 1.8 4.9 11.9 1.0 4151 3085
postlet_factorsite3 -4.0 1.8 -7.4 -0.5 1.0 3415 2988
postlet_factorsite4 0.9 2.4 -3.8 5.8 1.0 2116 2481
postlet_factorsite5 2.7 2.9 -2.8 8.5 1.0 2800 2725
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_watched 0.4 0.0 0.3 0.4 1.0 5925 2720
sigma_postlet 9.4 0.5 8.5 10.7 1.0 2718 1818
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
rescor(watched,postlet) -0.2 0.2 -0.5 0.1 1.0 1794
Tail_ESS
rescor(watched,postlet) 1762
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).