You are here

Calculation of BIC discrepancy?

2 posts / 0 new
Last post
slukow's picture
Offline
Joined: 02/25/2015 - 14:31
Calculation of BIC discrepancy?

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

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
Drop the rows

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.

mzData #data set, has some rows of missing values
selVars # character vector of names of variables you're analyzing
notAllMissing <-apply(mzData[,selVars], 1, function(x){!all(is.na(x))})
# logical vector with one element per row of mzData
# TRUE if not all of selVars are missing for that row
 
mzData <- mzData[notAllMissing,]
# Keep only the rows of mzData that are not all missing for selVars

Cheers,
Mike Hunter