Hello - I'm having some trouble with discrepancies in the calculation of BIC in comparing output from Mx with OpenMx on both univariate and bivariate ACE models. I am currently using OpenMx 2.0.1 and R version 3.1.2. Below is the univariate script I used. It seems from the output that the problem may be coming from the fact that mxData objects are pulling number of observations from the total number of rows in mzData and dzData, which includes some rows in which both variables are missing. Do rows in which both values are missing have to be removed ahead of time to allow number of observations to be correct (and subsequently BIC), is there another option I could use, or is this perhaps not what is causing discrepancy in the BIC values I see between Mx and OpenMx.
Also, unrelated to BIC, I've noticed in the standardization steps that some scripts use path estimates post multiplied by the iSD matrix and others have path estimates pre multiplied by iSD. I know in the bivariate case this effectively creates column standardization in post multiplication and row standardization in pre-multiplication. Therefore, the a21 element is either divided by the SD for variable 2 or the SD for variable 1, respectively. The convention is to report standardized path estimates, but is there a convention on which of the two options is more appropriate?
Thanks for your guidance,
Sarah
Vars <- 'nwidrw'
nv <- 1 # number of variables
ntv <- nv*2 # number of total variables
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")
selVars
mzData <- subset(data, zyg==1, selVars)
dzData <- subset(data, zyg==2, selVars)
pathA <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE,
values=.5, label="a11", name="a" )
pathC <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE,
values=.5, label="c11", name="c" )
pathE <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE,
values=.5, label="e11", name="e" )
covA <- mxAlgebra( expression=a %% t(a), name="A" )
covC <- mxAlgebra( expression=c %% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )
covP <- mxAlgebra( expression=(A+C+E), name="V" )
Ident <- mxMatrix ( type="Iden", nrow=nv, ncol=nv, name="I")
diagStdev <-mxAlgebra ( expression=solve(sqrt(I*V)),name="iSD")
stA <-mxAlgebra( expression=a%%iSD, name="sta")
stC <-mxAlgebra( expression=c%%iSD, name="stc")
stE <-mxAlgebra( expression=e%*%iSD, name="ste")
h2 <-mxAlgebra( expression=A/V, name="h2")
c2 <-mxAlgebra( expression=C/V, name="c2")
e2 <-mxAlgebra( expression=E/V, name="e2")
rowVars <-rep('vars',nv)
colVars <-rep(c('A','C','E','SA','SC','SE'),each=nv)
estVars <-mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="Est", dimnames=list(rowVars,colVars))
meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values= 0, label=c("meanMZ1", "meanMZ2"), name="expMeanMZ" )
meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values= 0, label=c("meanDZ1", "meanDZ2"), name="expMeanDZ" )
covMZ <- mxAlgebra( expression=
rbind( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" )
covDZ <- mxAlgebra( expression=
rbind( cbind(A+C+E, 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ",
dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ",
dimnames=selVars )
funML <- mxFitFunctionML()
pars <-list( pathA, pathC, pathE, covA, covC, covE, covP, Ident, diagStdev, estVars, stA, stC, stE, h2, c2, e2 )
modelMZ <- mxModel( pars, meanMZ, covMZ, dataMZ, expMZ, funML,
name="MZ" )
modelDZ <- mxModel( pars, meanDZ, covDZ, dataDZ, expDZ, funML,
name="DZ" )
fitML <- mxFitFunctionMultigroup(c("MZ.fitfunction", "DZ.fitfunction"))
Conf_h2 <- mxCI ('h2[1,1]' )
Conf_c2 <- mxCI ('c2[1,1]' )
Conf_e2 <- mxCI ('e2[1,1]' )
AceModel <- mxModel( "ACE", pars, modelMZ, modelDZ,
fitML, Conf_h2, Conf_c2, Conf_e2)
AceFit <- mxRun(AceModel, intervals=TRUE)
AceSumm <- summary(AceFit)
AceSumm
Hi Sarah,
When calculating BIC and other fit statistics that depend on the number of rows, OpenMx grabs the number of rows directly from the data. If you don't want those rows included in the number of rows calculation then you have to drop them before giving OpenMx the data. Here's an example.
Cheers,
Mike Hunter