You are here

Code mistake or Bug?

5 posts / 0 new
Last post
smedland's picture
Offline
Joined: 08/04/2009 - 16:08
Code mistake or Bug?

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 - 22:37
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, C11C31 ]
[ C11
C21, C22^2 +C21^2, C21C31 ]
[ C11
C31, C21C31, 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 - 15:14
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 - 15:24
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 - 16:08
Thanks everyone - coding

Thanks everyone - coding after midnight is dangerous
sarah