Login

News

Multiple Moderators

Posted on
No user picture. benruk Joined: 12/09/2019
Forums
umx
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.

Replied on Wed, 01/08/2020 - 12:54
Picture of user. AdminRobK Joined: 01/24/2014

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.

Log in or register to post comments
Replied on Wed, 01/22/2020 - 10:45
No user picture. benruk Joined: 12/09/2019

Wow, thank you very much for your answer, which surely will help me a lot.
Log in or register to post comments
  Login

News