For my MITER research, I have been working with a faculty member in the University of Minnesota’s Institute of Child Development. I have been involved on a couple of projects with this group. One project involved longitudinal modeling of work competency (from childhood to adulthood) using structural equation modeling (SEM). This project examined how childhood competence in externalizing behaviors/conduct, academics, and socializing behaviors predict work as an adult. If you’re interested in this work, it has just been accepted and should be in press soon. (Feel free to email me for a copy of this manuscript and I’ll see if I can get a copy).
To do the SEM, we initially tried the sem package in R. This proved to be to basic for us as it was unable to do a multi-group analysis. So we settled on MPlus. Even in MPlus, we were unable to run the models we really wanted to because the number of integration points were too many and MPlus would just error out. We dealt with it the best we could and I think our manuscript is strong and transparent about our computational problems.
After this paper was submitted, I become more interested in examining what sort of alternatives exist. Is there a Bayesian book on SEM? Well apparently there is. In the Bayesian framework, latent constructs are just dealt with as missing data problems, making them more computationally solvable. Additionally, using Bayesian you work directly from your observed data and not from a covariance matrix. Because of this you can create more highly complex, developmentally appropriate models that are infeasible in traditional software. Even Bollen is writing about this.
To help you get started. I am going to paste the code that I have used to replicate our models in WinBUGS. Unfortunately I can not share the data. But if you have specific data you can share then I’d be happy to modify the code as appropriate.
model {
for(i in 1:N){
##########################
# measurement equation model
##########################
for(j in 1:P){
y[i,j]~dnorm(mu[i,j],psi[j])
}
##########################
# alp[i] corresponds to the intercepts
##########################
mu[i,1]<-xi[i,1]+alp[1] ## Ext1
mu[i,2]<-lam[1]*xi[i,1]+alp[2]
mu[i,3]<-lam[2]*xi[i,1]+alp[3]
mu[i,4]<-xi[i,2]+alp[4] ## Soc1
mu[i,5]<-lam[3]*xi[i,2]+alp[5]
mu[i,6]<-lam[4]*xi[i,2]+alp[6]
mu[i,7]<-xi[i,3]+alp[7] ## Aca1
mu[i,8]<-lam[5]*xi[i,3]+alp[8]
mu[i,9]<-lam[6]*xi[i,3]+alp[9]
mu[i,10]<-xi[i,4]+alp[10] ## Ext2
mu[i,11]<-lam[7]*xi[i,4]+alp[11]
mu[i,12]<-lam[8]*xi[i,4]+alp[12]
mu[i,13]<-xi[i,5]+alp[13] ## Soc2
mu[i,14]<-lam[9]*xi[i,5]+alp[14]
mu[i,15]<-lam[10]*xi[i,5]+alp[15]
mu[i,16]<-xi[i,6]+alp[16] ## Aca2
mu[i,17]<-lam[11]*xi[i,6]+alp[17]
mu[i,18]<-lam[12]*xi[i,6]+alp[18]
mu[i,19]<-xi[i,7]+alp[19] ## Work2
mu[i,20]<-lam[13]*xi[i,7]+alp[20]
mu[i,21]<-lam[14]*xi[i,7]+alp[21]
##########################
# structural equation model
##########################
xi[i,1:3]~dmnorm(u[1:3],phi[1:3,1:3]) ## factor loadings are multivariate normal
xi[i,4]~dnorm(ext2[i],psd)
xi[i,5]~dnorm(soc2[i],psd)
xi[i,6]~dnorm(aca2[i],psd)
xi[i,7]~dnorm(work2[i],psd)
##########################
# continuity Paths
##########################
ext2[i]<-gam[1]*xi[i,1]
soc2[i]<-gam[2]*xi[i,2]
aca2[i]<-gam[3]*xi[i,3]
work2[i]<-gam[4]*xi[i,1] + gam[5]*xi[i,1] + gam[6]*xi[i,1]
} #end of i
##########################
# priors on intercepts
##########################
for(j in 1:21){alp[j]~dnorm(0.0, 1.0)}
##########################
# priors on loadings and coefficients
##########################
lam[1]~dnorm(0.8,psi[1])
lam[2]~dnorm(0.8,psi[2])
lam[3]~dnorm(0.8,psi[3])
lam[4]~dnorm(0.8,psi[4])
lam[5]~dnorm(0.8,psi[5])
lam[6]~dnorm(0.8,psi[6])
lam[7]~dnorm(0.8,psi[7])
lam[8]~dnorm(0.8,psi[8])
lam[9]~dnorm(0.8,psi[9])
lam[10]~dnorm(0.8,psi[10])
lam[11]~dnorm(0.8,psi[11])
lam[12]~dnorm(0.8,psi[12])
lam[13]~dnorm(0.8,psi[13])
lam[14]~dnorm(0.8,psi[14])
for(j in 1:6){gam[j]~dnorm(0.5, psd)}
##########################
#priors on precisions
##########################
for(j in 1:P){
psi[j]~dgamma(9.0, 4.0)
sgm[j]<-1/psi[j]
}
psd~dgamma(9.0, 4.0)
sgd<-1/psd
phi[1:3,1:3]~dwish(R[1:3,1:3], 5)
phx[1:3,1:3]<-inverse(phi[1:3,1:3])
} #end of model
I’d be happy to answer any questions about the code. Note this is just the model, not the starting values. I usually let WinBUGS generate them for me. Bad practice I know, but as long as my chains mix I am not worried.
