You are here

Adding covariates to a model

5 posts / 0 new
Last post
koak's picture
Offline
Joined: 05/10/2012 - 23:13
Adding covariates to a model

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)

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
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?

koak's picture
Offline
Joined: 05/10/2012 - 23:13
The mxCompare results in a

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

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
mxCompare(b,a) instead of mxCompare(a,b)

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()

Ryne's picture
Offline
Joined: 07/31/2009 - 15:12
The order for mxCompare is

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.