Code mistake or Bug?

4 replies [Last post]
smedland's picture
Offline
Joined: 08/04/2009

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

Mike Cheung's picture
Offline
Joined: 10/08/2009
Hi, From mxMatrix(

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

neale's picture
Offline
Joined: 07/31/2009
Agreed, it may be easier to

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

mspiegel's picture
Offline
Joined: 07/31/2009
The mxAlgebraObjective's are

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.

smedland's picture
Offline
Joined: 08/04/2009
Thanks everyone - coding

Thanks everyone - coding after midnight is dangerous
sarah