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
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 ntv nthresh1

Strab[,c(32,33)]
summary(Strab)

str(Strab)

Strab$Var_1 Strab$Var_2

#CREATE DEFINITION VARIABLE FOR EACH TWIN

Strab$Age_1
Strab$Age_2
Strab$Sex_1
Strab$Sex_2

defVars

Vars selVars

mzData dzData

# Fit ACE Model with RawData and Matrices Input, ONE overall Threshold

# ---------------------------------------------------------------------

univACEOrdModel
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
univACEOrdSumm
univACEOrdSumm

LL_ACE
LL_ACE

# Generate ACE Output

# -----------------------------------------------------------------------

parameterSpecifications(univACEOrdFit)

expectedMeansCovariances(univACEOrdFit)

tableFitStatistics(univACEOrdFit)

# Generate Table of Parameter Estimates using mxEval

pathEstimatesACE
varComponentsACE
rownames(pathEstimatesACE)
colnames(pathEstimatesACE)
rownames(varComponentsACE)
colnames(varComponentsACE)
pathEstimatesACE

varComponentsACE

# Fit AE model

# ---------------------------------------------------------------------

univAEOrdModel
mxModel(univACEOrdFit$ACE,

mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=0, label="c11", name="c" ) # drop c at 0

)

)

univAEOrdFit
univAEOrdSumm
univAEOrdSumm

LL_AE
LL_AE

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.

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