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(IV)), 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