I am interested in building multilevel structural equation model using international large-scale student survey data. The model is hierarchical, with response items from students, teachers and school admins, used to predict student achievement. Furthermore, the schools are nested within countries and different years. The response items on student, teacher and school level is mainly ordinal (Likert-scale), with a few categorical background variables and continuous achievement variables at student level. From this data I want to build a multivariate
multilevel model with
- latent predictors, and covariates, with ordinal indicators at each level,
- within-level and cross-level latent moderation/interaction,
- dependent latent variables with continuous indicators, and
- ideally, a way to account for missing data.
I started out using OpenMx due to its flexibility with multilevel models and integration in R (also Linux support, and being freely available open-source software). However, due to limitations in the software, such as lack of support ordinal data and WLS in multilevel models, I’ve been advised that this may only be achieved using Bayesian methods.
However, I am not very familiar with Bayesian statistics or the software used. I looked into OpenBUGS and Stan, which both can be used with R in Linux. However, the syntax is unfamiliar to me, and I can find little information on mulitlevel SEM in Stan (this is the only example I’ve come across).
I could really use some help/advise on how to (if at all possible) fit multilevel SEMs in Stan, with the aformentioned requirements. Are there any articles with descriptive methods, documentation, tutorials, or examples on how to build a model fitting the description above?
Below I’ve included an example model in OpenMx, I was hoping somebody could help me translate:
MPlus: Three-level growth model with a continuous outcome and one covariate on each of the three levels
# https://www.statmodel.com/usersguide/chapter9.shtml
Sys.setenv(OMP_NUM_THREADS=parallel::detectCores())
library(OpenMx)
library(dplyr)
library(magrittr)
#mxOption(NULL, "Number of Threads", 8L)
options(width=120)
ex923 <- suppressWarnings(try(read.table("ex9.23.dat")))
#if (is(ex923, "try-error")) ex923 <- read.table("ex9.23.dat")
colnames(ex923) <- c(paste0('y',1:4), 'x', 'w', 'z', 'level2', 'level3')
ex923$level2 <- as.integer(ex923$level2)
ex923$level3 <- as.integer(ex923$level3)
level3Model <- mxModel(
'level3Model', type='RAM',
latentVars=c(paste0('y',1:4), 'ib3', 'sb3', 'z'),
productVars=c("ib3×sb3"),
mxPath(from=c("ib3","sb3"), to="ib3×sb3", arrows=1, free=FALSE, values=1),
mxPath(from="ib3×sb3", to=paste0('y',1:4), arrows=1, free=TRUE, values=.5), #, labels="b1"),
mxPath('ib3', paste0('y',1:4), free=FALSE, values=1),
mxPath('sb3', paste0('y',1:4), free=FALSE, values=0:3),
mxPath(c('ib3','sb3'), arrows=2, connect="unique.pairs", values=c(1,0,1)),
mxPath('one', 'z', free=FALSE, labels="data.z"),
mxPath('z', c('ib3','sb3')),
mxPath('one', c('ib3','sb3')),
mxData(ex923[!duplicated(ex923$level3),], 'raw', primaryKey='level3')
)
level2Model <- mxModel(
'level2Model', type='RAM', level3Model,
latentVars=c(paste0('y',1:4), 'ib2', 'sb2', 'w'),
productVars=c("ib2×sb2"),
mxPath(from=c("ib2","sb2"), to="ib2×sb2", arrows=1, free=FALSE, values=1),
mxPath(from="ib2×sb2", to=paste0('y',1:4), arrows=1, free=TRUE, values=.5), #, labels="b1"),
mxPath('ib2', paste0('y',1:4), free=FALSE, values=1),
mxPath('sb2', paste0('y',1:4), free=FALSE, values=0:3),
mxPath(c('ib2','sb2'), arrows=2, connect="unique.pairs", values=c(1,0,1)),
mxPath('one', 'w', free=FALSE, labels="data.w"),
mxPath('w', c('ib2','sb2')),
mxPath(paste0('y',1:4), arrows=2),
mxPath(paste0('level3Model.y', 1:4), paste0('y',1:4), free=FALSE, values=1,
joinKey="level3"),
mxData(ex923[!duplicated(ex923$level2),], 'raw', primaryKey='level2')
)
withinModel <- mxModel(
'withinModel', type='RAM', level2Model,
manifestVars=c(paste0('y',1:4)),
latentVars=c('iw','sw', 'x'),
productVars=c("iw×sw"),
mxPath(from=c("iw","sw"), to="iw×sw", arrows=1, free=FALSE, values=1),
mxPath(from="iw×sw", to=paste0('y',1:4), arrows=1, free=TRUE, values=.5), #, labels="b1"),
mxPath('iw', paste0('y',1:4), free=FALSE, values=1),
mxPath('sw', paste0('y',1:4), free=FALSE, values=0:3),
mxPath(paste0('y',1:4), arrows=2),
mxPath(c('iw','sw'), arrows=2, connect="unique.pairs", values=c(1,0,1)),
mxPath('one', 'x', free=FALSE, labels="data.x"),
mxPath('x', c('iw','sw')),
mxPath(paste0('level2Model.y', 1:4), paste0('y',1:4), free=FALSE, values=1,
joinKey="level2"),
mxData(ex923, 'raw'),
mxFitFunctionWLS()
)
withinModel <- mxRun(withinModel)
summary(withinModel)
Data file: ex9.23.dat