Checking for the assumption of equal variances
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)
For next time, please include
Log in or register to post comments
It's pretty odd that you
Log in or register to post comments
In reply to It's pretty odd that you by Ryne
Thanks for that thought. Our
Log in or register to post comments
In reply to Thanks for that thought. Our by koak
Then it appears you've found
Log in or register to post comments
In reply to Then it appears you've found by Ryne
Normality
Log in or register to post comments
In reply to Normality by neale
Thank you so, so much. This
Much indebted!
Log in or register to post comments