Sat, 08/14/2010 - 22:57

Hi

I am teaching tomorrow and I have a werid glich in my code that I can't resolve

I am modeling the correlations between 3 variables in a very straight forward script. I then want to drop one of the correlations. the difference in fits between the full model and the reduced model is huge and the correlation is not being set to zero

Where have I gone wrong - I can't spot it typo(?)

thanks

Sarah

data(twinData)

selVars <- c('wt1','ht1','bmi1')

nv <- 3

factorData <- subset(twinData, zyg<=5, selVars)

SatModel <- mxModel("saturated",

mxModel("Cor",

mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.5, name="Chol" ),

mxAlgebra( Chol %*% t(Chol), name="expCov" ),

mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=c(60,1.5,20), name="expMean" ),

mxData( observed=factorData, type="raw" ),

mxFIMLObjective( covariance="expCov", means="expMean", dimnames=selVars)),

mxAlgebra( Cor.objective , name="modelfit" ),

mxAlgebraObjective("modelfit")

)

SatFit <- mxRun(SatModel)

SatSumm <- summary(SatFit)

SatSumm

Drop1Model <- mxModel("drop1",

mxModel("Cor",

mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c( TRUE,TRUE,TRUE,TRUE,FALSE,TRUE), name="Chol" ),

mxAlgebra( Chol %*% t(Chol), name="expCov" ),

mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=c(60,1.5,20), name="expMean" ),

mxData( observed=factorData, type="raw" ),

mxFIMLObjective( covariance="expCov", means="expMean", dimnames=selVars)),

mxAlgebra( Cor.objective , name="modelfit" ),

mxAlgebraObjective("modelfit")

)

Drop1Fit <- mxRun(Drop1Model)

Drop1Summ <- summary(Drop1Fit)

Drop1Summ

Drop1Fit$Cor$Chol

Drop1Fit$Cor$expCov

Hi,

From mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c( TRUE,TRUE,TRUE,TRUE,FALSE,TRUE), name="Chol" ), Chol is defined as

[ C11, 0, 0 ]

[ C21, C22, 0 ]

[ C31, 0 , C33 ].

Then the "expCov", which is Chol %*% t(Chol), is

[ C11^2, C11*C21, C11*C31 ]

[ C11*C21, C22^2 +C21^2, C21*C31 ]

[ C11*C31, C21*C31, C33^2 +C31^2 ].

Thus, the covariance and correlation between 'ht1' and 'bmi1' is not zero but C21*C31.

Would it be easier to model it with just a covariance matrix? Than you can directly fix the covariance between 'ht1' and 'bmi1' at zero.

Regards,

Mike

Agreed, it may be easier to specify the model as a covariance matrix directly, but sometimes this is hazardous as the covariance matrix can easily become non-positive definite.

Two comments:

i) simply setting a matrix element to free=FALSE does not set it to zero, simply to whatever was in there originally. In the present case, a complete re-specification of the reduced model does not include starting values for any of the parameters, which may cause some problems in the search for the solution. I expect that this is the cause of the really bad fit. Can be easier to take first solution as input for the submodel, and tweak it with a mxMatrix command only.

ii) Agree with Mike Cheung that the covariance between 2&3 won't be zero. You could use a non-linear constraint, something like mxConstraint(expCov[3,2]==0,name="onlyAllowSolutionsWhereCov23isZero")

The mxAlgebraObjective's are redundant in the model specifications. In the two models, the mxAlgebraObjective() each reference a single algebra, each of which references a single objective function. You may wish to eliminate the outer model, and just execute the inner model that contains the mxFIMLObjective() function.

Thanks everyone - coding after midnight is dangerous

sarah