MLM in Bayes
# -------------------------------------------------
# get the L1 vars
# -------------------------------------------------
### receptive vocab (DV): PPVT
PPVT <- c("PPVTEPF", "PPVTEPS")
### prereading (DV): letter-word identification from WJIII WJ1_SSEPF: WJIII ENG PK F LETTER-WORD ID STAND SCORE
preread <- c("WJ1_SSEPF", "WJ1_SSEPS")
### early math (DV): applied problems subtest of WJIII
math <- c("WJ10_SSEPF", "WJ10_SSEPS")
### controls
l1controls <- c("ICPSR_CLASS_ID", "CHGENP", "POORP")
l1vars <- c(l1controls, PPVT, preread, math)
# -------------------------------------------------
# get the L2 vars
# -------------------------------------------------
### teacher education
TeachEd <- c("T_EDUCP", "TMAJORP")
### classroom quality: ECERTOTP: ECERS-R: PREK: TOTAL (ITEMS 1-10,12-36) PK
ClassQual <- c("ECERTOTP")
### controls
l2controls <- c("ICPSR_CLASS_ID", "STATE", "RATIOPSP", "CLS_SIZEPS", "HOURSWKP", "T_RACEP", "WHITEP", "PERPOORP", "LOCATIONP", "FULLDAYP")
l2vars <- c(l2controls, TeachEd, ClassQual)
# -------------------------------------------------
# prep the data file for jags
# -------------------------------------------------
str(df)
jags.dat <- df[order(df$ICPSR_CLASS_ID), ] # order by L2 units (classroom ID).
jags.dat <- as.list(jags.dat[, l1vars[2:9]]) # selecting only the L1 vars to be used in the regressions
jags.dat$CLASS <- as.numeric(as.ordered(df$ICPSR_CLASS_ID)) # save L2 unit counter. note that this must be numeric.
str(jags.dat) #jags chokes on factors
jags.dat$CHGENP <- as.numeric(jags.dat$CHGENP) - 1
jags.dat$POORP <- as.numeric(jags.dat$POORP) - 1
# some weirdness in the datafile...
jags.dat$PPVTEPF <- as.numeric(jags.dat$PPVTEPF)
jags.dat$WJ1_SSEPF <- as.numeric(jags.dat$WJ1_SSEPF)
jags.dat$WJ10_SSEPF <- as.numeric(jags.dat$WJ10_SSEPF)
summary(df[,l2vars])
# L2 vars should be the length of the number of L2 units! I have 704 classrooms.
jags.dat$RATIOPSP <- as.vector(by(as.numeric(df$RATIOPSP), df$ICPSR_CLASS_ID, mean))
jags.dat$CLS_SIZEPS <- as.vector(by(as.numeric(df$CLS_SIZEPS), df$ICPSR_CLASS_ID, mean))
jags.dat$HOURSWKP <- as.vector(by(as.numeric(df$HOURSWKP), df$ICPSR_CLASS_ID, mean))
jags.dat$WHITEP <- as.vector(by(as.numeric(df$WHITEP), df$ICPSR_CLASS_ID, mean))
jags.dat$PERPOORP <- as.vector(by(as.numeric(df$PERPOORP), df$ICPSR_CLASS_ID, mean))
jags.dat$LOCATIONP <- as.vector(by(as.numeric(df$LOCATIONP), df$ICPSR_CLASS_ID, mean)) - 1
jags.dat$FULLDAYP <- as.vector(by(as.numeric(df$FULLDAYP), df$ICPSR_CLASS_ID, mean)) - 1
#jags.dat$T_EDUCP <- as.vector(by(as.numeric(df$T_EDUCP), df$ICPSR_CLASS_ID, mean)) # this needs to be dummy coded :(
#jags.dat$TMAJORP <- as.vector(by(as.numeric(df$TMAJORP), df$ICPSR_CLASS_ID, mean)) # this needs to be dummy coded :(
jags.dat$ECERTOTP <- as.vector(by(as.numeric(df$ECERTOTP), df$ICPSR_CLASS_ID, mean))
# need to add L1 n, and L2 n
jags.dat$N.child <- nrow(df)
jags.dat$N.class <- length(unique(df$ICPSR_CLASS_ID))
# check out the structure of the data file
str(jags.dat)
# -------------------------------------------------
# write the jags model
# -------------------------------------------------
jags.model <- function(){
for (i in 1:N.child){
PPVTEPS[i] ~ dnorm (mu[i], tau)
mu[i] <- b.class[CLASS[i]] + b.PPVTEPF * PPVTEPF[i] + b.CHGENP * CHGENP[i] + b.POORP * POORP[i]
}
b.PPVTEPF ~ dnorm(0, 0.01)
b.CHGENP ~ dnorm(0, 0.01)
b.POORP ~ dnorm (0, 0.01)
for (j in 1:N.class){
b.class[j] ~ dnorm(b.class.hat[j], tau.class)
b.class.hat[j] <- b.class0 + b.RATIOPSP * RATIOPSP[j] +
b.CLS_SIZEPS * CLS_SIZEPS[j] + b.HOURSWKP * HOURSWKP[j] + b.WHITEP * WHITEP[j] + b.PERPOORP * PERPOORP[j] + b.LOCATIONP * LOCATIONP[j] + b.FULLDAYP * FULLDAYP[j] + b.ECERTOTP * ECERTOTP[j]
}
b.class.hat.mu <- mean(b.class.hat[])
b.class0 ~ dnorm(0,.01)
b.RATIOPSP ~ dnorm(0,.01)
b.CLS_SIZEPS ~ dnorm(0,.01)
b.HOURSWKP ~ dnorm(0,.01)
b.WHITEP ~ dnorm(0,.01)
b.PERPOORP ~ dnorm(0,.01)
b.LOCATIONP ~ dnorm(0,.01)
b.FULLDAYP ~ dnorm(0,.01)
b.ECERTOTP ~ dnorm(0,.01)
tau ~ dgamma(1,1)
sigma.class <- 1 / tau.class
tau.class ~ dgamma(1,1)
}
# -------------------------------------------------
# set the parameters to track and run the model
# -------------------------------------------------
params <- c(".PPVTEPF", "b.CHGENP", "b.POORP", "b.RATIOPSP", "b.CLS_SIZEPS", "b.HOURSWKP", "b.WHITEP", "b.PERPOORP", "b.LOCATIONP", "b.FULLDAYP", "b.ECERTOTP", "b.class")
library(R2jags); set.seed(123)
fit <- jags(data = jags.dat, inits = NULL, parameters.to.save = params,
model.file = jags.model, n.chains = 2, n.iter = 2000, n.burnin = 1000)
# Not working yet. :)