Attachment | Size |
---|---|
Example3.R | 3.69 KB |
I'm running a measurement model (two time point common factor model) and I am getting an error message I've never seen before when I drop a parameter. I talked to AdminRobK about this, and he suggested I post it to the forums. The error reads:
Error in if (class(fit) == "try-error" || is.na(fit$output$minimum) || : missing value where TRUE/FALSE needed
In addition: Warning message:
In is.na(fit$output$minimum) : is.na() applied to non-(list or vector) of type 'NULL'
This only happens when I include the 'intervals=T' option. If this option is removed, it runs fine. Script is below (and attached).
My system information:
OpenMx version: 2.3.1 [GIT v2.3.1]
R version: R version 3.2.3 (2015-12-10)
Platform: x86_64-apple-darwin13.4.0
MacOS: 10.11.2
Default optimiser: NPSOL
require(OpenMx)
mzData <- cbind(1+rnorm(100, mean=.5, sd=.5),1+rnorm(100, mean=.5, sd=.5),1+rnorm(100, mean=.5, sd=.5),1+rnorm(100, mean=.5, sd=.5))
colnames(mzData) <- c("t1_v1", "t1_v2", "t2_v1", "t2_v2")
dzData <- cbind(1+rnorm(100, mean=.5, sd=.5),1+rnorm(100, mean=.5, sd=.5),1+rnorm(100, mean=.5, sd=.5),1+rnorm(100, mean=.5, sd=.5))
colnames(dzData) <- c("t1_v1", "t1_v2", "t2_v1", "t2_v2")
nv <- 2 # number of variables
ntv <- nv*2 # number of total variables
selVars <- c("t1_v1", "t1_v2", "t2_v1", "t2_v2")
nf <- 1
laAs <- paste("as",1:nv,1:nv,sep="")
laCs <- paste("cs",1:nv,1:nv,sep="")
laEs <- c("es_1_1", "es_2_2")
laFl <- "load"
svMe <- c(.5,.5,.5,.5)
laMe <- paste(selVars,"Mean",sep="")
laMeMZ <- paste(selVars,"MeanMZ",sep="")
laMeDZ <- paste(selVars,"MeanDZ",sep="_")
pathAl <- mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=.5, labels="al_1_1", lbound=.00001, name="al" )
pathCl <- mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=.5, labels="cl_1_1", lbound=.00001, name="cl" )
pathEl <- mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=.5, labels="el_1_1", lbound=.00001, name="el" )
covarLP <- mxAlgebra( expression= al %% t(al) + cl %% t(cl) + el %*% t(el), name="CovarLP" )
varLP <- mxAlgebra( expression= diag2vec(CovarLP), name="VarLP" )
unit <- mxMatrix( type="Unit", nrow=nf, ncol=1, name="Unit")
varLP1 <- mxConstraint( expression=VarLP == Unit, name="varLP1")
pathF <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.5, labels=laFl, name="fl", lbound=.00001 )
pathEs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.5, labels=laEs, lbound=.00001, name="es" )
covA <- mxAlgebra( expression=fl %&% (al %% t(al)), name="A" )
covC <- mxAlgebra( expression=fl %&% (cl %% t(cl)), name="C" )
covE <- mxAlgebra( expression=fl %&% (el %% t(el)) + es %% t(es), name="E" )
covP <- mxAlgebra( expression=A+C+E, name="V" )
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, label=c("meanv1", "meanv2", "meanv1", "meanv2"), name="expMean" )
covMZ <- mxAlgebra( expression= rbind( cbind(V , A+C), cbind(A+C , V)), name="expCovMZ" )
covDZ <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+C), cbind(0.5%x%A+C ,V)), name="expCovDZ" )
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars )
funML <- mxFitFunctionML()
pars <- list( pathAl, pathCl, pathEl, covarLP, varLP, unit, pathF, pathEs, covA, covC, covE, covP, matI, invSD)
modelMZ <- mxModel( "MZ", pars, covMZ, meanG, dataMZ, expMZ, funML )
modelDZ <- mxModel( "DZ", pars, covDZ, meanG, dataDZ, expDZ, funML )
minus2ll <- mxAlgebra( MZ.objective+ DZ.objective, name="minus2sumloglikelihood" )
obj <- mxFitFunctionAlgebra( "minus2sumloglikelihood" )
ComAceModel <- mxModel( "ComACE", pars, varLP1, modelMZ, modelDZ, minus2ll, obj, mxCI(c('al', 'cl', 'el', 'fl', 'es')))
ComAceFit <- mxRun(ComAceModel, intervals=T)
ComAceFit <- mxTryHard(ComAceModel, intervals=T)
ComAceSum <- summary(ComAceFit)
ComAceSum
ComAeFit <- mxModel(ComAceFit, name="ComAeFit" )
ComAeFit <- omxSetParameters(ComAeFit, labels="cl_1_1", values=0, free=F)
ComAeFit <- mxRun(ComAeFit, intervals=T)
ComAeFit <- mxTryHard(ComAeFit, intervals=T)
ComAeSum <- summary(ComAeFit)
ComAeSum
I modified the end parts of your script as follows, to diagnose the problem:
This revealed the problem.
mxTryHard()
handles requests for confidence intervals by first making attempts to find an acceptable MLE. Then, from that solution, it gives the model a custom compute plan that goes straight to estimating confidence limits. You can see a compute plan like this for yourself by queryingComAceFit$compute
. This compute plan is not a member of class 'MxComputeDefault'. In order to be compatible with custom compute plans,mxTryHard()
does not mess with them.mxRun()
, which is called bymxTryHard()
, checks to see if the model being run has a custom compute plan. If it does, and that MxCompute object has a value ofFALSE
in its '.persist' slot, it clears that object and creates a default compute plan for the model. If, on the other hand, the '.persist' slot containsTRUE
, then it leaves the MxCompute object alone.The essential problem is that the custom compute plan for finding confidence intervals that
mxTryHard()
creates should have its '.persist' slot set toFALSE
, but it doesn't.Edit: On second thought, the compute plan created by
mxTryHard()
can't have '.persist' set toFALSE
, because that would causemxRun()
to overwrite it with a default compute plan.Double edit: The compute plan should probably have .'persist' set to
FALSE
when it is returned frommxTryHard()
as part of the resulting MxModel.Anyhow, your AE model never does its primary optimization to find the MLE. The error gets thrown because, apparently, the 'minimum' element of a model's output is never created if there is no primary optimization.
mxTryHard()
expects to be able to find this element but it can't, leading to a fatal error.This should be an easy bug to repair. In the meantime, you can work around the problem by running
ComAceFit@compute <- NULL
before creating your AE model.