Checking for the assumption of equal variances

6 replies [Last post]
koak's picture
Offline
Joined: 05/10/2012

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.

1. 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)

2. 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)?

3. 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)

Ryne's picture
Offline
Joined: 07/31/2009
It's pretty odd that you

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.

koak's picture
Offline
Joined: 05/10/2012
Thanks for that thought. Our

Thanks for that thought. Our assignment to twin 1 and 2 was not random as it was based on their birth order...... any thoughts?

Ryne's picture
Offline
Joined: 07/31/2009
Then it appears you've found

Then it appears you've found an effect of birth order. Is that relevant to your research?

neale's picture
Offline
Joined: 07/31/2009
Normality

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.

koak's picture
Offline
Joined: 05/10/2012
Thank you so, so much. This

Thank you so, so much. This has fixed the problem totally!

Much indebted!

mspiegel's picture
Offline
Joined: 07/31/2009
For next time, please include

For next time, please include your script as an attachment. It will be easier to read the replies to the post.