Hello
I have been checking whether the equal variances assumption can be applied to my data for twins 1 and 2 and for MZ and DZ twins. It appears that I can't make this assumption.
-
If I am not able to constrain the variances for twins 1 and 2, how do I deal with this? Do I just leave the model so that the variances for each twin are included separately in creating the covariance (ie as in the model below)? (script below)
-
If I come up with above (ie a significant difference in the -2LL for constraining of variances for twins 1 and 2 compared with the last model), how is it best to look at whether or not MZ/DZ variances can be constrained? Do we still just compare with the constraining of variances for Twin 1/Twin 2? Or go back to compare with the constrainining of the means (ie model twinSatSub2 below)?
-
If I needed to have the variances for MZ/DZ twins different how do I do this? Is it again like with the model used to set this up (see below) of keeping the MZ and DZ variances specified separately?
Many, many thanks!
Script
DASS <- read.csv("DASS_Data.csv", header=TRUE)
mzData <- as.data.frame(subset(DASS, Zygosity==0 , names(DASS)))
dzData <- as.data.frame(subset(DASS, Zygosity==1 , names(DASS)))
stMean<-mean(mzData$DASSA, na.rm=T)+ sd(mzData$DASSA, na.rm=T)
stVar<-var(mzData$DASSA, na.rm=T)+ sd(mzData$DASSA, na.rm=T)
stCovMZ<-stVar.4
stCovDZ<-stVar.2
nv <- 1
selVars<-c("DASSA","DASSB")
twinSatSphereModel <- mxModel("twinSatSphere",
mxModel("MZ",
mxMatrix( type="Full", nrow=nv, ncol=2, free=TRUE,
values=stMean, label=c("MZm1", "MZm2"), name="interceptMZ" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=-0.058, label="bsex", name="b_sex" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label="data.Sex", name="obs_sex" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=-.124, label="bAge", name="b_age" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label=c("data.AgeA", "data.AgeB"), name="obs_age" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=-.052, label="bedu", name="b_edu" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label=c("data.EducationA", "data.EducationB"), name="obs_edu" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=1.041, label="bmarital", name="b_marital" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label=c("data.MaritalA", "data.MaritalB"), name="obs_marital" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=10.091, label="bpsych", name="b_psych" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label=c("data.PsychexcludeA", "data.PsychexcludeB"), name="obs_psych" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=14.538, label="bsphere", name="b_sphere" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label=c("data.SphereOnlyA", "data.SphereOnlyB"), name="obs_sphere" ),
mxAlgebra( expression= interceptMZ + b_sex %x% obs_sex + b_age %x% obs_age + b_edu %x% obs_edu + b_marital %x% obs_marital + b_psych %x% obs_psych + b_sphere %x% obs_sphere, name="expMeanMZ" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=stVar,
label=c("VMZ1"),name="VarMZtw1" ),
mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=stVar,
label=c("VMZ2"),name="VarMZtw2"),
mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=stCovMZ, name="CovMZ"),
mxAlgebra( expression= rbind ( cbind(VarMZtw1, t(CovMZ)),
cbind(CovMZ, VarMZtw2) ), name="expCovMZ" ),
mxData(observed=mzData, type="raw"),
mxFIMLObjective(covariance="expCovMZ", means="expMeanMZ", dimnames=selVars) ),
mxModel("DZ",
mxMatrix( type="Full", nrow=nv, ncol=2, free=TRUE,
values=stMean, label=c("DZm1", "DZm2"), name="interceptDZ" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=-0.058, label="MZ.bsex", name="b_sex" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label="data.Sex", name="obs_sex" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=-.124, label="MZ.bAge", name="b_age" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label=c("data.AgeA", "data.AgeB"), name="obs_age" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=-.052, label="MZ.bedu", name="b_edu" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label=c("data.EducationA", "data.EducationB"), name="obs_edu" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=1.041, label="MZ.bmarital", name="b_marital" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label=c("data.MaritalA", "data.MaritalB"), name="obs_marital" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=10.091, label="MZ.bpsych", name="b_psych" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label=c("data.PsychexcludeA", "data.PsychexcludeB"), name="obs_psych" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=14.538, label="MZ.bsphere", name="b_sphere" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label=c("data.SphereOnlyA", "data.SphereOnlyB"), name="obs_sphere" ),
mxAlgebra( expression= interceptDZ + b_sex %x% obs_sex + b_age %x% obs_age + b_edu %x% obs_edu + b_marital %x% obs_marital + b_psych %x% obs_psych + b_sphere %x% obs_sphere, name="expMeanDZ" ),
mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=stVar,
label=c("VDZ1"), name="VarDZtw1"),
mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=stVar,
label=c("VDZ2"), name="VarDZtw2"),
mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=stCovDZ, name="CovDZ"),
mxAlgebra( expression= rbind ( cbind(VarDZtw1, t(CovDZ)),
cbind(CovDZ, VarDZtw2) ), name="expCovDZ" ),
mxData(observed=dzData, type="raw"),
mxFIMLObjective(covariance="expCovDZ", means="expMeanDZ", dimnames=selVars)),
mxAlgebra(MZ.objective + DZ.objective, name="min2LL"),
mxAlgebraObjective("min2LL"))
twinSatSphereFit <- mxRun(twinSatSphereModel)
(twinSatFitSphereSum <-summary(twinSatSphereFit))
Specify and Run Sub-Model 1 equating means across twin order within zygosity group
twinSatModelSub1 <- (twinSatSphereModel)
twinSatModelSub1$MZ$interceptMZ <- mxMatrix(type="Full", nrow=1, ncol=2, free=T, values=stVar, label="mMZ", name="interceptMZ") #testing whether twin 1 mean MZ equals twin 2 meanMZ;
twinSatModelSub1$DZ$interceptDZ <- mxMatrix(type="Full", nrow=1, ncol=2, free=T, values=stVar, label="mDZ", name="interceptDZ")
twinSatFitSub1 <- mxRun(twinSatModelSub1)
twinSatModelSub1Summ <- summary(twinSatFitSub1)
twinSatModelSub1Summ
tableFitStatistics(twinSatSphereFit, twinSatFitSub1)
Specify and Run Sub-Model 2 equating means across twin order AND zygosity group
twinSatModelSub2 <-(twinSatModelSub1)
twinSatModelSub2$MZ$interceptMZ <- mxMatrix(type="Full", nrow=1, ncol=2, free=T, values=stMean, label="mean", name="interceptMZ")
twinSatModelSub2$DZ$interceptDZ <- mxMatrix(type="Full", nrow=1, ncol=2, free=T, values=stMean, label="mean", name="interceptDZ")
twinSatFitSub2 <- mxRun(twinSatModelSub2)
twinSatModelSub2Summ <- summary(twinSatFitSub2)
twinSatModelSub2Summ
tableFitStatistics(twinSatFitSub1, twinSatFitSub2)
tableFitStatistics(twinSatSphereFit, twinSatFitSub2)
Specify and Run Sub-Model 3 equating Variances across twin order in MZ and DZ group
twinSatModelSub3 <-(twinSatModelSub2)# getting the model that best fit the data
twinSatModelSub3$MZ$VarMZtw1 <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=stVar, label=c("VMZ"),name="VarMZtw1")
twinSatModelSub3$MZ$VarMZtw2 <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=stVar, label=c("VMZ"),name="VarMZtw2")
twinSatModelSub3$DZ$VarDZtw1 <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=stVar, label=c("VDZ"),name="VarDZtw1")
twinSatModelSub3$DZ$VarDZtw2 <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=stVar, label=c("VDZ"),name="VarDZtw2")
twinSatFitSub3 <- mxRun(twinSatModelSub3)
twinSatModelSub3Summ <- summary(twinSatFitSub3)
twinSatModelSub3Summ
tableFitStatistics(twinSatFitSub2,twinSatFitSub3)
Specify and Run Sub-Model 4 equating Variances across twin order AND zygosity group
twinSatModelSub4 <-(twinSatModelSub3)
twinSatModelSub4$MZ$VarMZtw1 <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=twicestVar, label=c("V"),name="VarMZtw1")
twinSatModelSub4$MZ$VarMZtw2 <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=twicestVar, label=c("V"),name="VarMZtw2")
twinSatModelSub4$DZ$VarDZtw1 <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=twicestVar, label=c("V"),name="VarDZtw1")
twinSatModelSub4$DZ$VarDZtw2 <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=twicestVar, label=c("V"),name="VarDZtw2")
twinSatFitSub4 <- mxRun(twinSatModelSub4)
twinSatModelSub4Summ <- summary(twinSatFitSub4)
twinSatModelSub4Summ
tableFitStatistics(twinSatFitSub3,twinSatFitSub4)
For next time, please include your script as an attachment. It will be easier to read the replies to the post.
It's pretty odd that you can't constrain variances across twins, as which twin is 1 and which is 2 is usually random. If the assignment is not random (i.e., opposite sex twins always put the male as 1 and the female as 2), then you've discovered a potentially interesting effect. If the assignment truly is random, then this is just an "unhappy randomization", and you should be able to re-assign twin numbers. One of the more experienced BG folk on this forum should probably weigh in on common practice, however.
Thanks for that thought. Our assignment to twin 1 and 2 was not random as it was based on their birth order...... any thoughts?
Then it appears you've found an effect of birth order. Is that relevant to your research?
In my experience, there are few effects of birth order in twins. What usually causes failure of the equal variances of T1 & T2 assumption is non-normality of the data. I recommend you (@koak) test this.
Thank you so, so much. This has fixed the problem totally!
Much indebted!