Hello,
I'm quite new to twin modelling and have come across the following problem.
Before I run my proper model, I'm trying to account for any existing sex differences on the phonetypes of interest.
Am running the following model for each trait:
Univariate sex-limitation ####RG is FREE
#
unSexACEModel <- mxModel("unSexACE",
mxModel("ACE",
# ace path coefficients for males and females
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="am11", name="am"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="cm11", name="cm"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.5, label="em11", name="em"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="af11", name="af"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="cf11", name="cf"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.5, label="ef11", name="ef"),
# ACE variance components for males and females
mxAlgebra(am %*% t(am), name="Am"),
mxAlgebra(cm %*% t(cm), name="Cm"),
mxAlgebra(em %*% t(em), name="Em"),
mxAlgebra(af %*% t(af), name="Af"),
mxAlgebra(cf %*% t(cf), name="Cf"),
mxAlgebra(ef %*% t(ef), name="Ef"),
# Total variation
mxAlgebra(Am+Cm+Em, name="Vm"),
mxAlgebra(Af+Cf+Ef, name="Vf"),
# Estimate genetic (or environmental) correlation in opposite-sex pairs
mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=0.5, label="Rg", lbound=0, ubound=0.5, name="rg"),
mxMatrix("Full", nrow=1, ncol=1, free=FALSE, values=1, label="Rc", lbound=0, ubound=1, name="rc"),
# Standard deviations for males and females
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( expression=solve(sqrt(I*Vm)), name="isdm"),
mxAlgebra( expression=solve(sqrt(I*Vf)), name="isdf"),
# Standardized path estimates
mxAlgebra(am %*% isdm, name="stam"),
mxAlgebra(cm %*% isdm, name="stcm"),
mxAlgebra(em %*% isdm, name="stem"),
mxAlgebra(af %*% isdf, name="staf"),
mxAlgebra(cf %*% isdf, name="stcf"),
mxAlgebra(ef %*% isdf, name="stef"),
# Standardized variance components
mxAlgebra(Am/Vm, name="h2m"),
mxAlgebra(Cm/Vm, name="c2m"),
mxAlgebra(Em/Vm, name="e2m"),
mxAlgebra(Af/Vf, name="h2f"),
mxAlgebra(Cf/Vf, name="c2f"),
mxAlgebra(Ef/Vf, name="e2f"),
# Expected means vectors: male-male pairs; fem-fem pairs; male-fem opposite-sex pairs
mxMatrix("Full", nrow=1, ncol=nv, free=TRUE, values= .6, label="meanM", name="Mm"),
mxMatrix("Full", nrow=1, ncol=nv, free=TRUE, values= .6, label="meanF", name="Mf"),
mxAlgebra(cbind(Mm,Mm), name="expMeanM"),
mxAlgebra(cbind(Mf,Mf), name="expMeanF"),
mxAlgebra(cbind(Mm,Mf), name="expMeanO"),
# Expected covariance matrix in MZ, for males and females
mxAlgebra(
rbind(
cbind(Am+Cm+Em, Am+Cm),
cbind(Am+Cm, Am+Cm+Em)
),
name="expCovMZm"
),
mxAlgebra(
rbind(
cbind(Af+Cf+Ef, Af+Cf),
cbind(Af+Cf, Af+Cf+Ef)
),
name="expCovMZf"
),
# Expected covariance matrix in DZ, for males and females
mxAlgebra(
rbind(
cbind(Am+Cm+Em, (0.5 %x% Am)+Cm),
cbind((0.5 %x% Am)+Cm, Am+Cm+Em)
),
name="expCovDZm"
),
mxAlgebra(
rbind(
cbind(Af+Cf+Ef, (0.5 %x% Af)+Cf),
cbind((0.5 %x% Af)+Cf, Af+Cf+Ef)
),
name="expCovDZf"
),
# Expected covariance matrix in opposite sex pairs, note use of rg and rc
mxAlgebra(
rbind(
cbind(Am+Cm+Em, (am%*%rg%*%t(af) + cm%*%rc%*%t(cf))),
cbind((af%*%rg%*%t(am) + cf%*%rc%*%t(cm)), Af+Cf+Ef)
),
name="expCovOS"
)
),
mxModel("MZm",
mxData(observed=mzmData , type="raw"),
mxFIMLObjective(covariance="ACE.expCovMZm", means="ACE.expMeanM", dimnames=selVars)
),
mxModel("DZm",
mxData(observed=dzmData , type="raw"),
mxFIMLObjective(covariance="ACE.expCovDZm", means="ACE.expMeanM", dimnames=selVars)
),
mxModel("MZf",
mxData(observed=mzfData , type="raw"),
mxFIMLObjective(covariance="ACE.expCovMZf", means="ACE.expMeanF", dimnames=selVars)
),
mxModel("DZf",
mxData(observed=dzfData , type="raw"),
mxFIMLObjective(covariance="ACE.expCovDZf", means="ACE.expMeanF", dimnames=selVars)
),
mxModel("DZo",
mxData(observed=dzoData , type="raw"),
mxFIMLObjective(covariance="ACE.expCovOS", means="ACE.expMeanO", dimnames=selVars)
),
mxAlgebra(MZm.objective + DZm.objective + MZf.objective + DZf.objective + DZo.objective, name="minus2sumll"),
mxAlgebraObjective("minus2sumll"),
mxCI(c('ACE.h2m', 'ACE.c2m', 'ACE.e2m', 'ACE.h2f', 'ACE.c2f', 'ACE.e2f', 'ACE.rg'))
)
unSexACEFit<-mxRun(unSexACEModel)
unSexACEFit<-mxRun(unSexACEModel, intervals=T)
univACESumm <- summary(unSexACEFit)
summary(unSexACEFit)
Submodel 1: Test of qualitative sex differences
Genetic correlation fixed to .5 in opposite-sex pairs
unSexACESub1Model <- mxModel(unSexACEFit, name="unSexACESub1",
mxModel(unSexACEFit$ACE,
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=0, label="Rg", lbound=0, ubound=0.5, name="rg" ),
mxCI(c('ACE.h2m', 'ACE.c2m', 'ACE.e2m', 'ACE.h2f', 'ACE.c2f', 'ACE.e2f', 'ACE.rg'))
))
unSexACESub1Fit <- mxRun(unSexACESub1Model)
unSexACESub1Fit <- mxRun(unSexACESub1Model, intervals=T)
univACESumm <- summary(unSexACESub1Model)
summary(unSexACESub1Model)
Submodel 2: Test of quantitative sex differences
Equate males and female path coefficients
unSexACESub2Model <- mxModel(unSexACEFit, name="unSexACESub2",
mxModel(unSexACEFit$ACE,
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="a11", name="am" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="c11", name="cm" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="e11", name="em" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="a11", name="af" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="c11", name="cf" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="e11", name="ef" )
)
)
unSexACESub2Fit <- mxRun(unSexACESub2Model)
univACESumm <- summary(unSexACESub2Model)
summary(unSexACESub2Model)
Submodel 3: Test both quantitative and qualitative sex differences
unSexACESub3Model <- mxModel(unSexACEFit, name="unSexACESub3",
mxModel(unSexACEFit$ACE,
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="a11", name="am" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="c11", name="cm" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="e11", name="em" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="a11", name="af" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="c11", name="cf" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="e11", name="ef" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=0.5, label="Rg", name="rg" )
),
mxCI(c('ACE.h2m', 'ACE.c2m', 'ACE.e2m', 'ACE.h2f', 'ACE.c2f', 'ACE.e2f', 'ACE.rg'))
)
unSexACESub3Fit <- mxRun(unSexACESub3Model)
unSexACESub3Fit <- mxRun(unSexACESub3Model, intervals=T)
univACESumm <- summary(unSexACESub3Model)
summary(unSexACESub3Model)
mxCompare(unSexACEFit, list(unSexACESub1Fit, unSexACESub2Fit, unSexACESub3Fit))
The model doesn't throw any errors, however, I'm getting back this weird output
base comparison ep minus2LL df AIC diffLL diffdf p
1 unSexACE 9 22697.06 8288 6121.058 NA NA NA
2 unSexACE unSexACESub1 8 22697.06 8289 6119.061 0.002734385 1 0.95829655
3 unSexACE unSexACESub2 6 22701.43 8291 6119.430 4.371689555 3 0.22402482
4 unSexACE unSexACESub3 5 22707.08 8292 6123.080 10.021437535 4 0.04006811
It appears that both full heterogeneity model as well as the qualitative differences model have the same -2LL value, which should not be the case since we are dropping a parameter in the letter (I am told).
What is it that I'm doing wrong? Are my starting value set incorrectly? I tried various values and the same output is repeated.
Help will be much appreciated!
Thank you.