# Bayes CFA

#### Lazar

##### Phineas Packard
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
Code:
###Bugs CFA model 2 factors. All cross loadings zero###
#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])
}
#latent 2
for (t in 4:6){
y[i,t] ~dnorm(condmn[i,t], invsig2[t])
}
fscore[i,1:2]~dmnorm(mn.fs[], sig.fs[,])
}
mn.fs<-0
mn.fs<-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){
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:
Code:
path <- 'C:/Users/30016475/Dropbox/Projects_Research/PVsimulation'
setwd(path)

library(simsem);library(R2OpenBUGS); library(MBESS); library(lavaan)
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
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
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

##### Phineas Packard
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

##### Phineas Packard
Corrected bug script is:
Code:
###Bugs CFA model 2 factors. All cross loadings zero###
#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])
}
#latent 2
for (t in 4:6){
y[i,t] ~dnorm(condmn[i,t], invsig2[t])
}
fscore[i,1:2]~dmnorm(mn.fs[], siginv.fs[,])
}
mn.fs<-0
mn.fs<-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){
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.

Code:
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

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

Last edited:

#### Lazar

##### Phineas Packard
Thanks Dason. I tried inverse(sig.fs[,]) but I get:
Code:
empty slot not allowed in variable name error pos 617
in openbugs.

Not such trouble in jags.

#### Lazar

##### Phineas Packard
Code:
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 