Hello,
I'm building a two-level SEM in openBUGS with a mixture of ordered
categorical, dichotomous and continuous variables. Some of these
variables contain missing values at random. OpenBUGS is able to handle
missing values in the continuous variables but crashes if there are
missing values in the ordered categorical or dichotomous variables.
I used the latent variable approach described in Lee 2007 (Structural
equation models a Bayesian approach, ch 6 and 7) to handle the
categorical and dichotomous variables. Each categorical variable is
assumed to come from a underlying continuous normal variable. In order
to link the order categorical variable (z) (which contains 3 levels
k=1,2,3) to the underlying continuous one (y) I defined some threshold
values (ths) as suggested by Lee: for k=1,2,3 z equal k if
ths_(k-1)<y<=ths_(k). The dichotomous variables were treated in a
similar way, but, in this case, there are only two levels, and the
threshold is set to the value 0. To implement the transformation from
categorical to continuous variable in openBUGS, I used the censoring
function I(lower, upper).
I can run the model in openBUGS (3.2.3) through R2OpenBUGS (R version
3.3.1) only if I do not include missing data in the categorical
variables. The problem is, I believe, that the censoring function cannot
handle missing vales. I tried to specify the missingness mechanism for
the categorical variables, so I created a data frame (rho) of the same
dimensions of the categorical variable containing 0 and 1. The value of
1 indicates a missing data point in that position of the data, 0
otherwise. Then I added the missingness mechanism as follow:
for (j in 1:9){
rho[kk[g]+i,j]~dbern(ber[kk[g]+i,j])
logit(ber[kk[g]+i,j])<-b[1]+b[2]*z[kk[g]+i,1]+b[3]*z[kk[g]+i,2]+
b[4]*z[kk[g]+i,3]+b[5]*z[kk[g]+i,4]+b[6]*z[kk[g]+i,5]+
b[7]*z[kk[g]+i,6]+b[8]*z[kk[g]+i,7]+b[9]*z[kk[g]+i,8]+
b[10]*z[kk[g]+i,9]
}
9 is the number of categorical variables with missing data in the model.
rho is a data frame of 0 and 1s, z are the categorical data, the row
specification kk[g]+i is to indicate observation from individual g. It
is a two-level SEM including for the fact that some individuals in the
population have repeated observations and the errors are correlated.
This code produces an error:
"OpenBUGS version 3.2.3 rev 1012 model is syntactically correct
data loaded
variable z is not defined
model must have been compiled but not updated to be able to change RN
generator"
I also include the full model (see below). y[kk[g]+i,j] are the
underlying continuous and normally distributed variables.
I(thd[j,z[kk[g]+i,j]], thd[j,z[kk[g]+i,j]+1]) supplies threshold values
(thd) to the censoring function for the ordered categorical variables
and I(low[z[kk[g]+i,j]+1], high[z[kk[g]+i,j]+1]) supplies threshold for
the dichotomous variable.
Lee provides an example but no code to implement it. I am not sure what
I am doing wrong. Do you have any idea? Any suggestion would be really
appreciated!
Thank you in advance for your time
Best wishes, Michela
model{ for(g in 1:339){ # number of individuals for(i in 1:N[g]){ #
number of observations per individuals
#missingness mechanism
for (j in 1:9){
rho[kk[g]+i,j]~dbern(ber[kk[g]+i,j])
logit(ber[kk[g]+i,j])<-
b[1]+b[2]*z[kk[g]+i,1]+b[3]*z[kk[g]+i,2]+b[4]*z[kk[g]+i,3]+b[5]*z[kk[g]+i,4]+b[6]*z[kk[g]+i,5]+b[7]*z[kk[g]+i,6]+b[8]*z[kk[g]+i,7]+b[9]*z[kk[g]+i,8]+b[10]*z[kk[g]+i,9]
}
#measurement model
#ordered categorical variables
for (j in 1:7){
y[kk[g]+i,j]~dnorm(u[kk[g]+i,j],1)I(thd[j,z[kk[g]+i,j]],
thd[j,z[kk[g]+i,j]+1])
ephat[kk[g]+i,j]<-y[kk[g]+i,j]-u[kk[g]+i,j] # errors,
epsilon hat
}
# dichotomous variables
for (j in 8:10){
y[kk[g]+i,j]~dnorm(u[kk[g]+i,j],1)I(low[z[kk[g]+i,j]+1],
high[z[kk[g]+i,j]+1])
ephat[kk[g]+i,j]<-y[kk[g]+i,j]-u[kk[g]+i,j] # errors,
epsilon hat
}
#continuous variables
for (j in 11:20){
z[kk[g]+i,j]~dnorm(u[kk[g]+i,j],psi[j])
ephat[kk[g]+i,j]<-z[kk[g]+i,j]-u[kk[g]+i,j] # errors,
epsilon hat
}
u[kk[g]+i,1]<-
a[1,1]*fc[kk[g]+i,1]+a[1,2]*fc[kk[g]+i,2]+a[1,3]*fc[kk[g]+i,3]+pi[g,1]+eta[g,i]
u[kk[g]+i,2]<-
a[2,1]*fc[kk[g]+i,1]+a[2,2]*fc[kk[g]+i,2]+a[2,3]*fc[kk[g]+i,3]+lb[1]*pi[g,1]+lw[1]*eta[g,i]
u[kk[g]+i,3]<-
a[3,1]*fc[kk[g]+i,1]+a[3,2]*fc[kk[g]+i,3]+a[3,3]*fc[kk[g]+i,3]+lb[2]*pi[g,1]+lw[2]*eta[g,i]
u[kk[g]+i,4]<-
a[4,1]*fc[kk[g]+i,1]+a[4,2]*fc[kk[g]+i,2]+a[4,3]*fc[kk[g]+i,3]+pi[g,4]+xi[g,i,3]
u[kk[g]+i,5]<-
a[5,1]*fc[kk[g]+i,1]+a[5,2]*fc[kk[g]+i,2]+a[5,3]*fc[kk[g]+i,3]+lb[3]*pi[g,4]+lw[3]*xi[g,i,3]
u[kk[g]+i,6]<-
a[6,1]*fc[kk[g]+i,1]+a[6,2]*fc[kk[g]+i,2]+a[6,3]*fc[kk[g]+i,3]+pi[g,5]+xi[g,i,4]
u[kk[g]+i,7]<-
a[7,1]*fc[kk[g]+i,1]+a[7,2]*fc[kk[g]+i,2]+a[7,3]*fc[kk[g]+i,3]+lb[4]*pi[g,5]+lw[4]*xi[g,i,4]
# dichotomous variables
u[kk[g]+i,8]<-
a[8,1]*fc[kk[g]+i,1]+a[8,2]*fc[kk[g]+i,2]+a[8,3]*fc[kk[g]+i,3]+lb[5]*pi[g,1]+lw[5]*eta[g,i]
u[kk[g]+i,9]<-
a[9,1]*fc[kk[g]+i,1]+a[9,2]*fc[kk[g]+i,2]+a[9,3]*fc[kk[g]+i,3]+pi[g,3]+xi[g,i,2]
u[kk[g]+i,10]<-
a[10,1]*fc[kk[g]+i,1]+a[10,2]*fc[kk[g]+i,2]+a[10,3]*fc[kk[g]+i,3]+lb[6]*pi[g,3]+lw[6]*xi[g,i,2]
# continuous
u[kk[g]+i,11]<-
a[11,1]*fc[kk[g]+i,1]+a[11,2]*fc[kk[g]+i,2]+a[11,3]*fc[kk[g]+i,3]+lb[7]*pi[g,1]+lw[7]*eta[g,i]
u[kk[g]+i,12]<-
a[12,1]*fc[kk[g]+i,1]+a[12,2]*fc[kk[g]+i,2]+a[12,3]*fc[kk[g]+i,3]+pi[g,2]+xi[g,i,1]
# lambda2 and lambda1set=1
u[kk[g]+i,13]<-
a[13,1]*fc[kk[g]+i,1]+a[13,2]*fc[kk[g]+i,2]+a[13,3]*fc[kk[g]+i,3]+lb[8]*pi[g,2]+lw[8]*xi[g,i,1]
u[kk[g]+i,14]<-
a[14,1]*fc[kk[g]+i,1]+a[14,2]*fc[kk[g]+i,2]+a[14,3]*fc[kk[g]+i,3]+lb[9]*pi[g,2]+lw[9]*xi[g,i,1]
u[kk[g]+i,15]<-
a[15,1]*fc[kk[g]+i,1]+a[15,2]*fc[kk[g]+i,2]+a[15,3]*fc[kk[g]+i,3]+lb[10]*pi[g,3]+lw[10]*xi[g,i,2]
u[kk[g]+i,16]<-
a[16,1]*fc[kk[g]+i,1]+a[16,2]*fc[kk[g]+i,2]+a[16,3]*fc[kk[g]+i,3]+lb[11]*pi[g,3]+lw[11]*xi[g,i,2]
u[kk[g]+i,17]<-
a[17,1]*fc[kk[g]+i,1]+a[17,2]*fc[kk[g]+i,2]+a[17,3]*fc[kk[g]+i,3]+lb[12]*pi[g,4]+lw[12]*xi[g,i,3]
u[kk[g]+i,18]<-
a[18,1]*fc[kk[g]+i,1]+a[18,2]*fc[kk[g]+i,2]+a[18,3]*fc[kk[g]+i,3]+lb[13]*pi[g,4]+lw[13]*xi[g,i,3]
u[kk[g]+i,19]<-
a[19,1]*fc[kk[g]+i,1]+a[19,2]*fc[kk[g]+i,2]+a[19,3]*fc[kk[g]+i,3]+lb[14]*pi[g,5]+lw[14]*xi[g,i,4]
u[kk[g]+i,20]<-
a[20,1]*fc[kk[g]+i,1]+a[20,2]*fc[kk[g]+i,2]+a[20,3]*fc[kk[g]+i,3]+lb[15]*pi[g,5]+lw[15]*xi[g,i,4]
xi[g,i,1:4]~dmnorm(ux[1:4],phi[1:4,1:4]) # ux=[0 0]^T is
fixed constant
eta[g,i]~dnorm(nu[g,i], psd) #psd e' semplicemente un par from
dgamma
nu[g,i]<- gam[1]*xi[g,i,1]+gam[2]*xi[g,i,2]+gam[3]*xi[g,i,3] +
gam[4]*xi[g,i,4]
dthat[g,i]<-eta[g,i]-nu[g,i] # questi sono gli errori delta
}# end of i
pi[g,1:5]~ dmnorm(uu[1:5],phip[1:5,1:5]) #distributions of omega2
set mean to zero not sure why
}# end of g
uu[1]<- 0.0 uu[2]<- 0.0 uu[3]<- 0.0 uu[4]<- 0.0 uu[5]<- 0.0
ux[1]<- 0.0 ux[2]<- 0.0 ux[3]<-0.0 ux[4]<-0.0
# priors on loadings and coefficients
a[1,1]~dnorm(0.0,4.0) a[2,1]~dnorm(-1.0,4.0)
a[3,1]~dnorm(1.0,4.0)
a[4,1]~dnorm(0.0,4.0) a[5,1]~dnorm(0.0,4.0)
a[6,1]~dnorm(-1.0,4.0)
a[7,1]~dnorm(-1.0,4.0) a[8,1]~dnorm(0.0,4.0) a[9,1]~dnorm(0.0,4.0)
a[10,1]~dnorm(-0.3,4.0) a[11,1]~dnorm(0.0,4.0)
a[12,1]~dnorm(1.0,4.0)
a[13,1]~dnorm(-0.5,4.0) a[14,1]~dnorm(0.0,4.0)
a[15,1]~dnorm(0.0,4.0)
a[16,1]~dnorm(0.0,4.0) a[17,1]~dnorm(0.1,4.0)
a[18,1]~dnorm(1.0,4.0)
a[19,1]~dnorm(0.5,4.0) a[20,1]~dnorm(0.0,4.0)
a[2,2]~dnorm(1.0,1) a[3,3]~dnorm(1.0,1) a[8,2]~dnorm(1.0,1)
a[1,2]<- 0.0 a[3,2]<- 0.0 a[4,2]<- 0.0 a[5,2]<- 0.0 a[6,2]<-
0.0 a[7,2]<- 0.0
a[9,2]<- 0.0 a[10,2]<- 0.0 a[11,2]<- 0.0 a[12,2]<- 0.0 a[13,2]<-
0.0 a[14,2]<- 0.0
a[15,2]<- 0.0 a[16,2]<- 0.0 a[17,2]<- 0.0 a[18,2]<- 0.0 a[19,2]<-
0.0 a[20,2]<- 0.0
a[1,3]<- 0.0 a[2,3]<- 0.0 a[4,3]<- 0.0 a[5,3]<- 0.0 a[6,3]<-
0.0 a[7,3]<- 0.0
a[8,3]<- 0.0 a[9,3]<- 0.0 a[10,3]<- 0.0 a[11,3]<- 0.0 a[12,3]<-
0.0 a[13,3]<- 0.0
a[14,3]<- 0.0 a[15,3]<- 0.0 a[16,3]<- 0.0 a[17,3]<- 0.0 a[18,3]<-
0.0 a[19,3]<- 0.0
a[20,3]<- 0.0
var.bw[1]<-4.0*psi[2] var.bw[2]<-4.0*psi[3] var.bw[3]<-4.0*psi[5]
var.bw[4]<-4.0*psi[7] var.bw[5]<-4.0*psi[8] var.bw[6]<-4.0*psi[10]
var.bw[7]<-4.0*psi[11] var.bw[8]<-4.0*psi[13] var.bw[9]<-4.0*psi[14]
var.bw[10]<-4.0*psi[15] var.bw[11]<-4.0*psi[16] var.bw[12]<-4.0*psi[17]
var.bw[13]<-4.0*psi[18] var.bw[14]<-4.0*psi[19] var.bw[15]<-4.0*psi[20]
lb[1]~dnorm(1.0,var.bw[1]) lb[2]~dnorm(1.0,var.bw[2])
lb[3]~dnorm(1.0,var.bw[3])
lb[4]~dnorm(1.0,var.bw[4]) lb[5]~dnorm(1.0,var.bw[5])
lb[6]~dnorm(1.0,var.bw[6])
lb[7]~dnorm(1.0,var.bw[7]) lb[8]~dnorm(1.0,var.bw[8])
lb[9]~dnorm(1.0,var.bw[9])
lb[10]~dnorm(1.0,var.bw[10]) lb[11]~dnorm(1.0,var.bw[11])
lb[12]~dnorm(1.0,var.bw[12])
lb[13]~dnorm(1.0,var.bw[13]) lb[14]~dnorm(1.0,var.bw[14])
lb[15]~dnorm(1.0,var.bw[15])
lw[1]~dnorm(0.0,var.bw[1]) lw[2]~dnorm(0.0,var.bw[2])
lw[3]~dnorm(0.8,var.bw[3])
lw[4]~dnorm(2.6,var.bw[4]) lw[5]~dnorm(3.0,var.bw[5])
lw[6]~dnorm(1.0,var.bw[6])
lw[7]~dnorm(0.5,var.bw[7]) lw[8]~dnorm(-0.2,var.bw[8])
lw[9]~dnorm(0.0,var.bw[9])
lw[10]~dnorm(1.0,var.bw[10]) lw[11]~dnorm(0.0,var.bw[11])
lw[12]~dnorm(0.0,var.bw[12])
lw[13]~dnorm(0.9,var.bw[13]) lw[14]~dnorm(0.3,var.bw[14])
lw[15]~dnorm(0.0,var.bw[15])
var.gam<-4.0*psd
gam[1]~dnorm(0.1,var.gam) gam[2]~dnorm(-0.1,var.gam)
gam[3]~dnorm(0.0,var.gam) gam[4]~dnorm(0.0,var.gam)
# priors on precisions
for(j in 1:20){psi[j]~dgamma(10.0,4.0)
ivpsi[j]<-1/psi[j]}
psd~dgamma(10.0,4.0)
ivpsd<-1/psd
phi[1:4,1:4]~dwish(R0[1:4,1:4],5)
phx[1:4,1:4]<-inverse(phi[1:4,1:4])
phip[1:5,1:5]~dwish(R1[1:5,1:5],5)
php[1:5,1:5]<-inverse(phip[1:5,1:5])
#priors missingness mechanism
b[1]~dnorm(0,0.1) b[2]~dnorm(0,0.1) b[3]~dnorm(0,0.1)
b[4]~dnorm(0,0.1) b[5]~dnorm(0,0.1) b[6]~dnorm(0,0.1)
b[7]~dnorm(0,0.1) b[8]~dnorm(0,0.1) b[9]~dnorm(0,0.1)
}# end of model
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|