Wed, 12/01/2010 - 03:12

How does one make the variance estimates a function of the definition variable sex? I am working on a univariate script with two groups, MZ and DZ, with opposite sex twins included in the DZ group. I have tried to apply principles from the moderator script from Boulder, but I'm not sure if that is OK for the definition variable case. The script below fails due to inconformable arrays, and I am working on getting the matrices in the right format, but I need to know whether I am on the right track here. Very grateful for any help to get some steps further.

mxModel("DZ", mxData( observed=dzData, type="raw" ),

# Matrix for definition variable

mxMatrix( type="Full", nrow=2, ncol=1, free=F, label=c("data.ASEX","data.ASEX"), name="DZDefVars"),

# Algebra for making the means a function of the definition variable sex

mxAlgebra( expression=twinADE.expMean + twinADE.beta %*% DZDefVars, name="expMeanDZMod"),

# Algebra for making the variance estimates a function of the definition variable sex

mxAlgebra( expression=(twinADE.a+ DZDefVars%*%twinADE.aI) %*% t(twinADE.a+ DZDefVars%*%twinADE.aI), name="A" ),

mxAlgebra( expression=(twinADE.d+ DZDefVars%*%twinADE.dI) %*% t(twinADE.d+ DZDefVars%*%twinADE.dI), name="D" ),

mxAlgebra( expression=(twinADE.e+ DZDefVars%*%twinADE.eI) %*% t(twinADE.e+ DZDefVars%*%twinADE.eI), name="E" ),

# Algebra for expected variance/covariance matrix and expected mean vector in DZ

mxFIMLObjective( covariance="twinADE.expCovDZ", means="expMeanDZMod", dimnames=selVars )

.

Hi

the issue is that your DZDefVars variable contains the sex of both twins - which is what you want for your means correction. but for the variance correction you need to consider the sex of each twin individually.

try something like this (this is a traditional sex-limitation approach to the problem - there are more elegant ways to set up these types of models but this is more explicit)

mxMatrix( type="Full", nrow=1, ncol=1, free=F, label=c("data.ASEX"), name="sex1"),

mxMatrix( type="Full", nrow=1, ncol=1, free=F, label=c("data.BSEX"), name="sex2"),

mxAlgebra( ( (twinADE.a+ sex1%*%twinADE.aI) %*% t(twinADE.a+ sex1%*%twinADE.aI)+

(twinADE.d+ sex1%*%twinADE.dI) %*% t(twinADE.d+ sex1%*%twinADE.dI)+

(twinADE.e+ sex1%*%twinADE.eI) %*% t(twinADE.e+ sex1%*%twinADE.eI) ), name="Var1" ),

mxAlgebra( ( (twinADE.a+ sex2%*%twinADE.aI) %*% t(twinADE.a+ sex2%*%twinADE.aI)+

(twinADE.d+ sex2%*%twinADE.dI) %*% t(twinADE.d+ sex2%*%twinADE.dI)+

(twinADE.e+ sex2%*%twinADE.eI) %*% t(twinADE.e+ sex2%*%twinADE.eI) ), name="Var2" ),

mxAlgebra( ( .5%x%(twinADE.a+ sex1%*%twinADE.aI) %*% t(twinADE.a+ sex2%*%twinADE.aI)+

.25%x%(twinADE.d+ sex1%*%twinADE.dI) %*% t(twinADE.d+ sex2%*%twinADE.dI) ), name="Cov12" ),

mxAlgebra( rbind( cbind(Var1,Cov12) , cbind(Cov12,Var2) ), name="expCovDZ")

Thank you - looking forward to try this out.