Objective function returned an infinite value

2 replies [Last post]
pgseye's picture
Offline
Joined: 10/13/2009

Hi again,

I'm trying to run an ordinal model with age and sex as covariates. The ACE model seems to run Ok, but then I receive the following error with an AE model:

> univAEOrdFit <- mxRun(univAEOrdModel, intervals=T)
Running univAEOrd
Error: The job for model 'univAEOrd' exited abnormally with the error message: Objective function returned an infinite value.
In addition: Warning message:
In runHelper(model, frontendStart, intervals, silent, suppressWarnings, :
Not calculating confidence intervals because of error status.

I'm even less confident on how to set up ordinal models, so it could just be something obvious to the trained eye in the way I've amended the code.

Thanks a lot,

Paul

The (abbreviated script is):

nv <- 1 # number of variables per twin
ntv <- nv*2 # number of variables per pair
nthresh1 <- 1 # number of max thresholds

Strab[,c(32,33)] <- mxFactor(Strab[,c(32,33)], c(0 : nthresh1))
summary(Strab)
str(Strab)

Strab$Var_1<-Strab$StraightorEsoT.1
Strab$Var_2<-Strab$StraightorEsoT.2

#CREATE DEFINITION VARIABLE FOR EACH TWIN
Strab$Age_1 <- Strab$Age.1
Strab$Age_2 <- Strab$Age.2
Strab$Sex_1 <- Strab$SexCode.1
Strab$Sex_2 <- Strab$SexCode.2

defVars <- c('Age_1', 'Age_2', 'Sex_1', 'Sex_2')

Vars <- 'Var'
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="_")

mzData <- subset(Strab, Zygosity=='MZ', c(selVars,defVars))
dzData <- subset(Strab, Zygosity=='DZ', c(selVars,defVars))

# Fit ACE Model with RawData and Matrices Input, ONE overall Threshold
# ---------------------------------------------------------------------
univACEOrdModel <- mxModel("univACEOrd",
mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="a11", name="a" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="c11", name="c" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="e11", name="e" ),
# Matrices A, C, and E compute variance components
mxAlgebra( expression=a %*% t(a), name="A" ),
mxAlgebra( expression=c %*% t(c), name="C" ),
mxAlgebra( expression=e %*% t(e), name="E" ),
# Algebra to compute total variances and standard deviations (diagonal only)
mxAlgebra( expression=A+C+E, name="V" ),
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( expression=solve(sqrt(I*V)), name="sd"),
# Constraint on variance of ordinal variables
mxConstraint( V == I, name="Var1"),
# Matrix & Algebra for expected means vector
mxMatrix( type="Zero", nrow=1, ncol=nv, name="M" ),
mxAlgebra( expression= cbind(M,M), name="expMean" ),
mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.8, label='thre',name="T" ),
mxAlgebra( expression= cbind(T,T), dimnames=list('th1',selVars), name="expThre" ),

mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= 0, label=c("betaAge","betaSex"), name="beta"),
# Algebra for expected variance/covariance matrix in MZ
mxAlgebra( expression= rbind ( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" ),
# Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix
mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )
),

mxModel("MZ", mxData( observed=mzData, type="raw" ),
mxMatrix( type="Full", nrow=2, ncol=2, free=F, label=c("data.Age_1","data.Sex_1","data.Age_1","data.Sex_2"), name="MZDefVars"),
mxAlgebra( expression=ACE.expThre + ACE.beta %*% MZDefVars,dimnames=list('th1',selVars), name="expThreMZ"),
mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", thresholds="expThreMZ", dimnames=selVars ) ),

mxModel("DZ", mxData( observed=dzData, type="raw" ),
mxMatrix( type="Full", nrow=2, ncol=2, free=F, label=c("data.Age_1","data.Sex_1","data.Age_1","data.Sex_2"), name="DZDefVars"),
mxAlgebra( expression=ACE.expThre + ACE.beta %*% DZDefVars,dimnames=list('th1',selVars), name="expThreDZ"),
mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", thresholds="expThreDZ", dimnames=selVars ) ),

mxAlgebra( expression=MZ.objective + DZ.objective, name="sumll" ),
mxAlgebraObjective("sumll"),
mxCI(c('ACE.A', 'ACE.C', 'ACE.E'))
)

univACEOrdFit <- mxRun(univACEOrdModel, intervals=T)
univACEOrdSumm <- summary(univACEOrdFit)
univACEOrdSumm
LL_ACE <- mxEval(objective, univACEOrdFit)
LL_ACE

# Generate ACE Output
# -----------------------------------------------------------------------
parameterSpecifications(univACEOrdFit)
expectedMeansCovariances(univACEOrdFit)
tableFitStatistics(univACEOrdFit)

# Generate Table of Parameter Estimates using mxEval
pathEstimatesACE <- print(round(mxEval(cbind(ACE.a,ACE.c,ACE.e), univACEOrdFit),4))
varComponentsACE <- print(round(mxEval(cbind(ACE.A/ACE.V,ACE.C/ACE.V,ACE.E/ACE.V, ACE.expThre), univACEOrdFit),4))
rownames(pathEstimatesACE) <- 'pathEstimates'
colnames(pathEstimatesACE) <- c('a','c','e')
rownames(varComponentsACE) <- 'varComponents'
colnames(varComponentsACE) <- c('a^2','c^2','e^2', c('Th_Tw1','Th_Tw2'))
pathEstimatesACE
varComponentsACE

# Fit AE model
# ---------------------------------------------------------------------
univAEOrdModel <- mxModel(univACEOrdFit, name="univAEOrd",
mxModel(univACEOrdFit$ACE,
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=0, label="c11", name="c" ) # drop c at 0
)
)
univAEOrdFit <- mxRun(univAEOrdModel, intervals=T)
univAEOrdSumm <- summary(univAEOrdFit)
univAEOrdSumm
LL_AE <- mxEval(objective, univAEOrdFit)
LL_AE

neale's picture
Offline
Joined: 07/31/2009
Hi Without the Strab

Hi

Without the Strab dataframe I can't debug it thoroughly, but at a guess I would say that just fixing c to zero without adjusting any of the other parameters has left the model in a region where the variance is too small, which in turn has put one or more of the thresholds so far out that the likelihood of an observation falling beyond it is (within computational precision) equal to zero, which in turn has created an infinite value when the logarithm of the likelihood of zero is computed.

To test this hypothesis, try something like

mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.707, label="a11", name="a" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.707, label="e11", name="e" ),

in addition to the fixing c part. If that still fails, try also resetting the regressions on the covariates with:
mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= 0, label=c("betaAge","betaSex"), name="beta"),

Good luck & let us know how it goes.

pgseye's picture
Offline
Joined: 10/13/2009
Thanks Mike, That seemed to

Thanks Mike,

That seemed to fix it for that model (I didn't need to reset beta) - that's great.

Then I tried modelling the other main strabismus phenotype (exotropia - divergent turn) and this time I get the infinite value error arising in the ACE model. I tried adjusting the starting values of the path coefficients, but am still getting the problem. I'm guessing there might be something wrong with my changes to the script.

If you had a chance to look I'd really appreciate it.

Thanks,

Paul

AttachmentSize
StraightsvsExoT_AgeSex.R 14.29 KB
Strab.RData_.txt 29.68 KB