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.
Hi,
Not sure I see the problem, just glancing at the table: all the comparisons are nested, with df decreasing sequentially...
Hello Tim,
Thanks for your reply.
It has been pointed out to me that the -2LL value for the main model one (I compare the submodels to it) is exactly the same for the sub model in which I test for qualitative sex differences in my sample. Both lines - as in your table - lines 1 and 2 show value of 22697.06.
I've been advised to change the starting values to obtain differential -2LL values, despite I tried that I see no change in them.
It is very difficult for me to interpret this output. And it is very possible, due to being so new to twin modelling, that I simply am lacking the experience in seeing the right differences between models.
What would be your interpretation of the above output?
Thank you.
The mxCompare function often does not round things in a way that is informative.
The LL are not the same (as you can see in the LLdif column). There's an ever-so-slight loss of likelihood: Which can be expected when a useless parameter is dropped.
Starting values are very unlikely to alter this: And it would be disturbing if (in the absence of a code Red or other warning about failure to converge) if they did.
Thank you so much for this, I very much appreciate your advice.
I can move forward!
Hi,
I tried adapting the script in this thread to fit an ADE model. When I get to the submodels, I am running into error messages:
> unSexADESub2Model <- mxModel(unSexADEFit, name="unSexADESub2",mxModel(unSexADEFit$ADE,
+ 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="k11", name="kf" ),
+ mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="l11", name="lf" ),
+ mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="m11", name="mf" )
+ ))
Error: You cannot insert an entity named 'ADE' into a model named 'ADE'
How do I fix this?
Thanks!
You've found a bug in mxModel(), where the 'name' argument is being processed after all the other arguments, and it should be processed before the other arguments. Are you certain that you want to create a new submodel within unSexADEFit? It appears that you want to modify the contents of the unSexADEFit, instead of creating a new submodel. But if it is your intent to create a new submodel, then use the following as a workaround:
Thanks for your speedy response. I am still having trouble with the script. I am trying to test a univariate ADE model where I a) set limitations for quantitative and qualitative sex differences and b) run a sub model by dropping D.
I think the first step is to address the sex limitations. I am running into errors when I get to the sex limitation submodels, including:
Error: unexpected symbol in:
"
unSexADESub1Fit"
> unSexADESub1Fit <- mxRun(unSexADESub1Model, intervals=T)
Running unSexADESub1
Warning message:
In model 'unSexADESub1' NPSOL returned a non-zero status code 1. The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. NPSOL was terminated because no further improvement could be made in the merit function (Mx status GREEN).
> univADESumm <- summary(unSexADESub1Model)
> summary(unSexADESub1Model)
AND
>
> unSexADESub2Fit <- mxRun(unSexADESub2Model)
Running unSexADESub2
Error: The free parameter 'a11' has been assigned multiple starting values! See matrix 'untitled7.am' at location (1, 1) and matrix 'a' at location (1, 1)
> univADESumm <- summary(unSexADESub2Model)
> summary(unSexADESub2Model)
Here is my script thus far. Thank you in advance for any suggestions of help you have to offer!
unSexADEModel <- mxModel("unSexADE",
mxModel("ADE",
ADE Model
Matrices declared to store a, d, and e Path Coefficients
TRANSLATING THE CODE ## Three matrices are set up to hold the
path coefficients for each of the sources of variance considered
in the model. They are 1x1 matrices with one element to be
estimated, which is given a .6 starting value. As the matrix
contains only one element, we assigned one label ("a11").
Its name is different from the name of the matrix ("a") which
in turn is different from the name for the R object (pathA) that
will store the matrix with all its arguments.
pathA <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="a11", name="a" )
pathD <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="d11", name="d" )
pathE <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="e11", name="e" )
pathK <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="k11", name="k" )
pathL <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="l11", name="l" )
pathM <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="m11", name="m" )
Matrices generated to hold A, C, and E computed Variance Components
TRANSLATING THE CODE ## Three new matrices are generated by
multiplying the matrix of the path coefficient with its transpose
to generate the variance component, which will be used to construct
the expected covariance matrix
covA <- mxAlgebra( expression=a %% t(a), name="A" )
covD <- mxAlgebra( expression=d %% t(d), name="D" )
covE <- mxAlgebra( expression=e %% t(e), name="E" )
covK <- mxAlgebra( expression=k %% t(k), name="K" )
covL <- mxAlgebra( expression=l %% t(l), name="L" )
covM <- mxAlgebra( expression=m %% t(m), name="M" )
Algebra for expected Mean and Variance/Covariance Matrices in
MZ & DZ twins
TRANSLATING THE CODE ## By using one label for the 1x2 matrix,
the two elements are constrained to be equal. This matrix will be
used in the MZ and the DZ model, thereby constraining the means
to be constrained to be equal across zygosity.
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values= 20, label="mean", name="expMean" )
TRANSLATING THE CODE ## The expected covariance matrix is
constructed by combining the elements using rbind and cbind functions.
The first element contains the expectation for the variance of twin 1
(sum of the three variance components). Note that the expectation for
the variance of twin 2 (element 2,2) is identical. The expectation
for the covariance is the off-diagonal element.
covMZm <- mxAlgebra( expression=
rbind( cbind(A+D+E , A+D),
cbind(A+D , A+D+E)), name="expCovMZm" )
covMZf <- mxAlgebra( expression=
rbind( cbind(K+L+M , K+M),
cbind(K+L , K+L+M)), name="expCovMZf" )
covDZm <- mxAlgebra( expression=
rbind( cbind(A+D+E, 0.5%x%A+0.25%x%D),
cbind(0.5%x%A+0.25%x%D , A+D+E)), name="expCovDZm" )
covDZf <- mxAlgebra( expression=
rbind( cbind(K+L+M, 0.5%x%K+0.25%x%L),
cbind(0.5%x%K+0.25%x%L , K+L+M)), name="expCovDZf" )
covDZmf <- mxAlgebra( expression=
rbind( cbind(A+D+E, 0.5%x%A+0.25%x%D),
cbind(0.5%x%K+0.25%x%L , K+L+M)), name="expCovDZmf" )
covDZfm <- mxAlgebra( expression=
rbind( cbind(K+L+M, 0.5%x%K+0.25%x%L),
cbind(0.5%x%A+0.25%x%D , A+D+E)), name="expCovDZfm" )
Data objects for Multiple Groups
dataMZf <- mxData( observed=mzfData, type="raw" )
dataDZf <- mxData( observed=dzfData, type="raw" )
dataMZm <- mxData( observed=mzmData, type="raw" )
dataDZm <- mxData( observed=dzmData, type="raw" )
dataDZmf <- mxData( observed=dzmfData, type="raw" )
dataDZfm <- mxData( observed=dzfmData, type="raw" )
Objective objects for Multiple Groups
objMZf <- mxFIMLObjective( covariance="expCovMZf", means="expMean", dimnames=selVars )
objDZf <- mxFIMLObjective( covariance="expCovDZf", means="expMean", dimnames=selVars )
objMZm <- mxFIMLObjective( covariance="expCovMZm", means="expMean", dimnames=selVars )
objDZm <- mxFIMLObjective( covariance="expCovDZm", means="expMean", dimnames=selVars )
objDZmf <- mxFIMLObjective( covariance="expCovDZmf", means="expMean", dimnames=selVars )
objDZfm <- mxFIMLObjective( covariance="expCovDZfm", means="expMean", dimnames=selVars )
Combine Groups
TRANSLATING THE CODE ## All the objects used in the algebra
for the expected covariance matrices now have to be included
into both the MZ and DZ models, so they are combined in a list
named 'pars', which is then included in the 'modelMZ', 'modelDZ'
and 'AdeModel' objects.
parsM <- list( pathA, pathD, pathE, covA, covD, covE )
parsF <- list( pathK, pathL, pathM, covK, covL, covM )
modelMZf <- mxModel( parsF, meanG, covMZf, dataMZf, objMZf, name="MZf" )
modelDZf <- mxModel( parsF, meanG, covDZf, dataDZf, objDZf, name="DZf" )
modelMZm <- mxModel( parsM, meanG, covMZm, dataMZm, objMZm, name="MZm" )
modelDZm <- mxModel( parsM, meanG, covDZm, dataDZm, objDZm, name="DZm" )
modelDZmf <- mxModel( parsM, parsF, meanG, covDZmf, dataDZmf, objDZmf, name="DZmf" )
modelDZfm <- mxModel( parsF, parsM, meanG, covDZfm, dataDZfm, objDZfm, name="DZfm" )
minus2ll <- mxAlgebra( MZf.objective+ DZf.objective+ MZm.objective+ DZm.objective+ DZmf.objective+ DZfm.objective, name="m2LL" )
obj <- mxAlgebraObjective( "m2LL" )
unSexADEModel <- mxModel( "ADE", parsM, parsF, modelMZm, modelMZf, modelDZm, modelDZf, modelDZmf, modelDZfm,
minus2ll, obj )
unSexADEFit<-mxRun(unSexADEModel)
unSexADEFit<-mxRun(unSexADEModel, intervals=T)
univADESumm <- summary(unSexADEFit)
summary(unSexADEFit)
Submodel 1: Test of qualitative sex differences
Genetic correlation fixed to .5 in opposite-sex pairs
unSexADESub1Model <- mxRename(unSexADEFit, "unSexADESub1")
unSexADESub1Model <- mxModel(unSexADESub1Model, mxModel(unSexADEFit$ADE,
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=0, label="Rg", lbound=0, ubound=0.5, name="rg" ),
mxCI(c('ADE.a2m', 'ADE.d2m', 'ADE.e2m', 'ADE.k2f', 'ADE.l2f', 'ADE.m2f', 'ADE.rg')
))
unSexADESub1Fit <- mxRun(unSexADESub1Model)
unSexADESub1Fit <- mxRun(unSexADESub1Model, intervals=T)
univADESumm <- summary(unSexADESub1Model)
summary(unSexADESub1Model)
Submodel 2: Test of quantitative sex differences
Equate males and female path coefficients
unSexADESub2Model <- mxRename(unSexADEFit, "unSexADESub2")
unSexADESub2Model <- mxModel(unSexADESub2Model, mxModel(unSexADESub1Fit$ADE,
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="a11", name="am" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="c11", name="cm" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="e11", name="em" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="k11", name="kf" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="l11", name="lf" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="m11", name="mf" )
))
unSexADESub2Fit <- mxRun(unSexADESub2Model)
univADESumm <- summary(unSexADESub2Model)
summary(unSexADESub2Model)
Submodel 3: Test both quantitative and qualitative sex differences
unSexADESub3Model <- mxRename(unSexADEFit, "unSexADESub3")
unSexADESub3Model <-mxModel(unSexADESub3Model, mxModel(unSexADESub2Fit$ADE,
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="a11", name="am" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="c11", name="cm" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="e11", name="em" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="k11", name="af" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="l11", name="cf" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="m11", name="ef" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=0.5, label="Rg", name="rg" )
),
mxCI(c('ADE.a2m', 'ADE.c2m', 'ADE.e2m', 'ADE.k2f', 'ADE.l2f', 'ADE.m2f', 'ADE.rg'))
)
unSexADESub3Fit <- mxRun(unSexADESub3Model)
unSexADESub3Fit <- mxRun(unSexADESub3Model, intervals=T)
univADESumm <- summary(unSexADESub3Model)
summary(unSexADESub3Model)
mxCompare(unSexADEFit, list(unSexADESub1Fit, unSexADESub2Fit, unSexADESub3Fit))
It looks like the previous expression above the
mxRun()
is not written properly, the "+" sign suggests that the previous line is extended onto the next line: