dataAnalysis.Rnw

\documentclass{article}
\usepackage{amsmath,amssymb}
\usepackage{hyperref}
\usepackage{bm}

\newcommand{\yti}{Y_{Ti}}
\newcommand{\yci}{Y_{Ci}}
\newcommand{\uti}{U_{Ti}}
\newcommand{\uci}{U_{Ci}}
\newcommand{\etat}{\eta_T}
\newcommand{\etati}{\eta_{Ti}}

\newcommand{\mti}{\bar{m}_{Ti}}
\newcommand{\byt}{\bm{Y_T}}
\newcommand{\byc}{\bm{Y_C}}
\newcommand{\bmt}{\bm{\bar{m}_T}}
\newcommand{\bmi}{\bm{m}_i}
\newcommand{\bsi}{\bm{s}_i}

\newcommand{\EE}{\mathbb{E}}

\title{Data Analysis for 'The Role of Mastery Learning'}

<<prelim,include=FALSE>>=
library(knitr)
library(xtable)
opts_chunk$set(
echo=TRUE,results='asis',warning = FALSE,message = FALSE,cache=TRUE,error = FALSE
)

if(!exists('runMods')) runMods <- FALSE

@


\begin{document}
This document includes code to produce all of the results and run all
of the models reported in ``The Role of Mastery Learning in
Intelligent Tutoring Systems: Principal Stratification on a Latent Variable.''

The auxilliary files sourced here are available at our github repository,
\url{https://github.com/adamSales/ctaiAdvance}.

First, load in and transform the (pre-imputed) data:
<<data>>=
load('data/HSdata.RData')
load('data/advanceData.RData')
@
We'll use the \texttt{R} package \texttt{rstan} to run the models:
<<packages>>=
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
@

\section{Data Description (Section \ref{sec:data})}

This code produces the missigness information from Table 1, summarizing the student level data:
<<missignessTable>>=
miss <- NULL
for(i in c('race','sex','spec','xirt')) miss <- rbind(miss,
 c(sum(is.na(covs[[i]])),mean(is.na(covs[[i]])),error[i,'error']))
miss <- as.data.frame(miss)
miss$`Error Type` <- c('PFC','PFC','PFC','SRMSE')
rownames(miss) <- c('Race/Ethnicity','Sex','Special Education','Pretest')
names(miss)[1:3] <- c('# Missing','% Missing','Imputation Error')
miss[,2] <- as.integer(round(miss[,2]*100))
miss[,1] <- as.integer(miss[,1])
miss['Pretest','Imputation Error'] <- sqrt(miss['Pretest','Imputation Error'])/sd(covs$xirt,na.rm=TRUE)


print(xtable::xtable(miss))
@
This code produces the covariate balance information:
<<covBal>>=
covBal <- NULL
for(i in c('race','sex','spec')){
    covBal <- rbind(covBal,c(i,NA,NA,NA,NA))

    for(ll in levels(dat[[i]])){
        covBal <- rbind(covBal,c(NA,ll,round(c(mean(dat[[i]]==ll),mean(dat[[i]][dat$treatment==1]==ll),mean(dat[[i]][dat$treatment==0]==ll)),2)))

    }
}
colnames(covBal) <- c('Covariate','Category','Overall Percent','Percent of Treated','Percent of Control')
print(xtable(covBal),floating=FALSE,include.rownames=FALSE)
@
The overall p-value for balance is:
<<overalBal,results='markup'>>=
library(RItools) ## using development version
balMod <- balanceTest(treatment~poly(xirt,2)+spec+race+sex+strata(pair)+cluster(schoolid2),data=dat,report='chisquare.test')
print(balMod$overall['pair',])
@

\section{PS Model with $\bar{m}_T$}

Here we estimate the model in Section
\ref{sec:principalStratification} stratifying on $\bar{m}_T$.

First, we create the datasets:
<<obsData,dependson='data'>>=
source('R/prelimMbar.r')
@

The model is encoded in the file \texttt{psmodObs.stan}.
It may be summarized as follows.
The model for $\bar{m}_T$ is:
\begin{equation}\label{eq:mbarUsage}
\bar{m}_{Ti}=\alpha^U_s+\bm{x}_i^T\bm{\beta^U}+\epsilon^{Ui}_i+\epsilon^{Ut}_{t[i]}
\end{equation}
where $\alpha^U_s$ is a separate intercept for each state, and $\bm{x}_i$ is a vector of covariates: dummy variables for
racial/ethnic category, a dummy variable for sex, dummy variables for
special education category, and linear and quadratic terms for
pretest.
The normally-distributed errors $\epsilon^{Ui}$ and $\epsilon^{Ut}_{t[i]}$ vary at the
individual and teacher levels.

The model for $Y$ is
\begin{equation}
Y_i=\alpha^Y_p+\bm{x}_i^T\bm{\beta^Y}+a_1\bar{m}_{Ti}+Z_i(b_0+b_1*\bar{m}_{Ti})+\epsilon^{Yi}_i+\epsilon^{Yt}_{t[i]}+\epsilon^{Ys}_{s[i]}
\end{equation}
where $\alpha^Y_p$ is a separate intercept for each randomization
block $p$, $Z_i$ is a dummy variable for treatment status,
$\epsilon^{Ys}_{s[i]}$ is a normally distributed error at the school
level, and the rest of the variables are analogous to those in
(\ref{eq:mbarUsage}).
We run the model with the \texttt{stan} command from \texttt{rstan}:
<<mbarModel,dependson='obsData',eval=runMods>>=
mbarMod <- stan('R/psmodObs.stan',data=sdatObs,seed=613)
@
<<loadMbarModel,eval=!exists('mbarModel')>>=

@

Figure \ref{fig:mbarModel} can be replicated with the following code:

<<mbarModelFig>>=
library(tikzDevice) ## allows latex code in figure
options( tikzLatexPackages = c(
getOption( "tikzLatexPackages" ),
"\\usepackage{amsmath,amsfonts}"
))

draw <- 1000

samps <- extract(mbarMod)
plotDatObs <- with(sdatObs,data.frame(Y=c(YtO,YtM,Yc),mbar=c(MbarTO,samps$MbarTM[draw,],samps$MbarC[draw,]),Z=c(rep(1,nstudTO),rep(1,nstudTM),rep(0,nstudC))))
plotDatObs$treat <- ifelse(plotDatObs$Z==1,'Treatment','Control')
plotDatObs$slope <- ifelse(plotDatObs$treat=='Control',samps$a1[draw],samps$a1[draw]+samps$b1[draw])
plotDatObs$int <- ifelse(plotDatObs$treat=='Control',samps$a0[draw],samps$a0[draw]+samps$b0[draw])

plotDatObs <- within(plotDatObs, int <- int-( mean(int+slope*mbar)-mean(plotDatObs$Y)))
plotDatObs <- plotDatObs[order(plotDatObs$treat),]
plotDatObs$treat2 <- plotDatObs$treat

tikz(file = "figure/mbarModel.tex",
  standAlone = T,
  width  = 6, height  = 3)

ggplot(plotDatObs,aes(mbar,Y,fill=treat,group=treat,alpha=treat,color=treat))+geom_point(size=2)+
    geom_abline(aes(intercept=int,slope=slope,linetype=treat2),color='black',size=2,alpha=1)+scale_alpha_discrete(range=c(0.4,.8))+
    scale_colour_manual(values=c('red','blue'))+
    labs(group=NULL,fill=NULL,alpha=NULL)+xlab('$\\bar{m}_T$')+
    ylab('Posttest Score')+theme(legend.position='top',text=element_text(size=20))+
    guides(color = guide_legend(title=NULL,override.aes=list(alpha=1),keywidth=3),linetype=guide_legend(title=NULL,keywidth=1))#override.aes=list(size=2)))
dev.off()
setwd('figure'); tools::texi2dvi('mbarModel.tex', pdf = T, clean = T); setwd('..')
@
\includegraphics{figure/mbarModel.pdf}

To save memory, save and delete the $\bar{m}_T$ model:
<<mbarSave>>=
save(mbarMod,file='output/mbarMod.RData')
@
<<mbarDelete,results='hide'>>=
rm(mbarMod); gc()
@


\section{The Main PS Model}
This section reproduces our paper's main model, described in Section \ref{sec:themodel}.

The data for the main model (similar to the $\bar{m}$ model but
including student-section level mastery data) relies on a secondary
file (available at github):
<<prelimMain,dependson='data',cache=FALSE,eval=(nrow(dat)!=5308)>>=
source('R/prelimStan.r')
@

\newpage
Since this is the main model, we will include full stan code in this
online supplement:
\begin{verbatim}
<<mainModel,echo=FALSE,results='asis'>>=
cat(readLines('R/psmod.stan'),sep='\n')
@
\end{verbatim}

This code runs the model:
<<runmainModel,dependson='prelimMain',eval=runMods,cache=FALSE>>=
main <- stan('R/psmod.stan',data =sdat,warmup=1500,chains=10,iter=5000)
save(main,file='output/mainModel.RData')
@



<<saveSmallMain,eval=runMods,cache=FALSE>>=
draws <- extract(main)

### for "multImp" and "trtEff"
set.seed(613)
U <- draws$studEff
Usamp <- U[sample(1:nrow(U),1000),]

### for sampleSizeEta & etaDiff
draw <- 1000
U <- U[,sort(unique(sdat$studentM))]
eta <- U[draw,]
etasd <- apply(U,2,sd)

### for "usageModel"
sdEta <- sqrt(mean(apply(draws$studEff,1,var)))
Eeta <- colMeans(draws$studEff)

draws$studEff <- Usamp

summMain <- summary(main)[[1]]

save(draws,draw,eta,etasd,sdEta,Eeta,summMain,Usamp,file='output/smallMain.RData')
@

<<loadMain,echo=FALSE,eval=!exists('draws'),cache=FALSE>>=
load('output/smallMain.RData')
@

\section{Multiple Imputation Model Fit}

To give some intuition on how the model fitting worked, and to what
extent treatment effect moderation was discernable in this dataset
anyway, we re-fit the model using (something akin to) multiple
imputation.
First, extract 1000 random draws of $\eta_T$ (denoted as
\texttt{studEff} in the model code) from the fitted model.
Then, for each draw, fit a standard HLM interacting treatment with the
$\eta_T$ draw.
<<multImp,dependson='prelimMain'>>=
library(lme4)
set.seed(613)
## U <- extract(main,'studEff')[[1]]
## Usamp <- U[sample(1:nrow(U),1000),]
multImp <- apply(Usamp,1,
                 function(u) summary(
                                 lmer(Y~treatment*u+poly(xirt,2)+race+sex+spec+state+(1|schoolid2)+(1|teachid2),
                                      data=dat))$coef['treatment:u',])
@

Does this give similar answers to the main model?
<<multImpVsMain>>=
mean(multImp['Estimate',])

## use "Rubin's Rule" to estimate SE
sqrt(var(multImp['Estimate',])+mean(multImp['Std. Error',]^2))

## now for the main model:
## summary(main,par='b1',probs=c())[[1]]
summMain['b1',c(1:3,9:10)]
@

\section{Figures Comparing PS with $\bar{m}_T$ to PS with $\eta$}

Figure \ref{fig:sampleSizeMbar}:
<<sampleSizeMbar>>=
## smart jittering:
datObs$mbarJ <- datObs$mbar
datObs$nsecJ <- datObs$nsec
tab <- with(datObs,table(mbar,nsec))
mult <- which(tab>1,arr.ind=TRUE)
ms <- sort(unique(datObs$mbar))
ns <- sort(unique(datObs$nsec))
for(i in 1:nrow(mult)){
    w <- which(datObs$mbar==ms[mult[i,'mbar']] & datObs$nsec==ns[mult[i,'nsec']])
    s <- length(w)
    if(s>1){
        width=min(s*0.002,0.01)
        height=min(s*0.2,2)
        datObs$nsecJ[w] <- datObs$nsecJ[w]+runif(s,-width,width)
        datObs$mbarJ[w] <- datObs$mbarJ[w]+runif(s,-width,width)
    }
}



tikz(file='figure/mbarSampleSize.tex',
     standAlone=T,
     width=3,height=3)
ggplot(datObs,aes(mbarJ,nsecJ))+geom_point()+xlab('$\\bar{m}$')+ylab('$n_{sec}$')+theme(text=element_text(size=20))
dev.off()
setwd('figure'); tools::texi2dvi('mbarSampleSize.tex', pdf = T, clean = T); setwd('..')
@
\includegraphics{figure/mbarSampleSize.pdf}

Figure \ref{fig:etaSampleSize}:
<<sampleSizeEta>>=
## draw <- 1000
sdatLat <- sdat
nsec <- as.vector(table(sdatLat$studentM))
## etaDraws <- extract(main,'studEff')[[1]][,sort(unique(sdatLat$studentM))]
## eta <- etaDraws[draw,]
## etasd <- apply(etaDraws,2,sd)

plotDat <- data.frame(nsec=nsec,eta=eta,etasd=etasd)

tikz(file='figure/etaSampleSize.tex',
     standAlone=T,
     width=3,height=3)
ggplot(plotDat,aes(eta,nsec,size=1/etasd))+geom_point()+ylab(NULL)+#ylab('$n_{sec}$')+
    labs(size='$1/\\text{SE}(\\eta_T)$')+scale_size(range=c(.5,2))+guides(size=FALSE)+xlab('$\\eta_T$')+
    theme(text=element_text(size=20))+
    theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())#+ggtitle('One Posterior Draw')#+xlab('$\\mathbb{E}\\eta$')
dev.off()
setwd('figure'); tools::texi2dvi('etaSampleSize.tex', pdf = T, clean = T); setwd('..')

@
\includegraphics{figure/etaSampleSize.pdf}

Figure \ref{fig:mbarDiff}:
<<mbarDiff>>=
secDiff <- -draws$secEff[draw,]
sss <- secDiff[sdatLat$sec]
mDiff <- aggregate(sss,list(stud=sdatLat$studentM),mean)
mbar <- aggregate(sdatLat$grad,list(sdatLat$studentM),mean)

mbarDiffDat <- data.frame(mbar=mbar$x,mDiff=mDiff$x)


tikz(file='figure/mbarDiff.tex',
     standAlone=T,
     width=3,height=3)
ggplot(mbarDiffDat,aes(mbar,mDiff))+geom_point()+xlab('$\\bar{m}$')+ylab('Avg. Sec. Difficulty')+
    theme(text=element_text(size=20))
dev.off()
setwd('figure'); tools::texi2dvi('mbarDiff.tex', pdf = T, clean = T); setwd('..')
@
\includegraphics{figure/mbarDiff.pdf}

Figure \ref{fig:etaDiff}:
<<etaDiff>>=
plotDat$mDiff <- mDiff$x
tikz(file='figure/etaDiff.tex',
     standAlone=T,
     width=3,height=3)
ggplot(plotDat,aes(eta,mDiff,size=1/etasd))+geom_point()+#ylab('$Avg. Section Difficulty$')+
    labs(size='$1/\\text{SE}(\\eta_T)$')+scale_size(range=c(.5,2))+guides(size=FALSE)+xlab('$\\eta_T$')+
    theme(text=element_text(size=20))+
    theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())#+ggtitle('One Posterior Draw')#+xlab('$\\mathbb{E}\\eta$')
dev.off()
setwd('figure'); tools::texi2dvi('etaDiff.tex', pdf = T, clean = T); setwd('..')
@
\includegraphics{figure/etaDiff.pdf}

\section{Main Model Results}
Get all the MCMC draws from the main model:
% <<mainDraws>>=
% draws <- extract(main)
% @
\subsection{Predicting Mastery}
<<usageModel>>=
library(coefplot)

## coefs <- summary(main,'betaU')[[1]]
coefs <- summMain[grep('betaU',rownames(summMain)),]
rownames(coefs) <- colnames(sdat$X)

## sqrt of average (over the draws) of the variance of eta
## sdEta <- sqrt(mean(apply(draws$studEff,1,var)))

coefs <- coefs/apply(sdat$X,2,sd)/sdEta

cpdf <- data.frame(Value=coefs[3:7,1],Coefficient=c('Black/\n Multiracial','Hispanic/\n Native American','Male','Special Ed.','Gifted'),HighInner=coefs[3:7,'75%'],LowInner=coefs[3:7,'25%'],HighOuter=coefs[3:7,'97.5%'],LowOuter=coefs[3:7,'2.5%'],Model='model')

coefplot.data.frame(cpdf,title=NULL,xlab=expression(hat(beta)[std]),ylab=NULL,lwdOuter=0.5,lwdInner=1.5)+theme(text=element_text(size=12))
ggsave('figure/usageCoef.pdf',width=6,height=3)

pdf('figure/pretestEta.pdf')
omar <- par()$mar
par(mar=omar+c(0,1,-1,0))
## Eeta <- colMeans(draws$studEff)
plot(dat$xirt,Eeta/sdEta,col=ifelse(dat$treatment==1,'blue','red'),
     xlab='Pretest (std)',ylab=expression(paste('E[',eta[T],'|x]')),cex.lab=2)
X <- scale(model.matrix(~poly(xirt,2),data=dat)[,-1])
xpred <- (X[,1]*mean(draws$betaU[,1])+X[,2]*mean(draws$betaU[,2]))/sdEta
set.seed(613)
samp <- sample(1:4000,100)
for(ss in samp){
    xpredS <- (X[,1]*draws$betaU[ss,1]+X[,2]*draws$betaU[ss,2])/sdEta
    lines(sort(dat$xirt),xpredS[order(dat$xirt)],col=adjustcolor('pink',0.5))
}
lines(sort(dat$xirt),xpred[order(dat$xirt)],lwd=2)
legend('bottomright',legend=c('Treated','Control (Imputed)','Model (Avg.)','Model (draws)'),col=c('blue','red','black','pink'),pch=c('o','o','.','.'),lwd=c(0.01,0.01,2,2))
dev.off()

@

\subsection{CTAI Treatment Effects}

<<trtEff,results='hide'>>=
pdf('figure/treatmentEffects.pdf')
set.seed(613)
## samp <- sample(1:nrow(draws$studEff),1000)
## Usamp <- draws$studEff[samp,]
studEff95 <- quantile(Usamp,c(0.025,0.975))
Usamp[Usamp<studEff95[1] | Usamp>studEff95[2]] <- NA
trtEff <- sweep(sweep(Usamp,1,draws$b1[samp],'*'),1,draws$b0[samp],'+')

plot(1,1,type='n',ylim=range(trtEff,na.rm=TRUE),
      xlim=studEff95,xlab=expression(eta[T]),ylab='Treatment Effect',cex.lab=2)
sapply(samp,function(rr) invisible(abline(draws$b0[rr],draws$b1[rr],col=adjustcolor('red',.3))))
abline(mean(draws$b0),mean(draws$b1),lwd=3)
legend('topright',legend=c('Mean Est. Effect','MCMC Draws'),lty=1,col=c('black','red'))
dev.off()

pdf('figure/potentialOutcomes.pdf')
a0 <- rnorm(length(draws$a1),mean(sdat$Y[sdat$Z==0]),sd(sdat$Y[sdat$Z==0])/sqrt(sum(sdat$Z==0)))
a1 <- draws$a1
b0 <- draws$b0
b1 <- draws$b1


xx <- seq(studEff95[1],studEff95[2],length=100)
Yc <- outer(a1,xx)
Yc <- sweep(Yc,1,a0,'+')

YcUp <- apply(Yc,2,function(x) quantile(x,0.975))
YcDown <- apply(Yc,2,function(x) quantile(x,0.025))


Yt <- outer(a1+b1,xx)
Yt <- sweep(Yt,1,a0+b0,'+')
YtUp <- apply(Yt,2,function(x) quantile(x,0.975))
YtDown <- apply(Yt,2,function(x) quantile(x,0.025))


curve(mean(a0)+mean(a1)*x,from=min(xx), to=max(xx),lwd=2,col='red',xlab=expression(eta[T]),ylab=expression(paste('E[',Y[Z],'|',eta[T],']',sep='')),ylim=range(c(YtDown,YcDown,YtUp,YcUp)))

curve(mean(a0)+mean(b0)+(mean(b1)+mean(a1))*x,add=TRUE,lwd=2,col='blue')
polygon(c(xx,rev(xx)),c(YcUp,rev(YcDown)),col=adjustcolor('red',0.1))
polygon(c(xx,rev(xx)),c(YtUp,rev(YtDown)),col=adjustcolor('blue',0.1))

 legend('topleft',legend=c(expression(Y[C]),expression(Y[T])),col=c('red','blue'),lwd=2)
 dev.off()
@
\includegraphics{figure/treatmentEffects.pdf}
\includegraphics{figure/potentialOutcomes.pdf}

\subsection{Fake Models}
<<>>=
load('../advance/fittedModels/noEffect.RData')

plotFake <- function(fakeMod,expr){
 set.seed(613)
 draws <- extract(fakeMod)
 samp <- sample(1:nrow(draws$studEff),1000)
 Usamp <- draws$studEff[samp,]
 studEff95 <- quantile(Usamp,c(0.025,0.975))
 Usamp[Usamp<studEff95[1] | Usamp>studEff95[2]] <- NA
 trtEff <- sweep(sweep(Usamp,1,draws$b1[samp],'*'),1,draws$b0[samp],'+')

 plot(1,1,type='n',ylim=range(trtEff,na.rm=TRUE),
      xlim=studEff95,xlab=expression(eta[T]),ylab='Treatment Effect',cex.lab=2)
 sapply(samp,function(rr) invisible(abline(draws$b0[rr],draws$b1[rr],col=adjustcolor('red',.3))))
 abline(mean(draws$b0),mean(draws$b1),lwd=3)
 curve(expr,add=TRUE,lwd=3,lty=2,from=studEff95[1]-1,to=studEff95[2]+1)
 legend('topright',legend=c('Mean Est. Effect','MCMC Draws','True Effect'),lty=c(1,1,2),col=c('black','red','black'))
}
@


\end{document}