Bayes CFA

[Lazar]

Hi Folks,

I have several questions about Bayes CFA in Open bugs that I hope you might be able to help me out with.

First the setup:
Here is a basic two factor CFA set up for the HS data that is available in several R packages

###Bugs CFA model 2 factors. All cross loadings zero###
#Loadings constrained to be positive
#Mildly informative priors used  on litem intercepts
#n = number of observations
#t = equals number of indicator items

model{
  for (i in 1:N){
      #latent 1
      for (t in 1:3){
        y[i,t] ~ dnorm(condmn[i,t], invsig2[t])
        condmn[i,t] <-mu[t] + fload[t]*fscore[i,1]
      }
      #latent 2
      for (t in 4:6){
        y[i,t] ~dnorm(condmn[i,t], invsig2[t])
        condmn[i,t] <-mu[t] + fload[t]*fscore[i,2]
      }
  fscore[i,1:2]~dmnorm(mn.fs[], sig.fs[,])
  }
  mn.fs[1]<-0
  mn.fs[2]<-0
  sig.fs[1,1]<-1
  sig.fs[2,2]<-1
  sig.fs[1,2]<-phi
  sig.fs[2,1]<-phi
  #Prior distribution
  phi ~ dunif(-1,1)
  for (t in 1:6){
    fload[t] ~ dnorm(0,1.0E-3)I(0,)
    invsig2[t] <- 1/psi[t]
    psi[t] ~ dunif(0,400)
    mu[t] ~ dnorm(20,.05)
  }
}

and here is the R script I am using to run it:

path <- 'C:/Users/30016475/Dropbox/Projects_Research/PVsimulation'
setwd(path)

#Load relavent libraries
library(simsem);library(R2OpenBUGS); library(MBESS); library(lavaan)
#Load HS data
data(HS.data)

#Run two factor model in Lavaan as a bench mark to check Open BUGS results 
#CFA using ML
Model <- '
  L1 =~ visual + cubes + flags
  L2 =~ paragrap + sentence + wordm
'
fit<- sem(Model, data=HS.data, std.lv=TRUE)
summary(fit)


#Open bugs run of same two factor model
#Load items into a matrix 
y<-as.matrix(HS.data[,c('visual', 'cubes', 'flags', 'paragrap', 'sentence', 'wordm')])
#Number of cases
N<-nrow(HS.data)
#No inits  (I will rely on defaults) because I am lazy....maybe this is my problem

#Data and parameters to monitor
data<-list('N', 'y')
params<-c('fscore', 'sig.fs', 'fload', 'mu')#Note fscore are the plausible values. Can delete to reduce size of output
#Bugs call
out <- bugs(data, parameters.to.save=params, inits=NULL,
            model.file='C:/Users/30016475/Dropbox/Projects_Research/PVsimulation/CFAsimple.txt',
            debug=TRUE, n.iter=1000)
#check fit
all(out$summary[,"Rhat"]<1.1)
#Check results
out$summary

Ok I am going to start with an easy question first.
The results I get for the bugs run and lavaan are close for the most part but further away than I would have guessed given I am mostly using uninformative priors. The real problem is that the correlation is correct is size but in the wrong direction. Like I have accidentally multiplied it by -1. Can anyone see the mistake I am making?

 

[Lazar]

Ok author of the paper I was using as a template said there was an error in the code. His suggested fix was to use the inverse function BUT this function exists in JAGS but not openbugs. Does anyone know the open bugs equivalent?

 

[Lazar]

Corrected bug script is:

###Bugs CFA model 2 factors. All cross loadings zero###
#Loadings constrained to be positive
#Mildly informative priors used  on litem intercepts
#n = number of observations
#t = equals number of indicator items

model{
  for (i in 1:N){
      #latent 1
      for (t in 1:3){
        y[i,t] ~ dnorm(condmn[i,t], invsig2[t])
        condmn[i,t] <-mu[t] + fload[t]*fscore[i,1]
      }
      #latent 2
      for (t in 4:6){
        y[i,t] ~dnorm(condmn[i,t], invsig2[t])
        condmn[i,t] <-mu[t] + fload[t]*fscore[i,2]
      }
  fscore[i,1:2]~dmnorm(mn.fs[], siginv.fs[,])
  }
  mn.fs[1]<-0
  mn.fs[2]<-0
  siginv.fs <- inverse(sig.fs)
  sig.fs[1,1]<-1
  sig.fs[2,2]<-1
  sig.fs[1,2]<-phi
  sig.fs[2,1]<-phi
  #Prior distribution
  phi ~ dunif(-1,1)
  for (t in 1:6){
    fload[t] ~ dnorm(0,1.0E-3)I(0,)
    invsig2[t] <- 1/psi[t]
    psi[t] ~ dunif(0,400)
    mu[t] ~ dnorm(20,.05)
  }
}

 

[Dason]

I don't have bugs installed right now but if the issue is the inverse you can easily just do the inverse directly.

siginv.fs[1,1] <- 1/(1-phi*phi)
siginv.fs[2,2] <- 1/(1-phi*phi)
siginv.fs[1,2] <- (-phi)/(1-phi*phi)
siginv.fs[2,1] <- (-phi)/(1-phi*phi)

but the manual implies that it might work if you just use the following syntax

siginv.fs <- inverse(sig.fs[,])

 

[Lazar]

Thanks Dason. I tried inverse(sig.fs[,]) but I get:

empty slot not allowed in variable name error pos 617

in openbugs.

Not such trouble in jags.

 

[Lazar]

siginv.fs[1,1] <- 1/(1-phi*phi)
siginv.fs[2,2] <- 1/(1-phi*phi)
siginv.fs[1,2] <- (-phi)/(1-phi*phi)
siginv.fs[2,1] <- (-phi)/(1-phi*phi)

Works fine