You are here

Multiple Moderators

2 posts / 0 new
Last post
benruk's picture
Offline
Joined: 12/09/2019 - 12:01
Multiple Moderators

Hello,
some days ago I read the paper by Shaun Purcell (2002) "Variance Components Models for Gene-Environment Interaction in Twin Analysis" and I'm especially interested in modelling a three-way interaction (GxExE). So, for example, the genetic effect may be moderated by the socioeconomic status and by age, so that the additive genetic variance component is:

$$
(a+\beta_{Age}M_{Age}+\beta_{SES}M_{SES}+\beta_{AgexSES}M_{Age}M_{SES})^{2}
$$

Since my goal is to model such a three-way interaction, my question is what software I could use to do so. Since Purcell used for his analysis the Mx-Software, it should be possible to model such a three-way interaction with the R package "OpenMx". My question is: Could I use the "umx"-package as well to model a GxExE interaction? In the "umx"-manual, I have read that it is possible to model a two-way interaction, but I am not sure, if it could be extended to model a three-way interaction.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
answer

As to whether or not you can do the moderation model you describe in umx, my guess is the answer is "no", but tbates (umx's maintainer) would be the best person to answer that question. Nonetheless, the model you describe can certainly be done in OpenMx--I have done it myself. Below is a snippet of the OpenMx-specific syntax I used:

defnvars <- mxMatrix(type="Full",nrow=1,ncol=5,free=F,
                     labels=c("data.age1","data.sex1","data.age2","data.sex2","data.ses"),name="xdef")
 
Beta0mz <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=110,labels="beta0mz",name="Beta0mz")
Beta0dz <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=110,labels="beta0mz",name="Beta0dz") #Same by zyg.
Beta0ado <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=110,labels="beta0ado",name="Beta0ado")
Beta0bio <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=110,labels="beta0bio",name="Beta0bio")
Beta0mix <- mxMatrix(type="Full",nrow=1,ncol=2,free=T,values=110,labels=c("beta0ado","beta0bio"),name="Beta0mix")
Betage <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=-1,labels="betage",name="Betage")
Betasex <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=-2,labels="betasex",name="Betasex")
BetasesA <- mxMatrix(type="Full",nrow=1,ncol=1,free=F,values=0,labels="betasesA",name="BetasesA")
BetasesB <- mxMatrix(type="Full",nrow=1,ncol=1,free=F,values=0,labels="betasesB",name="BetasesB")
 
Mzmeans <- mxAlgebra(cbind(Beta0mz+(Betage*xdef[1,1])+(Betasex*xdef[1,2])+(BetasesB*xdef[1,5]), 
                           Beta0mz+(Betage*xdef[1,3])+(Betasex*xdef[1,4])+(BetasesB*xdef[1,5])), name="Mzmeans",
                     dimnames=list(NULL,c("fsiq1","fsiq2")))
Dzmeans <- mxAlgebra(cbind(Beta0dz+(Betage*xdef[1,1])+(Betasex*xdef[1,2])+(BetasesB*xdef[1,5]), 
                           Beta0dz+(Betage*xdef[1,3])+(Betasex*xdef[1,4])+(BetasesB*xdef[1,5])), name="Dzmeans",
                     dimnames=list(NULL,c("fsiq1","fsiq2")))
Adomeans <- mxAlgebra(cbind(Beta0ado+(Betage*xdef[1,1])+(Betasex*xdef[1,2])+(BetasesA*xdef[1,5]), 
                            Beta0ado+(Betage*xdef[1,3])+(Betasex*xdef[1,4])+(BetasesA*xdef[1,5])), name="Adomeans",
                      dimnames=list(NULL,c("fsiq1","fsiq2")))
Biomeans <- mxAlgebra(cbind(Beta0bio+(Betage*xdef[1,1])+(Betasex*xdef[1,2])+(BetasesB*xdef[1,5]), 
                            Beta0bio+(Betage*xdef[1,3])+(Betasex*xdef[1,4])+(BetasesB*xdef[1,5])), name="Biomeans",
                      dimnames=list(NULL,c("fsiq1","fsiq2")))
Mixmeans <- mxAlgebra(cbind(Beta0mix[1,1]+(Betage*xdef[1,1])+(Betasex*xdef[1,2])+(BetasesA*xdef[1,5]), 
                            Beta0mix[1,2]+(Betage*xdef[1,3])+(Betasex*xdef[1,4])+(BetasesB*xdef[1,5])), name="Mixmeans",
                      dimnames=list(NULL,c("fsiq1","fsiq2")))
 
GammaA0 <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=sqrt(0.4*225),labels="gammaA0",name="GammaA0")
GammaC0 <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=sqrt(0.3*225),labels="gammaC0",name="GammaC0")
GammaE0 <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=sqrt(0.3*225),labels="gammaE0",name="GammaE0")
GammaT0 <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=sqrt(0.01*225),labels="gammaT0",name="GammaT0")
#Moderation by Age:
GammaA1 <- mxMatrix(type="Full",nrow=1,ncol=1,free=F,values=0,labels="gammaA1",name="GammaA1")
GammaC1 <- mxMatrix(type="Full",nrow=1,ncol=1,free=F,values=0,labels="gammaC1",name="GammaC1")
GammaE1 <- mxMatrix(type="Full",nrow=1,ncol=1,free=F,values=0,labels="gammaE1",name="GammaE1")
#Moderation by SES:
GammaA2 <- mxMatrix(type="Full",nrow=1,ncol=1,free=F,values=0,labels="gammaA2",name="GammaA2")
GammaC2 <- mxMatrix(type="Full",nrow=1,ncol=1,free=F,values=0,labels="gammaC2",name="GammaC2")
GammaE2 <- mxMatrix(type="Full",nrow=1,ncol=1,free=F,values=0,labels="gammaE2",name="GammaE2")
#Moderation by AgexSES:
GammaA3 <- mxMatrix(type="Full",nrow=1,ncol=1,free=F,values=0,labels="gammaA3",name="GammaA3")
GammaC3 <- mxMatrix(type="Full",nrow=1,ncol=1,free=F,values=0,labels="gammaC3",name="GammaC3")
#raa:
r <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=0.51,lbound=0,ubound=1,labels="raa",name="r")
 
Mzvar <- mxAlgebra(
  rbind(
    cbind(
      (GammaA0+(GammaA1*xdef[1,1])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,1]*xdef[1,5]))^2 +
        (GammaC0+(GammaC1*xdef[1,1])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,1]*xdef[1,5]))^2 +
        (GammaE0+(GammaE1*xdef[1,1])+(GammaE2*xdef[1,5]))^2 + GammaT0^2,
      (GammaA0+(GammaA1*xdef[1,1])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,1]*xdef[1,5]))*(GammaA0+(GammaA1*xdef[1,3])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,3]*xdef[1,5])) +
        (GammaC0+(GammaC1*xdef[1,1])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,1]*xdef[1,5]))*(GammaC0+(GammaC1*xdef[1,3])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,3]*xdef[1,5])) +
        GammaT0^2
    ),
    cbind(
      (GammaA0+(GammaA1*xdef[1,1])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,1]*xdef[1,5]))*(GammaA0+(GammaA1*xdef[1,3])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,3]*xdef[1,5])) +
        (GammaC0+(GammaC1*xdef[1,1])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,1]*xdef[1,5]))*(GammaC0+(GammaC1*xdef[1,3])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,3]*xdef[1,5])) +
        GammaT0^2,
      (GammaA0+(GammaA1*xdef[1,3])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,3]*xdef[1,5]))^2 +
        (GammaC0+(GammaC1*xdef[1,3])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,3]*xdef[1,5]))^2 +
        (GammaE0+(GammaE1*xdef[1,3])+(GammaE2*xdef[1,5]))^2 + GammaT0^2
  )),
  "Mzvar", dimnames=list(c("fsiq1","fsiq2"),c("fsiq1","fsiq2"))
)
 
Dzvar <- mxAlgebra(
  rbind(
    cbind(
      (GammaA0+(GammaA1*xdef[1,1])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,1]*xdef[1,5]))^2 +
        (GammaC0+(GammaC1*xdef[1,1])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,1]*xdef[1,5]))^2 +
        (GammaE0+(GammaE1*xdef[1,1])+(GammaE2*xdef[1,5]))^2 + GammaT0^2,
      (GammaA0+(GammaA1*xdef[1,1])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,1]*xdef[1,5]))*r*(GammaA0+(GammaA1*xdef[1,3])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,3]*xdef[1,5])) +
        (GammaC0+(GammaC1*xdef[1,1])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,1]*xdef[1,5]))*(GammaC0+(GammaC1*xdef[1,3])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,3]*xdef[1,5])) +
        GammaT0^2
    ),
    cbind(
      (GammaA0+(GammaA1*xdef[1,1])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,1]*xdef[1,5]))*r*(GammaA0+(GammaA1*xdef[1,3])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,3]*xdef[1,5])) +
        (GammaC0+(GammaC1*xdef[1,1])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,1]*xdef[1,5]))*(GammaC0+(GammaC1*xdef[1,3])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,3]*xdef[1,5])) +
        GammaT0^2,
      (GammaA0+(GammaA1*xdef[1,3])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,3]*xdef[1,5]))^2 +
        (GammaC0+(GammaC1*xdef[1,3])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,3]*xdef[1,5]))^2 +
        (GammaE0+(GammaE1*xdef[1,3])+(GammaE2*xdef[1,5]))^2 + GammaT0^2
    )),
  "Dzvar", dimnames=list(c("fsiq1","fsiq2"),c("fsiq1","fsiq2"))
)
 
Biovar <- mxAlgebra(
  rbind(
    cbind(
      (GammaA0+(GammaA1*xdef[1,1])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,1]*xdef[1,5]))^2 +
        (GammaC0+(GammaC1*xdef[1,1])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,1]*xdef[1,5]))^2 +
        (GammaE0+(GammaE1*xdef[1,1])+(GammaE2*xdef[1,5]))^2,
      (GammaA0+(GammaA1*xdef[1,1])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,1]*xdef[1,5]))*r*(GammaA0+(GammaA1*xdef[1,3])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,3]*xdef[1,5])) +
        (GammaC0+(GammaC1*xdef[1,1])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,1]*xdef[1,5]))*(GammaC0+(GammaC1*xdef[1,3])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,3]*xdef[1,5]))
    ),
    cbind(
      (GammaA0+(GammaA1*xdef[1,1])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,1]*xdef[1,5]))*r*(GammaA0+(GammaA1*xdef[1,3])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,3]*xdef[1,5])) +
        (GammaC0+(GammaC1*xdef[1,1])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,1]*xdef[1,5]))*(GammaC0+(GammaC1*xdef[1,3])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,3]*xdef[1,5])),
      (GammaA0+(GammaA1*xdef[1,3])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,3]*xdef[1,5]))^2 +
        (GammaC0+(GammaC1*xdef[1,3])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,3]*xdef[1,5]))^2 +
        (GammaE0+(GammaE1*xdef[1,3])+(GammaE2*xdef[1,5]))^2
    )),
  "Biovar", dimnames=list(c("fsiq1","fsiq2"),c("fsiq1","fsiq2"))
)
 
Adovar <- mxAlgebra(
  rbind(
    cbind(
      (GammaA0+(GammaA1*xdef[1,1])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,1]*xdef[1,5]))^2 +
        (GammaC0+(GammaC1*xdef[1,1])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,1]*xdef[1,5]))^2 +
        (GammaE0+(GammaE1*xdef[1,1])+(GammaE2*xdef[1,5]))^2,  
      (GammaC0+(GammaC1*xdef[1,1])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,1]*xdef[1,5])) *
        (GammaC0+(GammaC1*xdef[1,3])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,3]*xdef[1,5]))
    ),
    cbind(
      (GammaC0+(GammaC1*xdef[1,1])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,1]*xdef[1,5])) *
        (GammaC0+(GammaC1*xdef[1,3])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,3]*xdef[1,5])),
      (GammaA0+(GammaA1*xdef[1,3])+(GammaA2*xdef[1,5])+(GammaA3*xdef[1,3]*xdef[1,5]))^2 +
        (GammaC0+(GammaC1*xdef[1,3])+(GammaC2*xdef[1,5])+(GammaC3*xdef[1,3]*xdef[1,5]))^2 +
        (GammaE0+(GammaE1*xdef[1,3])+(GammaE2*xdef[1,5]))^2
    )),
  "Adovar", dimnames=list(c("fsiq1","fsiq2"),c("fsiq1","fsiq2"))
)
 
Mzmodel <- mxModel("MZtwins", Beta0mz, Betage, Betasex,
                   BetasesB,
                   GammaA0,GammaA1,GammaA2,GammaA3,GammaC0,GammaC1,GammaC2,GammaC3,GammaE0,GammaE1,GammaE2,GammaT0,
                   Mzmeans, Mzvar, mzdat, defnvars,
                   mxFIMLObjective(covariance="Mzvar",means="Mzmeans"))
Dzmodel <- mxModel("DZtwins", Beta0dz, Betage, Betasex, r,
                   GammaA0,GammaA1,GammaA2,GammaA3,GammaC0,GammaC1,GammaC2,GammaC3,GammaE0,GammaE1,GammaE2,GammaT0,
                   BetasesB,
                   Dzmeans, Dzvar, dzdat, defnvars,
                   mxFIMLObjective(covariance="Dzvar",means="Dzmeans"))
Biomodel <- mxModel("Biosibs", Beta0bio, Betage, Betasex,
                    GammaA0,GammaA1,GammaA2,GammaA3,GammaC0,GammaC1,GammaC2,GammaC3,GammaE0,GammaE1,GammaE2,GammaT0,
                    BetasesB, r,
                    Biomeans, Biovar, biodat, defnvars,
                    mxFIMLObjective(covariance="Biovar",means="Biomeans"))
Adomodel <- mxModel("Adosibs", Beta0ado, Betage, Betasex, BetasesA,
                    GammaA0,GammaA1,GammaA2,GammaA3,GammaC0,GammaC1,GammaC2,GammaC3,GammaE0,GammaE1,GammaE2,GammaT0,
                    Adomeans, Adovar, adodat, defnvars,
                    mxFIMLObjective(covariance="Adovar",means="Adomeans"))
Mixmodel <- mxModel("Mixsibs", Beta0mix, Betage, Betasex,
                    GammaA0,GammaA1,GammaA2,GammaA3,GammaC0,GammaC1,GammaC2,GammaC3,GammaE0,GammaE1,GammaE2,GammaT0,
                    BetasesA, BetasesB,
                    Mixmeans, Adovar, mixdat, defnvars,
                    mxFIMLObjective(covariance="Adovar",means="Mixmeans"))
 
Model_1 <- mxModel(
  "ACTE", Mzmodel, Dzmodel, Biomodel, Adomodel, Mixmodel,
  mxAlgebra(MZtwins.objective + DZtwins.objective + Adosibs.objective + Biosibs.objective + Mixsibs.objective,
            "AllObj"),
  mxAlgebraObjective("AllObj")
)
run_1 <- mxRun(Model_1)
( summ_1 <- summary(run_1) )

Note that this syntax was written for a sample that included non-twin full sibs and unrelated sibs reared together, meaning that twin-specific ('T') effects and the genetic resemblance of full sibs (raa) could be freely estimated. If you only have twins, you'll have to modify the syntax accordingly. Note also that this syntax was written for OpenMx v1.x; you'll probably want to replace the *Objectives with the appropriate mxFitFunction* and mxExpectation* objects.