h1

SEM in WinBUGS

October 24, 2009

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.

Advertisements

48 comments

  1. Hi,
    Can you kindly explain the model and the assumptions(distributions chosen)


    • The model involves three exogenous variables and four endogenous variables. The first part of the code refers to the the confirmatory factor analysis part of the model. The next part of the model says that these factors are distributed as multivariate normal. The section after that is the actual structural component of the model, where I have the exogenous variables predicting the endogenous variables. The remaining code refers to the priors on the intercepts and factor loadings (all with weak normally distributed priors) and finally the VCV matrix (using the inverse Wishart and gamma distributions, these are conjugate priors to multivariate normal distributions). For more information about the model see the Sik-Yum Lee’s book. This model is very similar to a model he describes.


  2. Hi, I am trying to estimate similar SEM model. I have draft of the code written would you mind to look at it? I would appreciate your help very much! Thanks in advance…


    • Sure I can take a look.


      • Thanks a lot!
        Do you have email where I can send you the code and explanation? If you do not wanna publish your email here just send it to me to stochl@centrum.cz.
        Also I can publish the code here if you prefer it. The model is quite complex.
        once again, thank you very much for your kind offer.
        Jan


      • hi
        can i send to you my model too look on it please.


  3. Ok, here below is my code (you will that i was inspired by your code:-)). The item complexity for several items is more than 1 (i.e I estimate bifactorial structure). I am interested in correlation between xi[i,9] and xi[i,7] and also between xi[i,8] and xi[i,7].
    I just provide 1 line of data since they are confidential (parkinson disease patients). Sorry for that.

    Thanks a lot for your help!
    Jan

    model {
    for(i in 1:N){
    #measurement equation model
    for(j in 1:27){
    y[i,j]~ dcat(Z[1:5]) #first 27 items are categorical, but it does not work with dcat(y[i,j,1:5]) # Added categorical variable as per our previous model but not sure whether alpha0 should be specified in the data (as below) or as a hyperprior with distribution.
    }
    for(j in 28:35){
    y[i,j]~dnorm(mu[i,j],psi[j]) # last 8 items are continuous
    }

    mu[i,1]<-xi[i,1]+alp[1] #Face factor
    mu[i,2]<-lam[1]*xi[i,1]+alp[2]
    mu[i,3]<-xi[i,2]+alp[3] #Tremor factor + right+left factors
    mu[i,4]<-lam[2]*xi[i,2]+lam[23]*xi[i,6]+alp[4]
    mu[i,5]<-lam[3]*xi[i,2]+lam[24]*xi[i,7]+alp[5]
    mu[i,6]<-lam[4]*xi[i,2]+lam[25]*xi[i,6]+alp[6]
    mu[i,7]<-lam[5]*xi[i,2]+lam[26]*xi[i,7]+alp[7]
    mu[i,8]<-lam[6]*xi[i,2]+lam[27]*xi[i,6]+alp[8]
    mu[i,9]<-lam[7]*xi[i,2]+lam[28]*xi[i,7]+alp[9]
    mu[i,10]<-xi[i,3]+alp[10] #Rigidity factor + right+left factors
    mu[i,11]<-lam[8]*xi[i,3]+lam[29]*xi[i,6]+alp[11]
    mu[i,12]<-lam[9]*xi[i,3]+lam[30]*xi[i,7]+alp[12]
    mu[i,13]<-lam[10]*xi[i,3]+lam[31]*xi[i,6]+alp[13]
    mu[i,14]<-lam[11]*xi[i,3]+lam[32]*xi[i,7]+alp[14]
    mu[i,15]<-xi[i,4]+lam[33]*xi[i,6]+alp[15] #Bradykinesia factor + right+left factors mu[i,16]<-lam[12]*xi[i,4]+lam[34]*xi[i,7]+alp[16]
    mu[i,17]<-lam[13]*xi[i,4]+lam[35]*xi[i,6]+alp[17]
    mu[i,18]<-lam[14]*xi[i,4]+lam[36]*xi[i,7]+alp[18]
    mu[i,19]<-lam[15]*xi[i,4]+lam[37]*xi[i,6]+alp[19]
    mu[i,20]<-lam[16]*xi[i,4]+lam[38]*xi[i,7]+alp[20]
    mu[i,21]<-lam[17]*xi[i,4]+xi[i,6]+alp[21] #lambda not here – identification of Right handedness
    mu[i,22]<-lam[18]*xi[i,4]+xi[i,7]+alp[22]#lambda not here – identification of Left handedness
    mu[i,23]<-xi[i,5]+alp[23] #Body bradykinesia factor
    mu[i,24]<-lam[19]*xi[i,5]+alp[24]
    mu[i,25]<-lam[20]*xi[i,5]+alp[25]
    mu[i,26]<-lam[21]*xi[i,5]+alp[26]
    mu[i,27]<-lam[22]*xi[i,5]+alp[27]
    mu[i,28]<-xi[i,8]+alp[28] #Right Handedness factor
    mu[i,29]<-lam[39]*xi[i,8]+alp[29]
    mu[i,30]<-lam[40]*xi[i,8]+alp[30]
    mu[i,31]<-lam[41]*xi[i,8]+alp[31]
    mu[i,32]<-xi[i,9]+alp[32] #Left Handedness factor
    mu[i,33]<-lam[42]*xi[i,9]+alp[33]
    mu[i,34]<-lam[43]*xi[i,9]+alp[34]
    mu[i,35]<-lam[44]*xi[i,9]+alp[35]
    #structural equation model
    xi[i,1:9]~dmnorm(u[1:9],phi[1:9,1:9])
    # xi[i,9]<-gam[1]*xi[i,7] #HERE I WANTED TO SPECIFY COVARIANCE BETWEEN xi (i,9) and xi(i,7) or at least regression
    # xi[i,8]<-gam[2]*xi[i,6] #HERE I WANTED TO SPECIFY COVARIANCE BETWEEN xi (i,8) and xi(i,6) or at least regression
    }# end of i

    for(i in 1:9){u[i]<-0.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])
    lam[15]~dnorm(0.8,psi[15])
    lam[16]~dnorm(0.8,psi[16])
    lam[17]~dnorm(0.8,psi[17])
    lam[18]~dnorm(0.8,psi[18])
    lam[19]~dnorm(0.8,psi[19])
    lam[20]~dnorm(0.8,psi[20])
    lam[21]~dnorm(0.8,psi[21])
    lam[22]~dnorm(0.8,psi[22])
    lam[23]~dnorm(0.8,psi[23])
    lam[24]~dnorm(0.8,psi[24])
    lam[25]~dnorm(0.8,psi[25])
    lam[26]~dnorm(0.8,psi[26])
    lam[27]~dnorm(0.8,psi[27])
    lam[28]~dnorm(0.8,psi[28])
    lam[29]~dnorm(0.8,psi[29])
    lam[30]~dnorm(0.8,psi[30])
    lam[31]~dnorm(0.8,psi[31])
    lam[32]~dnorm(0.8,psi[32])
    lam[33]~dnorm(0.8,psi[33])
    lam[34]~dnorm(0.8,psi[34])
    lam[35]~dnorm(0.8,psi[35])
    lam[36]~dnorm(0.8,psi[36])
    lam[37]~dnorm(0.8,psi[37])
    lam[38]~dnorm(0.8,psi[38])
    lam[39]~dnorm(0.8,psi[39])
    lam[40]~dnorm(0.8,psi[40])
    lam[41]~dnorm(0.8,psi[41])
    lam[42]~dnorm(0.8,psi[42])
    lam[43]~dnorm(0.8,psi[43])
    lam[44]~dnorm(0.8,psi[44])

    for(j in 1:2){gam[j]~dnorm(0.5, psd)}

    #priors on intercepts
    for(j in 1:P){alp[j]~dnorm(0.0, 1.0)}

    #priors on precisions
    for(j in 1:P){
    psi[j]~dgamma(10,9)
    sgm[j]<-1/psi[j]
    }
    psd~dgamma(10,8)
    sgd<-1/psd
    phi[1:9,1:9]~dwish(R[1:9,1:9], 30)
    phx[1:9,1:9]<-inverse(phi[1:9,1:9])
    Z[1:5] ~ ddirch(alpha0[1:5])
    } #end of model

    Data Set
    list(N=119, P=35, R=structure(
    .Data=c(8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 0.0, 8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 0.0, 0.0, 8.0, 0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 0.0, 0.0, 0.0, 8.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 0.0, 0.0, 0.0, 0.0, 8.0, 0.0, 0.0, 0.0,
    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.0, 0.0, 0.0,
    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.0, 0.0,
    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.0),
    .Dim=c(9,9)),
    y=structure(
    .Data=c(1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,33,19,496,12,29,21,389,11,
    .
    .
    .),
    .Dim=c(119,35)))


  4. It looks good to me. Are your estimates what you’re expecting and is this model converging?


    • Thanks for looking on it. Unfortunately it cannot be even compiled! One of the problem is with priors on loadings, namely

      lam[36]~dnorm(0.8,psi[36])
      lam[37]~dnorm(0.8,psi[37])
      lam[38]~dnorm(0.8,psi[38])
      lam[39]~dnorm(0.8,psi[39])
      lam[40]~dnorm(0.8,psi[40])
      lam[41]~dnorm(0.8,psi[41])
      lam[42]~dnorm(0.8,psi[42])
      lam[43]~dnorm(0.8,psi[43])
      lam[44]~dnorm(0.8,psi[44])

      since the loop for i is from 1 till 35 only (35 items). But this bi-factor structure implies more loadings….how can I fix it?
      The other problem is how to specify covariances between factors?
      Thanks, Jan


      • Desjardins, this is a very useful thread — I have been trying to do this for a while in an ecological problem and this helps. Could you post a graphical version of this model (an OpenBugs doodle equivalent) — that would help me make sure I am understanding all the linkages right.

        Jan, it looks to me as if you have set up the multivariate normal prior on your loadings xi[,1:9] correctly. If you fit the model, you can look at the elements of the covariance matrix phi[,] that correspond to the covariances between those xi’s. Just make sure you set up the Wishart prior so it doesn’t force the off-diagonal elements of the covariance matrix to 0.

        Note there is still no latent link level corresponding to {ext2, soc2, aca2, work2} in the original model; essentially the links between the exogenous and endogenous variables as I understand it. But maybe that is OK for your application.

        To fix the loop problem you mention, it looks like you could just change P to 44 in your data set.

        Andrew


      • I will write up a graphical version in Graphviz of my model and post this in the next day or so.


      • Much obliged!


      • Click here for the SEM model


  5. Thanks a lot Andrew!
    the model is now working after the corrections you suggested! So thanks a lot! I am very grateful for your useful help!


    • That’s great, glad to hear it. Out of curiosity, did your SEM model give you the inferences you needed?
      Andrew


  6. Dear Andrew, thanks a lot! The model is working after the corrections you suggested! I am very grateful for your useful help! Thank you!!!

    Jan


  7. Hi, I’m trying to solve SEM in winbugs. I’m using simulation data to build SEM model with 4 observed and 2 latent endogenous variables and 2 observed and 1 latent exogenous variables. And model specified correctly and data loaded. But it gives error (expected multivariate node) when model is compiling. Can you help me to solve this problem?
    And I wanna constrain some parameters like psi(3)=psi(4) and psi(5)=psi(6). How can I constrain these parameters? Thank you so much.

    model {
    for(i in 1:N){
    #measurement equation model
    for(j in 1:P){
    y[i,j]~dnorm(mu[i,j],psi[j])
    ephat[i,j]<-y[i,j]-mu[i,j]
    }
    mu[i,1]<-eta[i,1]
    mu[i,2]<-lam[1]*eta[i,1]
    mu[i,3]<-eta[i,2]
    mu[i,4]<-lam[2]*eta[i,2]
    mu[i,5]<-xi[i]
    mu[i,6]<-lam[3]*xi[i]

    #structural equation model
    xi[i]~dnorm(u[i],phi)
    eta[i,1:2]~dmnorm(nu[i,1:2],psd[1:2,1:2])
    #eta[i,1]~dnorm(nu[i,1],psd[1,1])
    #eta[i,2]~dnorm(nu[i,2],psd[2,2])
    nu[i,1]<-gam[1]*xi[i]
    nu[i,2]<-gam[2]*xi[i]+beta[1]*nu[i,1]
    dthat[i,1:2]<-eta[i,1:2]-nu[i,1:2]
    } #end of i

    #priors on loadings and coefficients
    lam[1]~dnorm(0.8,psi[2])
    lam[2]~dnorm(0.8,psi[4])
    lam[3]~dnorm(0.8,psi[6])
    beta[1]~dnorm(0.5,psd[2,2])
    gam[1]~dnorm(0.5, psd[1,1])
    gam[2]~dnorm(0.5,psd[2,2])

    #priors on precisions
    for(j in 1:P){
    psi[j]~dgamma(9.0, 4.0)
    sgm[j]<-1/psi[j]
    }
    psd[1:2,1:2]~dwish(R[1:2,1:2], 4.0)
    sgd[1:2,1:2]<-inverse(psd[1:2,1:2])
    phi~dgamma(9.0, 5.0)
    phx<-1/phi

    } #end of model

    Data
    list(N=100, P=6, u=c(0,0),
    y=structure(simdata),
    .Dim=c(100,6)),
    R=structure(
    .Data=c(1.0),.Dim=c(1,1)))
    initial value
    list(lam=c(0.0, 0.0, 0.0),
    psi=c(1.0, 1.0, 1.0, 1.0, 1.0, 1.0),
    psd=structure(
    .Data=c(0.2, 0.1,
    0.1, 0.9),
    .Dim=c(2,2)),
    phi=1.0,
    gam=c(1.0, 1.0))


  8. Is it possible to fit a multilevel Factor Analysis in WinBUGS?


    • I have never done it. But I would imagine it is definitely possible in WinBUGS.


  9. I am not quite conversant with this.
    My problem:
    I have 30 indicators of quality observed for about 500 health workers in 8 hospitals. the data were collected every six months for a period of 2.5 years.
    I would like use factor analysis or some similar method to group this indicators. the idea being that measuring one of such indicators in group would be sufficient. I postulate that the indicators will fall in 3 main groups: assessment, diagnosis and treatment.
    I also have alot of missing data as the indicators span 3 conditions. Thus for clinicians who did not see a patient with a certain diagnosis, they have a missing observation.
    Do you think the analysis is appropriate? Is SEM the way to go?


    • Hi Stephen,
      It certainly sounds like you could do a factor analysis with your data and you might run an EFA or CFA depending on your certainty with which you expect what indicators to fall into which group. I might also suggest that you consider doing a principal components analysis if the purpose is to reduce the number of indicators. So you might want to look into that. What is the ultimate goal of your analysis? Are you planning on doing latent multilevel modeling?


  10. Hi,
    I am trying to write winbugs codes for a latent variable. Could you please help me a little bit?The model is like this:
    p1 is an unobserved latent variable, y1 is an endogenous variable (for y2) and e1 is the residual that enters directly into equation for y2 to control the influence from latent variable p1. So b2 captures the effect of y1 on y2 when controlling p1. y1, y2 are observed.

    y1=a0+a1*p1+e1

    y2=b0+b1*p1+b2*e1

    Could you please give me some hints on how to model such a latent variable model in winbugs?

    Thanks very much


    • Hi Alice,
      I will quickly look at this tomorrow and let you know my thoughts.
      Chris


      • Hi,
        Do you have a chance to look at the model? I know you must be very busy and it takes time to write the code. I really appreciate your time and help on this! I am looking forward to hearing from you!

        Thanks a lot!
        Alice


      • Hi Alice,

        First, you’re going to have issues with model identification if you only have two manifest variables for your one latent variable. You need at least 3 manifest variables, see Mulaik 2010 Foundations of Factor Analysis.

        Second, you need to supply me with some data if you want me to reproduce your problem.

        Third, what have you tried already? Can you provide me with your syntax. Asking me to write your model for you with out specification of what you want your priors to be, etc., is a lot of work and can result in a very bad model.


      • Hi Chris,

        Here is my full model and code: I have 6 observed endogenous variables (y_1, y_2, y_3, y_4,y_5, y_6),1 latent (unobserved) exogenous variable p, and 4 other observed exogenous variables (x1,x2,x3,x4). Could you please help me check if the code fits the model? If you need more infor, please let me know.Thanks very much. –Alice

        I use latent variable p to estimate y_1, y_2, y_3:

        y_j= π_0j+π_1j*p + δ_j , j=1,2,3, δ_j is residuals

        I use several instrumental variables to estimate y_4 and y_5 and get instruments for y_4 and y_5:

        y_4inst= φ_10+ φ_11*latent+ φ_12*δ_1+φ_13*δ_2+φ_14*δ_3+φ_15*x1+φ_16*x2

        y_5inst= φ_20+ φ_21*latent+ φ_22*δ_1+φ_23*δ_2+φ_24*δ_3+φ_25*x1+φ_26*x2

        y_6:

        y_6= β_0+β_1*latent + β_2*δ_1+β_3*δ_2+β_4*δ_3+ β_5*y_4inst+ β_6*y_5inst+β_7*x3+β_8*x1+β_9*x4+ε

        Here is my winbugs code:
        model
        {
        for(i in 1: N) {

        # get y_1, y_2, y_3, y_4, y_5

        Y[1,i] <- y1[i]
        Y[2,i] <- y2[i]
        Y[3,i] <- y3[i]
        Y[4,i] <- y4[i]
        Y[5,i] <- y5[i]

        y6[i] ~ dnorm(u_y6[i], tau[1])

        # latent variable
        p[i] ~ dnorm(0, taul)

        #observed endogenous variable y_1, y_2, y_3, y_4 and y_5

        for (j in 1:5) { Y[j,i] ~ dnorm (Yinst[j,i], taun[j])}

        # get instruments for y_1, y_2 and y_3 and residuals
        for (j in 1:3)
        {Yinst[j,i] <- pi[j,1]+pi[j,2]*p[i]
        res[j,i] <- Y[j,i]-Yinst[j,i]
        }

        # instruments for y_4 and y_5

        Yinst[4,i] <- phi[1,1] + phi[1,2]*p[i]+phi[1,3]*res[1,i] +phi[1,4]*res[2,i] +phi[1,5]*res[3,i] + phi[1,6]*x1[i] +phi[1,7]*x2[i]

        Yinst[5,i] <- phi[2,1] + phi[2,2]*p[i]+phi[2,3]*res[1,i] +phi[2,4]*res[2,i] +phi[2,5]*res[3,i] + phi[2,6]*x1[i] +phi[2,7]*x2[i]

        # stage 2

        u_y6[i]<- beta[1]+beta[2]*p[i]+beta[3]*res[1,i]+beta[4]*res[2,i]+beta[5]*res[3,i] + beta[6]*Yinst[4,i]+ beta[7]*Yinst[5,i]+ beta[8]*x3[i] +beta[9]*x1[i] + beta[10]*x4[i]
        }

        # priors on regression parameters and precisions

        for (j in 1:1) {tau[j] ~ dgamma(0.01, 0.01)}

        for (j in 1:3) { pi[j,1] ~ dnorm(0, 1.0E-6)
        pi[j,2] ~ dnorm(0,1.0E-3)}
        for (j in 1:5) {taun[j] ~ dgamma(1, 1.0E-3)}
        taul ~ dgamma(0.01, 0.01)

        for (j in 1:2) { phi[j,1] ~ dnorm(0, 1.0E-6)
        phi[j,2] ~ dnorm(0,1.0E-3)
        phi[j,3] ~ dnorm(0,1.0E-3)
        phi[j,4] ~ dnorm(0,1.0E-3)
        phi[j,5] ~ dnorm(0,1.0E-3)
        phi[j,6] ~ dnorm(0,1.0E-3)
        phi[j,7] ~ dnorm(0,1.0E-3)}

        for (j in 1:10) {beta[j] ~ dnorm(0, 0.001)}

        }


      • I recommend you post your model to the BUGS mailing list. I don’t have time to look at the code. Sorry.


  11. Well, if you don’t have time to look at it, tell me directly at the beginning and I would be perfectly OK with that. Why did you post all those comments to make me feel that you would have helped me and looked at it? It just wasted both of our time.


    • While I understand your disappointment, I had intended on spending time reviewing your code and helping you. I did not realize that I would be so busy initially and I am now working on 6 projects concurrently. The fact that I referred you to a more definitive source on Bayesian analysis shouldn’t be received with overt hostility. The fact that you can easily copy/paste code and questions between these media should make this a non-issue. This blog is a hobby intended to help others. I am not your advisor or employer and thus have no obligation to helping you. Furthermore, I am not sure with your attitude how much success you’ll receive.


  12. I have many questions about my code. Syntactically its correct but when i choose load init or gen inits in winBUGS I always got message “Trap”. Can you help me please please please, how i can solve this problem? This is the code and the data:

    model{
    for(i in 1:N){
    #measurement equation model
    for(j in 1:P){
    y[i,j]~dnorm(mu[i,j],psi[j])I(thd[j,z[i,j]],thd[j,z[i,j]+1])
    ephat[i,j]<-y[i,j]-mu[i,j]
    }
    mu[i,1]<-eta5[i] # mu eta5
    mu[i,2]<-lam[1]*eta5[i]

    mu[i,3]<-eta4[i] # mu eta4
    mu[i,4]<-lam[2]*eta4[i]
    mu[i,5]<-lam[3]*eta4[i]

    mu[i,6]<-eta3[i] # mu eta3
    mu[i,7]<-lam[4]*eta3[i]
    mu[i,8]<-lam[5]*eta3[i]

    mu[i,9]<-eta2[i] # mu eta2
    mu[i,10]<-lam[6]*eta2[i]
    mu[i,11]<-lam[7]*eta2[i]
    mu[i,12]<-lam[8]*eta2[i]
    mu[i,13]<-lam[9]*eta2[i]

    mu[i,14]<-eta1[i] # mu eta1
    mu[i,15]<-lam[10]*eta1[i]
    mu[i,16]<-lam[11]*eta1[i]
    mu[i,17]<-lam[12]*eta1[i]

    mu[i,18]<-xi[i,1] # mu xi1
    mu[i,19]<-lam[13]*xi[i,1]
    mu[i,20]<-lam[14]*xi[i,1]
    mu[i,21]<-lam[15]*xi[i,1]
    mu[i,22]<-lam[16]*xi[i,1]
    mu[i,23]<-lam[17]*xi[i,2]
    mu[i,24]<-lam[18]*xi[i,2]
    mu[i,25]<-lam[19]*xi[i,2]
    mu[i,26]<-lam[20]*xi[i,2]

    mu[i,27]<-xi[i,2] #mu xi2
    mu[i,28]<-lam[21]*xi[i,2]
    mu[i,29]<-lam[22]*xi[i,2]
    mu[i,30]<-lam[23]*xi[i,2]

    mu[i,31]<-xi[i,3] #mu xi3
    mu[i,32]<-lam[24]*xi[i,3]
    mu[i,33]<-lam[25]*xi[i,3]

    #structural equation model
    xi[i,1:3]~dmnorm(u[1:3],phi[1:3,1:3])

    eta1[i]~dnorm(nu1[i],psd1)
    nu1[i]<-gam[1]*xi[i,1]+gam[3]*xi[i,2]+gam[4]*xi[i,3]
    dthat1[i]<-eta1[i]-nu1[i]

    eta2[i]~dnorm(nu2[i],psd2)
    nu2[i]<-gam[2]*xi[i,1]+gam[5]*xi[i,3]+beta[1]*eta1[i]
    dthat2[i]<-eta2[i]-nu2[i]

    eta3[i]~dnorm(nu3[i],psd3)
    nu3[i]<-beta[2]*eta1[i]+beta[3]*eta2[i]
    dthat3[i]<-eta3[i]-nu3[i]

    eta4[i]~dnorm(nu4[i],psd4)
    nu4[i]<-beta[4]*eta2[i]+beta[5]*eta3[i]
    dthat4[i]<-eta4[i]-nu4[i]

    eta5[i]~dnorm(nu5[i],psd5)
    nu5[i]<-beta[6]*eta4[i]
    dthat5[i]<-eta5[i]-nu5[i]

    }# end of i

    for(i in 1:3){u[i]<-0.0}

    #priors on loadings and coefficients
    #lamda
    var.lam[1]<-4.0*psi[2]

    var.lam[2]<-4.0*psi[4]
    var.lam[3]<-4.0*psi[5]

    var.lam[4]<-4.0*psi[7]
    var.lam[5]<-4.0*psi[8]

    var.lam[6]<-4.0*psi[10]
    var.lam[7]<-4.0*psi[11]
    var.lam[8]<-4.0*psi[12]
    var.lam[9]<-4.0*psi[13]

    var.lam[10]<-4.0*psi[15]
    var.lam[11]<-4.0*psi[16]
    var.lam[12]<-4.0*psi[17]

    var.lam[13]<-4.0*psi[19]
    var.lam[14]<-4.0*psi[20]
    var.lam[15]<-4.0*psi[22]
    var.lam[16]<-4.0*psi[22]
    var.lam[17]<-4.0*psi[23]
    var.lam[18]<-4.0*psi[24]
    var.lam[19]<-4.0*psi[25]
    var.lam[20]<-4.0*psi[26]

    var.lam[21]<-4.0*psi[28]
    var.lam[22]<-4.0*psi[29]
    var.lam[23]<-4.0*psi[30]

    var.lam[24]<-4.0*psi[32]
    var.lam[25]<-4.0*psi[33]

    for(i in 1:25){lam[i]~dnorm(0.8,var.lam[i])}

    #Gamma
    var.gam1<-4.0*psd1
    var.gam2<-4.0*psd2
    var.gam3<-4.0*psd3
    var.gam4<-4.0*psd4
    var.gam5<-4.0*psd5

    gam[1]~dnorm(0.6,var.gam1)
    gam[2]~dnorm(0.6,var.gam2)
    gam[3]~dnorm(0.4,var.gam3)
    gam[4]~dnorm(0.4,var.gam4)
    gam[5]~dnorm(0.4,var.gam5)

    #Beta
    var.bet1<-2*psd6
    beta[1]~dnorm(0.8,var.bet1)
    var.bet2<-2*psd7
    beta[2]~dnorm(0.4,var.bet2)
    var.bet3<-2*psd8
    beta[3]~dnorm(0.3,var.bet3)
    var.bet4<-2*psd9
    beta[4]~dnorm(0.3,var.bet4)
    var.bet5<-2*psd10
    beta[5]~dnorm(0.3,var.bet5)
    var.bet6<-2*psd11
    beta[6]~dnorm(0.3,var.bet6)

    #priors on precisions
    for(j in 1:P){
    psi[j]~dgamma(1,2)
    sgm[j]<-1/psi[j]
    }
    psd1~dgamma(1,2)
    psd2~dgamma(1,2)
    psd3~dgamma(1,2)
    psd4~dgamma(1,2)
    psd5~dgamma(1,2)

    psd6~dgamma(1,2)
    psd7~dgamma(1,2)
    psd8~dgamma(1,2)
    psd9~dgamma(1,2)
    psd10~dgamma(1,2)
    psd11~dgamma(1,2)

    sgd1<-1/psd1
    sgd2<-1/psd2
    sgd3<-1/psd3
    sgd4<-1/psd4
    sgd5<-1/psd5

    sgd6<-1/psd6
    sgd7<-1/psd2
    sgd8<-1/psd3
    sgd9<-1/psd4
    sgd10<-1/psd5
    sgd11<-1/psd6
    phi[1:3,1:3]~dwish(R[1:3,1:3], 30)
    phx[1:3,1:3]<-inverse(phi[1:3,1:3])

    } #end of model

    Data Set
    list(N=93, P=33,
    R=structure(
    .Data=c(8.0, 0.0, 0.0,
    0.0, 8.0, 0.0,
    0.0, 0.0, 8.0),
    .Dim=c(3,3)),
    thd=structure(
    .Data=c(-40,-1.51792916,-1.437635037,-0.315096782,0.521310428,40,
    -40,-1.51792916,-1.437635037,-0.315096782,0.521310428,40,
    -40,-1.609409261,-0.904762168,0.946122671,1.365668579,40,
    -40,-1.437635037,-0.752728794,0.904762168,1.18363106,40,
    -40,-1.609409261,-0.789007665,0.789007665,0.946122671,40,
    -40,-2.023605558,-1.365668579,0.616394157,0.989168627,40,
    -40,-2.298992304,-1.51792916,0.286893917,0.616394157,40,
    -40,-1.716768353,-0.826356044,0.946122671,1.300153433,40,
    -40,-2.298992304,-1.437635037,0.430727299,0.584119695,40,
    -40,-2.298992304,-1.365668579,0.343552575,0.521310428,40,
    -40,-1.716768353,-1.081286046,0.460494539,0.682973392,40,
    -40,-1.848596288,-1.51792916,0.490675589,0.789007665,40,
    -40,-1.848596288,-1.365668579,0.343552575,0.682973392,40,
    -40,-1.239788147,-0.752728794,0.946122671,1.300153433,40,
    -40,-1.300153433,-0.864894359,0.864894359,1.239788147,40,
    -40,-1.081286046,-0.789007665,0.989168627,1.300153433,40,
    -40,-1.18363106,-0.752728794,0.826356044,1.18363106,40,
    -40,-2.023605558,-0.752728794,1.034130266,1.365668579,40,
    -40,-1.51792916,-0.989168627,0.521310428,0.717414656,40,
    -40,-0.826356044,0.040440509,0.989168627,1.239788147,40,
    -40,-2.023605558,-0.864894359,0.460494539,0.826356044,40,
    -40,-0.946122671,0.040440509,1.034130266,1.300153433,40,
    -40,-1.716768353,-0.864894359,0.989168627,1.300153433,40,
    -40,-1.130977608,-0.343552575,0.789007665,0.989168627,40,
    -40,-2.298992304,-1.239788147,0.148788621,0.460494539,40,
    -40,-2.023605558,-1.300153433,0.37228936,0.490675589,40,
    -40,-0.752728794,-0.286893917,1.239788147,2.298992304,40,
    -40,-0.752728794,-0.401336977,1.130977608,1.609409261,40,
    -40,-1.034130266,-0.584119695,1.18363106,1.848596288,40,
    -40,-1.51792916,-1.081286046,0.717414656,1.081286046,40,
    -40,-1.848596288,-0.682973392,1.034130266,1.365668579,40,
    -40,-0.717414656,-0.067433552,1.18363106,1.239788147,40,
    -40,-1.609409261,-0.789007665,0.752728794,0.946122671,40),

    .Dim=c(33,6)),
    z=structure(
    .Data=c(3,3,3,3,4,3,3,3,3,4,4,4,4,2,2,3,2,2,3,1,4,3,3,2,4,5,3,2,3,4,3,2,2,
    5,5,5,5,3,5,5,5,5,5,5,5,5,5,5,3,3,5,5,3,5,5,3,5,5,5,5,5,5,5,5,5,5,
    4,4,3,4,3,3,4,3,4,4,5,3,4,3,3,4,3,3,3,2,4,4,3,3,4,3,2,2,3,2,2,3,3,
    3,3,4,3,3,3,4,4,3,3,3,3,4,2,3,4,4,4,2,2,4,2,3,4,4,4,4,4,3,4,3,3,3,
    4,4,4,4,4,4,4,3,5,5,5,5,4,4,4,3,3,3,3,2,4,3,3,3,3,3,4,4,4,4,4,3,3,
    4,4,3,3,3,3,5,3,5,5,5,5,5,3,3,3,3,3,3,1,5,1,3,3,5,5,3,3,3,3,3,1,3,
    4,4,3,3,3,4,4,3,4,4,4,4,3,4,4,4,4,2,4,1,2,2,4,4,4,4,4,4,4,4,4,1,4,
    3,3,5,3,4,4,4,5,3,4,4,3,4,2,3,3,4,5,5,5,5,5,4,5,5,5,3,4,4,4,3,2,3,
    3,3,2,1,2,2,3,2,2,2,2,3,2,3,3,3,2,3,3,2,3,2,2,3,3,3,2,2,2,2,3,2,2,
    5,5,4,5,5,3,4,4,5,5,4,4,4,4,4,4,3,3,5,2,4,3,4,3,5,5,4,4,4,4,5,4,5,
    1,1,1,1,2,1,3,2,1,2,1,2,1,3,4,3,3,2,2,2,1,1,3,2,4,5,3,3,4,5,2,1,2,
    3,3,3,4,3,3,2,2,4,3,3,4,3,4,3,3,3,2,1,2,1,1,3,3,2,1,3,3,4,5,2,3,5,
    4,4,5,3,5,5,5,3,5,5,5,5,5,3,5,3,5,5,5,3,5,1,3,3,5,5,3,3,3,3,5,3,5,
    3,3,3,3,3,3,3,3,3,3,2,3,3,3,3,3,3,3,3,3,3,2,3,3,3,2,2,2,3,3,3,2,3,
    5,5,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,3,3,3,1,3,3,3,1,3,
    3,3,3,2,2,3,3,3,3,3,2,3,2,3,3,3,3,2,3,2,3,2,3,3,3,3,2,2,2,3,2,3,3,
    3,3,3,3,3,4,4,3,3,3,2,3,3,3,3,3,3,3,3,3,3,3,3,2,3,3,2,2,2,3,2,1,3,
    5,5,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,3,
    5,5,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,3,3,3,1,3,3,3,1,3,
    5,5,3,3,3,3,3,3,3,3,3,3,3,2,3,3,3,2,2,3,3,3,3,1,3,3,2,3,3,3,3,1,3,
    5,5,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,3,3,3,3,1,3,3,2,3,3,3,3,3,3,
    1,1,3,3,3,2,2,3,3,2,3,3,3,3,3,1,1,3,3,1,2,1,2,3,1,1,1,1,1,1,3,1,3,
    3,3,3,3,4,3,5,4,3,3,3,4,4,4,4,3,3,3,5,2,2,3,3,2,4,4,4,4,4,4,5,5,5,
    3,3,5,3,5,3,5,3,3,3,3,3,3,1,1,1,1,3,5,3,3,3,5,3,5,5,4,3,3,3,2,3,4,
    3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,3,1,1,3,3,3,3,1,3,3,1,1,1,
    3,3,3,5,5,3,5,5,5,5,5,5,5,5,5,3,5,4,4,3,3,3,2,4,5,3,3,3,3,3,3,5,5,
    3,3,5,3,5,3,5,5,5,5,3,5,5,5,5,3,3,4,5,4,4,4,4,5,4,5,4,3,3,3,4,5,4,
    3,3,4,5,5,5,5,4,4,5,4,4,4,5,3,5,4,4,5,4,3,2,4,4,5,5,4,3,3,3,4,5,5,
    3,3,4,4,5,5,4,4,5,5,5,5,5,5,4,4,4,5,3,5,4,3,2,5,5,5,3,3,3,5,4,5,5,
    4,4,3,5,5,5,4,3,5,5,5,3,5,3,3,5,3,3,5,1,5,1,3,3,5,5,3,3,3,3,3,1,5,
    3,3,3,5,5,3,3,3,5,4,5,4,5,3,3,5,4,3,5,3,3,3,5,5,5,5,1,3,3,3,3,3,5,
    3,3,3,3,3,3,3,5,5,5,5,5,5,5,5,5,5,3,3,5,5,5,5,5,5,5,3,3,3,5,3,5,5,
    5,5,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,3,3,2,5,5,5,3,3,3,5,4,5,5,
    3,3,3,3,1,1,3,3,3,3,3,3,3,3,1,1,1,3,3,1,3,1,3,1,3,3,1,1,1,3,3,1,3,
    1,1,2,4,3,3,3,4,3,1,1,1,1,3,3,3,3,3,3,1,3,1,4,1,5,3,1,3,3,5,3,1,5,
    5,5,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,3,1,3,3,3,5,3,3,3,3,3,1,3,
    4,4,3,3,5,5,5,5,5,5,5,5,5,3,3,3,3,3,3,3,5,1,3,3,5,5,3,3,3,5,3,3,3,
    4,4,3,3,3,3,3,3,3,5,5,4,4,2,3,3,3,3,3,3,3,2,3,2,3,3,3,2,3,3,3,2,2,
    5,5,4,3,5,4,4,4,5,5,5,5,5,3,4,4,4,2,4,2,4,3,3,4,3,3,3,4,3,4,3,2,4,
    1,1,1,1,1,3,3,1,3,3,3,3,3,3,3,3,3,3,2,1,3,1,3,2,3,3,3,3,3,3,1,1,1,
    4,4,3,3,5,3,5,3,3,3,3,3,3,3,3,3,3,3,5,5,3,3,2,3,4,3,3,3,5,3,3,3,3,
    4,4,2,2,2,3,3,2,3,3,3,3,3,3,2,3,2,3,3,3,3,2,3,3,3,3,3,3,2,3,3,3,3,
    3,3,1,1,1,3,3,1,3,2,2,3,3,1,1,1,1,2,3,1,2,2,3,2,2,2,1,1,1,1,3,3,3,
    4,4,1,1,2,2,2,2,2,3,3,3,3,1,1,1,1,2,3,2,2,2,3,2,2,2,1,1,1,2,2,2,2,
    3,3,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,
    4,4,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,
    3,3,2,2,2,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,2,2,2,2,
    4,4,3,3,3,3,3,2,3,3,3,3,3,2,2,2,2,3,2,2,3,2,2,2,2,3,3,3,3,3,3,3,3,
    5,5,4,5,3,5,5,2,5,4,3,3,4,3,2,3,3,5,4,4,5,3,5,2,4,3,4,5,5,5,2,3,3,
    5,5,3,3,5,3,5,2,3,5,5,5,5,3,3,2,3,3,3,1,5,2,5,2,5,5,3,3,3,5,2,2,3,
    5,5,3,3,3,5,5,3,3,3,3,3,3,1,1,1,1,3,3,3,3,1,3,3,3,3,1,3,1,3,3,3,3,
    5,5,3,3,3,3,3,3,3,3,3,3,2,1,1,1,1,3,3,5,5,5,3,5,3,3,1,1,1,1,2,2,2,
    5,5,3,3,3,3,3,3,5,5,5,3,3,1,1,1,5,5,5,5,5,5,5,5,5,5,1,1,2,2,3,3,3,
    5,5,2,2,2,3,3,3,5,5,5,5,5,2,2,1,2,2,3,3,3,3,3,5,3,3,1,1,1,1,3,3,3,
    4,4,2,2,2,3,3,2,3,3,3,3,3,3,3,3,3,3,4,4,4,4,3,3,3,3,1,1,1,3,2,2,2,
    4,4,3,3,3,5,5,3,5,5,3,5,5,5,3,3,3,3,5,3,2,3,5,5,3,5,3,3,3,3,3,3,3,
    4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
    4,4,3,3,3,3,1,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,3,3,3,3,3,3,
    2,2,3,3,3,3,3,3,5,5,5,5,5,3,3,3,3,3,3,2,3,2,2,3,5,3,2,2,2,3,3,3,3,
    4,4,2,2,3,3,3,3,3,3,3,3,3,3,3,1,3,3,5,3,3,2,3,2,3,3,1,2,2,3,3,2,2,
    4,4,2,2,2,4,3,3,3,3,3,3,3,2,3,2,3,4,4,4,4,4,3,3,3,3,3,3,3,3,3,2,3,
    3,3,3,3,3,5,5,3,5,3,5,3,3,3,3,3,3,3,3,1,3,3,3,3,3,3,3,3,3,3,3,3,3,
    3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,3,3,2,3,2,3,5,5,5,3,3,3,3,3,3,3,
    5,5,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,5,1,3,2,3,3,3,3,3,3,3,3,3,3,3,
    5,5,3,3,3,3,3,1,3,3,3,3,3,3,3,3,3,1,1,2,3,2,1,3,3,3,3,3,3,3,3,1,2,
    3,3,3,2,2,2,3,3,3,3,3,3,3,3,3,2,3,1,1,1,2,2,1,1,5,5,3,3,3,3,2,3,2,
    3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,3,2,3,3,2,3,1,3,3,3,3,
    5,5,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,3,3,3,3,3,3,2,2,3,2,3,
    4,4,3,3,1,5,5,3,5,5,5,5,5,3,5,5,5,3,5,3,5,5,5,5,5,5,3,5,3,3,5,5,3,
    4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,5,
    5,5,5,5,5,3,5,3,5,5,3,5,5,3,5,3,5,3,3,1,5,3,1,5,5,5,3,3,3,5,5,3,5,
    4,4,3,3,3,3,3,2,3,3,3,3,3,2,2,2,2,3,5,3,5,2,3,3,5,3,2,1,2,3,2,1,3,
    4,4,3,2,3,4,3,3,3,2,3,3,3,3,2,3,2,2,1,1,3,2,3,1,3,3,2,1,2,3,3,1,3,
    4,4,3,2,3,4,3,3,3,3,3,3,3,3,3,3,3,2,1,1,3,2,3,1,3,3,2,1,2,3,3,2,3,
    3,3,5,5,3,5,5,5,5,5,3,3,3,3,3,5,5,3,3,3,3,2,3,2,5,5,3,1,2,3,3,3,3,
    3,3,3,3,3,5,5,3,5,5,4,4,4,4,3,3,3,3,3,3,5,2,5,3,3,3,3,3,3,3,3,3,3,
    4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,3,2,3,3,3,3,3,3,3,3,2,2,3,
    5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,3,5,5,5,3,5,3,5,5,1,1,
    5,5,2,2,2,3,5,3,3,5,3,3,5,3,3,5,5,3,3,5,3,5,3,3,5,3,3,3,3,3,3,5,3,
    4,4,3,3,3,3,5,3,3,3,5,3,5,3,5,3,5,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
    5,5,3,3,3,4,4,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
    5,5,2,2,2,3,3,3,2,3,2,1,3,3,3,3,3,3,3,2,3,2,3,1,3,3,2,1,3,3,3,3,2,
    5,5,3,3,3,3,3,3,3,3,1,1,3,3,3,3,3,3,2,2,2,2,3,2,2,2,2,2,3,3,2,2,1,
    1,1,3,3,3,4,5,3,4,5,5,5,3,3,3,3,3,3,5,2,2,1,2,2,5,5,1,2,3,3,3,3,3,
    1,1,2,3,2,3,3,2,2,3,2,3,3,2,3,3,2,3,3,2,2,1,3,1,5,5,1,1,1,2,3,3,3,
    5,5,3,3,3,3,3,2,3,3,3,3,3,3,3,3,3,2,3,1,3,1,3,2,3,3,3,3,3,3,3,3,3,
    5,5,3,2,3,3,3,3,5,5,4,5,5,4,4,3,4,4,5,2,5,4,3,2,5,4,2,5,4,4,5,2,4,
    4,4,3,3,3,3,3,3,3,3,3,3,3,1,3,1,1,3,3,1,3,2,2,2,3,3,3,1,3,3,1,1,3,
    4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,3,2,2,3,2,1,2,3,1,3,3,3,3,1,1,
    4,4,3,5,3,3,5,3,3,3,3,3,5,3,3,3,3,3,5,5,5,5,3,3,5,5,3,3,3,5,3,3,3,
    3,3,3,3,2,3,3,3,3,3,3,3,3,3,3,3,3,2,3,2,4,3,3,3,4,3,3,3,3,3,2,2,3,
    4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,3,2,5,3,3,3,4,3,1,3,3,3,2,1,2,
    5,5,3,1,3,3,3,1,3,3,1,3,1,1,3,3,1,3,3,3,3,3,3,3,3,3,1,3,1,3,3,2,5),
    .Dim=c(93,33)))

    Two different Initial Values
    list(
    lam=c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
    0.5, 0.5, 0.5, 0.5, 0.5),
    psi=c(1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
    1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0),
    psd1=0.5, psd2=0.5, psd3=0.5, psd4=0.5, psd5=0.5, psd6=0.5
    gam=c(1.0, 1.0, 1.0, 1.0, 1.0),
    beta=c(1.0, 1.0, 1.0, 1.0, 1.0, 1.0),
    phi=structure(
    .Data=c(1.0, 0.0, 0.0,
    0.0, 1.0, 0.0,
    0.0, 0.0, 1.0),
    .Dim=c(3,3)),
    xi=structure(
    .Data=c(0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0,
    0,0,0),
    .Dim=c(93,3)))


  13. Hi can I ask for the winbugs code for confirmatory factor analysis? Can we use the same SEM codes? here’s the data:
    5 varibales
    3 factors (Finance, Marketing, Policy)

    Finance Marketing Policy
    1 3 6 5
    2 7 3 3
    3 10 9 8
    4 3 9 7
    5 10 6 5


    • I tried to replicate your code. However, I have a problem.
      The winbugs programs says I need to define xi[i,1:3]~dmnorm(u[1:3] somehow. I didn’t see it in the model,

      I just wonder how you define xi[i,1:3]~dmnorm(u[1:3]?

      then did you define xi[i,1:3]~dmnorm(u[1:3] in the list?
      just let me know


  14. Hi, i’m trying to write winbugs code for a sem model with ordinal variables (2 latent factors and 8 observed variables -measurement model-, and 17 regressors)…only 2 regressors are continue.
    Is it possible to work with this software using ordinal variables? What do you recommend?

    Many thanks


  15. One thing that would be helpful is if you could include a path diagram for the code above – it would help me translate variable and parameter names to the code itself.


  16. Dear Desjardins

    I am using similar code like the one you posted .I am not sure of one thing though and I wanted to check with you. What do we do when we have one indicator variable that we want to include in a structure equation?..i.e. we have observed variables and also few latent ones(with multiple indicators) in our data. For the latent variables with multiple indicators it is clear we fix one of the factor loading in the measurement equation for model indientifaibility. But if we want to use the observed variables in the structural equation, we make them a one indicator variable and we are not sure what to fix when we write the measurement equation. I was inclined to proceed as usual by fixing the single factor loading as 1.
    For example if the measurement equation of a 3 indicator latent variable is given us shown below,

    mu[i,1]<-eta[i,1]+alp[1]
    mu[i,2]<-lam[1]*eta[i,1]+alp[2]
    mu[i,3]<-lam[2]*eta[i,1]+alp[3]

    Then the measurement equation for one indicator would be
    mu[i,4]<-eta[i,4]+alp[4]

    However,Someone cautioned me probably I need to fix error terms as well(as is done in LISREL).But I am not sure which one to fix in the code ? is it the psi[j] ? Can you please comment on the problem?
    Thanks in advance!

    With best regards

    l.mesfin


  17. dear sir
    i have question on previous commands, my question can i use correlation matrix in the list of data without depending on raw data in
    confirmatory factor analysis.
    regards


  18. dear sir
    i want to ask you about using CFA in winbugs .my question can i put covariance matrix in data not raw data??
    regards


  19. hi, may I ask the syntax of Gelman and Rubin Diagnostic to test Convergence and Posterior Predictive p-value to assess the goodness-of-fit of SEM Bayesian?
    Thank…


  20. hi
    i have a problem with winbugs program i can not doing anything because this phrase “no monitor set’ please help me .
    regards


  21. Reblogged this on thanoonalshakerchy80's Blog and commented:
    hi
    is there anyone help me i cannot post here in this website.
    regards


  22. i have a problem with my winbugs codes of sem can you help me to solve this problem so can i send to you my winbugs codes to your e.mail please.
    regards


  23. Dear Desjardins,

    I am using WinBUGS under R and i woold like to change this loop :
    for(i in 1:N){

    #measurement equation model

    for (j in 1:P){
    if(j==1 || j==3 || j==4 || j==5 || j>=10 )
    {
    y[i,j]~dnorm(mu[i,j],psi[j])*I(low[z[i,j]+1],high[z[i,j]+1])
    }
    else
    {
    y[i,j]~dnorm(mu[i,j],psi[j])
    ephat[i,j]<-y[i,j]-mu[i,j]
    }

    }

    because it is not accepted by WinBUGS ! Please do you any idea how it can be written in WinBUGS ?

    Thank you for your help

    Sincerely yours,
    jade


  24. Hello everyone, I am sorry but there is a mistake in the loop :

    for(i in 1:N){

    #measurement equation model

    for (j in 1:P){
    if(j==1 || j==3 || j==4 || j==5 || j>=10 )
    {
    y[i,j]~dnorm(mu[i,j],psi[j])I(low[z[i,j]+1],high[z[i,j]+1])
    }
    else
    {
    y[i,j]~dnorm(mu[i,j],psi[j])
    ephat[i,j]<-y[i,j]-mu[i,j]
    }

    }

    because it is not accepted by WinBUGS ! Please do you any idea how it can be written in WinBUGS ?

    Thank you for your help

    Sincerely yours,
    jade


  25. hi everyone
    i need your help to correct winbugs code is there anyone help me please.
    thanks alot


  26. I’m trying to decipher your code, so forgive me if I am missing something. I was wondering, under the section you labeled “continuity Paths”, should the last line include xi[i, 1], xi[i, 2], and xi[i, 3], instead of all three being xi[i, 1]?


  27. Hello. My name is Yoonjeong. I just started to learn Winbugs.
    Can I ask you a question?

    In multiple CFA analysis, we often constrain some parameters to be equal across groups. For example, suppose that I want to constrain the first factor loading between two groups. I named the first factor loading “lambda.g1.[1]” in first group and “lambda.g2[1]” in second group. How can I specify the equal constraint in Winbugs?


  28. hi

    i want to specify a new distribution for thresholds (inverse normal distribution),how can i do that in SEM model and how can i write this type of distribution in winbugs?

    many thanks in advance



Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: