I have a couple of questions relating to adding covariates to the twin model and testing for their relevance
When a nested model differs to a previous model, but is worse (ie has a higher -2LL) how do we assess its significance? This happened in adding in a variable as a covariate. As it makes the model worse do we then exclude it as a covariate? For example – with the script for Full DASS the covariate of Developmental when doing the covariate check after constraining the means was a larger -2LL producing a negative difference.
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")
twinSatPsychModel <- mxModel("twinSatPsych",
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" ),
mxAlgebra( expression= interceptMZ + b_sex %x% obs_sex + b_age %x% obs_age + b_marital %x% obs_marital + b_edu %x% obs_edu + b_psych %x% obs_psych, 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" ),
mxAlgebra( expression= interceptDZ + b_sex %x% obs_sex + b_age %x% obs_age + b_marital %x% obs_marital + b_edu %x% obs_edu + b_psych %x% obs_psych, 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"))
twinSatFitPsych <- mxRun(twinSatPsychModel)
(twinSatFitPsychSum <-summary(twinSatFitPsych))
twinSatDevelopModel <- mxModel("twinSatDevelop",
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=-.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=4.026, label="bdevelop", name="b_develop" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label=c("data.DevelopmentalA", "data.DevelopmentalB"),name="obs_develop" ),
mxAlgebra( expression= interceptMZ + b_sex %x% obs_sex + b_age %x% obs_age + b_marital %x% obs_marital + b_edu %x% obs_edu + b_psych %x% obs_psych + b_develop %x% obs_develop, 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" ),
matrix for the regresion coefficent - sex
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=4.026, label="MZ.bdevelop", name="b_develop" ),
mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
label=c("data.DevelopmentalA", "data.DevelopmentalB"),name="obs_develop" ),
mxAlgebra( expression= interceptDZ + b_sex %x% obs_sex + b_age %x% obs_age + b_marital %x% obs_marital + b_edu %x% obs_edu + b_psych %x% obs_psych + b_develop %x% obs_develop, 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"))
twinSatFitDevelop <- mxRun(twinSatDevelopModel)
(twinSatFitDevelopSum <-summary(twinSatFitDevelop))
tableFitStatistics(twinSatFitDevelop, twinSatFitPsych)
mxCompare(twinSatFitDevelop, twinSatFitPsych)