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)
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?
The mxCompare results in a negative -2LL. I am just not sure how to interpret that!
I haven't read the full post, but isn't the model where you set the covariate path to zero nested inside the model where this is free (not vice versa)
perhaps reverse the model-order in mxCompare()
The order for mxCompare is base, then comparison, where base is the most saturated model (most parameters, fewest degrees of freedom) and the list of comparison models has fewer parameters and more degrees of freedom. If you have a negative -2LL difference and a negative diffdf, then change the order.
If you added parameters but the fit got worse (i.e., diffdf is positive, diffLL is negative), then one or both of your models is either misspecified or didn't converge. If that's the case, post your output and we'll see if we can diagnose it.