Hi! guys,

I am doing the threshold model for ordinal variables. and I finished the variables for 3 levels. However, I need modify the script to another variable which contain more 5 levels. I sucess with saturated model wiithout increments It gives me 18 parameters : rMZ,rDZ,4 threshold for MZ1 4 threshold for MZ2,4 threshold for DZ1 and 4 threshold for DZ2,respectively. But I was stuck at ADE model. Some problems in increments.

My script like this:

Set Starting Values

thVals <-c(-0.226,-0.099,0.637,2.628) # start value for thresholds t1,t2,t3,t4.
corValsM <-.876 # start value for MZ correlations
corValsD <-.316 # start value for DZ correlations

SATURATED MODEL:

Matrices for expected Means & Thresholds (on liabilities) in MZ & DZ twin pairs

MeanG <-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMean" )
meanG <-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMean" )
threMZ <-mxMatrix(type="Full", nrow=nth, ncol=ntv, free=TRUE, values=thVals, labels=c("t1MZ1","t1MZ2","t2MZ1","t2MZ2","t3MZ1","t3MZ2","t4MZ1","t4MZ2"), name="expThreMZ" )
threDZ <-mxMatrix(type="Full", nrow=nth, ncol=ntv, free=TRUE, values=thVals, labels=c("t1DZ1","t1DZ2","t2DZ1","t2DZ2","t3DZ1","t3DZ2","t4DZ1","t4DZ2"), name="expThreDZ" )
corMZ <-mxMatrix(type="Stand", nrow=ntv, ncol=ntv, free=T, values=corValsM, lbound=-.99, ubound=.99, labels=c("rMZ"), name="expCorMZ")
corDZ <-mxMatrix(type="Stand", nrow=ntv, ncol=ntv, free=T, values=corValsD, lbound=-.99, ubound=.99, labels=c("rDZ"), name="expCorDZ")

Data objects for Multiple Groups

dataMZ <-mxData(mzDataBin, type="raw")
dataDZ <-mxData(dzDataBin, type="raw")

Objective objects for Multiple Groups

objMZ <-mxFIMLObjective( covariance="expCorMZ", means="expMean", dimnames=selVars, thresholds="expThreMZ" )
objDZ <-mxFIMLObjective( covariance="expCorDZ", means="expMean", dimnames=selVars, thresholds="expThreDZ" )

Combine Groups

groupMZ <-mxModel("MZ", corMZ, meanG, threMZ, dataMZ, objMZ )
groupDZ <-mxModel("DZ", corDZ, meanG, threDZ, dataDZ, objDZ )
minus2ll<-mxAlgebra( MZ.objective + DZ.objective, name="minus2sumloglikelihood" )
obj <-mxAlgebraObjective("minus2sumloglikelihood")
ciCor <-mxCI(c('MZ.expCorMZ[2,1]', 'DZ.expCorDZ[2,1]'))
ciThre <-mxCI( c('MZ.expThreMZ','DZ.expThreDZ' ))
twinSatModel <-mxModel( "twinSat", minus2ll, obj, groupMZ, groupDZ, ciCor, ciThre )

-----------------------------------------------------------------------

RUN SATURATED MODEL (Tetrachoric correlations)

-----------------------------------------------------------------------

twinSatFit <- mxRun(twinSatModel, intervals=F)
twinSatSumm <- summary(twinSatFit)
round(twinSatFit@output$estimate,4)
twinSatSumm

It works.

FOR ACEmodel:

athA <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.76, label="a11", name="a" )
pathC <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.23, label="c11", name="c" )
pathE <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, label="e11", name="e" )

Matrices generated to hold A, C, and E computed Variance Components

covA <- mxAlgebra( expression=a %% t(a), name="A" )
covC <- mxAlgebra( expression=c %% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )

Matrix & Algebra for expected means vector and expected thresholds

mean <- mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMean" )
Inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="Inc" )

thresh <- mxMatrix( type="Full", nrow=nth, ncol=1, free=TRUE,
values=c(-0.226,0.127,0.736,1.892),
lbound=c(-3,.001),
labels=c("t1","inc21","inc32","inc43"),
name="Thre")

threshold <- mxAlgebra( expression= cbind(Inc %% Thre,Inc %% Thre), name="expThre" )

Algebra to compute total variances and standard deviations (diagonal only)

covP <- mxAlgebra( expression=A+C+E, name="V" )

Algebras generated to hold Parameter Estimates and Derived Variance Components

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="Vars", dimnames=list(rowVars,colVars))

Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins

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" )

Constraint on variance of Binary variables

matUnv <- mxMatrix( type="Unit", nrow=nv, ncol=1, name="Unv1" )
var1 <- mxConstraint( expression=diag2vec(V)==Unv1, name="Var1" )

Data objects for Multiple Groups

dataMZ <- mxData( observed=mzDataBin, type="raw" )
dataDZ <- mxData( observed=dzDataBin, type="raw" )

Objective objects for Multiple Groups

objMZ <- mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="expThre" )
objDZ <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="expThre" )

Combine Groups

pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, Inc, matUnv, estVars )
modelMZ <- mxModel( pars, mean, thresh, threshold, covMZ, dataMZ, objMZ, name="MZ" )
modelDZ <- mxModel( pars, mean, thresh, threshold, covDZ, dataDZ, objDZ, name="DZ" )
minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj <- mxAlgebraObjective( "m2LL" )
AceModel <- mxModel( "ACE", pars, var1, modelMZ, modelDZ, minus2ll, obj )

tableFitStatistics(twinSatFit,aceFitConst)

it only gives me 5 parameters in ACE model: a11, c11,e11, t1, inc.

I think there must be some thing wrong with Increment.

Does anyone have some good idea about how to fix it? Thank you so much for that!