You are here

Checking for the assumption of equal variances

7 posts / 0 new
Last post
koak's picture
Offline
Joined: 05/10/2012 - 23:13
Checking for the assumption of equal variances

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)

mspiegel's picture
Offline
Joined: 07/31/2009 - 15:24
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.

Ryne's picture
Offline
Joined: 07/31/2009 - 15:12
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 - 23:13
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 - 15:12
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 - 15:14
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 - 23:13
Thank you so, so much. This

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

Much indebted!