## Adding covariates to a model

4 replies [Last post]
Offline
Joined: 05/10/2012

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

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)

Offline
Joined: 07/31/2009
doesn't mxCompare give you what you need?

I'd have thought that mxCompare() will give you the significance of your change in -2ll based on the change in df? Also can see the AIC change?

Offline
Joined: 05/10/2012
The mxCompare results in a

The mxCompare results in a negative -2LL. I am just not sure how to interpret that!

Offline
Joined: 07/31/2009