#' Purcell (2002) Bivariate GxE model: Suitable when twins differ on the moderator.
#'
#' GxE interaction models test the hypothesis that the strength of genetic
#' and environmental influences vary parametrically across levels of a measured environment.
#'
#' Whereas univariate [umxGxE()] models assume the twins share the moderator,
#' or have zero correlation on the moderator, [umxGxEbiv()] allows testing moderation in
#' cases where members of a twin pair differ on the moderator, (Purcell, 2002; van der Sluis et al., 2012).
#'
#' This is the same model we teach at Boulder.
#'
#' The following figure shows this bivariate GxE model as a path diagram (Twin 1 shown). Whereas
#' the univariate model incorporates the moderator in the means model, the bivariate model incorporates
#' the moderator as a first class variable, with its own ACE structure, shared pathways to the trait of interest,
#' and the ability to moderate both specific and shared A, C, and E, influences on the trait of interest.
#'
#'
#' ![](GxEbiv.png)
#'
#' Twin 1 and twin 2 A, C, and E latent traits are connected in the standard fashion, with the
#' covariance of the T1 and T2 latent genetic traits set to .5 for DZ and 1.0 for MZ pairs.
#' For the sake of clarity, C, and E paths are omitted here. These mirror those for A.
#'
#' @param name The name of the model (defaults to "GxEbiv")
#' @param selDVs The dependent variable (e.g. IQ)
#' @param selDefs The definition variable (e.g. socioeconomic status)
#' @param sep Expand variable base names, i.e., "_T" makes var -> var_T1 and var_T2
#' @param dzData The DZ dataframe containing the Twin 1 and Twin 2 DV and moderator (4 columns)
#' @param mzData The MZ dataframe containing the Twin 1 and Twin 2 DV and moderator (4 columns)
#' @param lboundACE If !NA, then lbound the main effects at this value (default = NA)
#' @param lboundM If !NA, then lbound the moderators at this value (default = NA)
#' @param dropMissingDef Whether to automatically drop missing def var rows for the user (gives a warning) default = FALSE
#' @param autoRun Whether to run the model (default), or just to create it and return without running.
#' @param tryHard Default ('no') uses normal mxRun. "yes" uses mxTryHard. Other options: "ordinal", "search"
#' @param optimizer Optionally set the optimizer (default NULL does nothing)
#' @return - GxEbiv [mxModel()]
#' @export
#' @md
#' @family Twin Modeling Functions
#' @seealso - [plot()], [umxSummary()], [umxReduce()]
#' @references
#' - Purcell, S. (2002). Variance components models for gene-environment interaction in twin analysis. \emph{Twin Research},
#' \strong{6}, 554-571. \doi{10.1375/twin.5.6.554}.
#'
#' - van der Sluis, S., Posthuma, D., & Dolan, C. V. (2012). A note on false positives and
#' power in G x E modelling of twin data. \emph{Behavior Genetics},
#' \strong{42}, 170-186. \doi{10.1007/s10519-011-9480-3}.
#'
#' @examples
#' require(umx)
#' data(twinData)
#' selDVs = "wt"
#' selDefs = "ht"
#' df = umx_scale_wide_twin_data(twinData, varsToScale = c("ht", "wt"), sep = "")
#' mzData = subset(df, zygosity %in% c("MZFF", "MZMM"))
#' dzData = subset(df, zygosity %in% c("DZFF", "DZMM", "DZOS"))
#'
#' \dontrun{
#' m1 = umxGxEbiv(selDVs = selDVs, selDefs = selDefs,
#' dzData = dzData, mzData = mzData, sep = "", dropMissingDef = TRUE)
#'
#' # Plot Moderation
#' umxSummaryGxEbiv(m1)
#' umxSummary(m1, location = "topright")
#' umxSummary(m1, separateGraphs = FALSE)
#' m2 = umxModify(m1, update = c("cBeta2_r1c1", "eBeta1_r1c1", "eBeta2_r1c1"), comparison = TRUE)
#'
#' # TODO: teach umxReduce to test all relevant hypotheses for umxGxEbiv
#' umxReduce(m1)
#' }
umxGxEbiv <- function(name = "GxEbiv", selDVs, selDefs, dzData, mzData, sep = NULL, lboundACE = 0, lboundM = NA, dropMissingDef = FALSE, autoRun = getOption("umx_auto_run"), tryHard = c("no", "yes", "ordinal", "search"), optimizer = NULL) {
tryHard = match.arg(tryHard)
nSib = 2
if(!is.null(optimizer)){umx_set_optimizer(optimizer)}
if(!is.null(sep)){
if(length(sep) > 1){
stop("sep should be just one word, like '_T'. I will add 1 and 2 afterwards... \n",
"i.e., you have to name your variables 'obese_T1' and 'obese_T2' etc.")
}
selDVs = umx_paste_names(selDVs, sep, 1:nSib)
selDefs = umx_paste_names(selDefs, sep, 1:nSib)
}
if(any(selDefs %in% selDVs)){
warning("selDefs was found in selDVs: You probably gave me all the variables in selDVs instead of just the DEPENDENT variable")
}
if(length(selDVs)/nSib != 1){
umx_msg(selDVs)
stop("selDVs list must be 1 variable... You tried ", length(selDVs)/nSib, " or perhaps you didn't set sep= ?")
}
if(length(selDefs) != 2){
umx_msg(selDefs)
warning("selDefs must be length = 2. Or perhaps you didn't set sep = ?")
}
umx_check_names(selDVs, mzData)
umx_check_names(selDVs, dzData)
message("selDVs: ", omxQuotes(selDVs))
umx_check(!umx_is_cov(dzData, boolean = TRUE), "stop", "Data must be raw for umxGxEbiv")
# TODO umxGxEbiv Check Defs are not correlated 1 or 0
obsMean = mean(colMeans(mzData[,selDVs], na.rm = TRUE)); # Just one average mean for all twins
nVar = length(selDVs)/nSib; # number of dependent variables ** per INDIVIDUAL ( so times-2 for a family)**
rawVar = diag(var(mzData[,selDVs], na.rm = TRUE))[1]
startMain = sqrt(c(.8, .0 ,.6) * rawVar)
# selVars is used to set the dimnames for the expectation and means ordering, so in this hand-assembled function order matters...
# 2021-01-19: was selVars = c(selDVs, selDefs)
# expMean variable order = m_mod, m_trait, m_mod, m_trait
# expectations have dimnames = selVars
selVars = c(selDefs[1], selDVs[1], selDefs[2], selDVs[2]) #MN
# drop any unused variables
dzData = dzData[ , selVars]
mzData = mzData[ , selVars]
mzData = xmu_data_missing(mzData, selVars = selDefs, dropMissingDef=dropMissingDef, hint="mzData")
dzData = xmu_data_missing(dzData, selVars = selDefs, dropMissingDef=dropMissingDef, hint="dzData")
model = mxModel(name,
mxModel("top",
umxMatrix("a11" , "Lower", nrow = 1, ncol = 1, free = TRUE, values = .6),
umxMatrix("c11" , "Lower", nrow = 1, ncol = 1, free = TRUE, values = .6),
umxMatrix("e11" , "Lower", nrow = 1, ncol = 1, free = TRUE, values = .6),
umxMatrix("a21" , "Lower", nrow = 1, ncol = 1, free = TRUE, values = .6),
umxMatrix("c21" , "Lower", nrow = 1, ncol = 1, free = TRUE, values = .6),
umxMatrix("e21" , "Lower", nrow = 1, ncol = 1, free = TRUE, values = .6),
umxMatrix("a22" , "Lower", nrow = 1, ncol = 1, free = TRUE, values = .6),
umxMatrix("c22" , "Lower", nrow = 1, ncol = 1, free = TRUE, values = .6),
umxMatrix("e22" , "Lower", nrow = 1, ncol = 1, free = TRUE, values = .6),
umxMatrix("aBeta1", "Lower", nrow = 1, ncol = 1, free = TRUE, values = .0),
umxMatrix("cBeta1", "Lower", nrow = 1, ncol = 1, free = TRUE, values = .0),
umxMatrix("eBeta1", "Lower", nrow = 1, ncol = 1, free = TRUE, values = .0),
umxMatrix("aBeta2", "Lower", nrow = 1, ncol = 1, free = TRUE, values = .0),
umxMatrix("cBeta2", "Lower", nrow = 1, ncol = 1, free = TRUE, values = .0),
umxMatrix("eBeta2", "Lower", nrow = 1, ncol = 1, free = TRUE, values = .0),
# Assemble Cholesky decomposition for twin 1 and twin 2
umxMatrix("PsAmz", "Symm", nrow = 4, ncol = 4, free = FALSE, values = c(1, 0, 1.0, 0, 1, 0, 1.0, 1, 0, 1)),
umxMatrix("PsAdz", "Symm", nrow = 4, ncol = 4, free = FALSE, values = c(1, 0, 0.5, 0, 1, 0, 0.5, 1, 0, 1)),
umxMatrix("PsC", "Symm", nrow = 4, ncol = 4, free = FALSE, values = c(1, 0, 1.0, 0, 1, 0, 1.0, 1, 0, 1)),
umxMatrix("expMean", "Full", nrow = 1, ncol = 4, free = TRUE, values = 0, labels = c("m_mod", "m_trait", "m_mod", "m_trait"))
),
mxModel("MZ",
# Matrices generated to hold A and E computed Variance Components
# This is a Cholesky decomposition of A for twin 1 and twin 2 (mz and dz)
# note that mod1 appears in the top left part (twin 1) and mod2 in the bottom right part (twin 2)
# Definition variables to create moderated paths (M -> T)
umxMatrix("mod1", "Full", nrow = 1, ncol = 1, free = FALSE, labels = paste0("data.", selDefs[1])),
umxMatrix("mod2", "Full", nrow = 1, ncol = 1, free = FALSE, labels = paste0("data.", selDefs[2])),
mxAlgebra(name = "chA",
rbind(
cbind(top.a11, 0, 0, 0),
cbind(top.a21 + (mod1 %x% top.aBeta1), top.a22 + (mod1 %x% top.aBeta2), 0, 0),
cbind(0, 0, top.a11, 0),
cbind(0, 0, top.a21 + (mod2 %x% top.aBeta1), top.a22 + (mod2 %x% top.aBeta2)))
),
mxAlgebra(name = "chC",
rbind(
cbind(top.c11, 0, 0, 0),
cbind(top.c21 + (mod1 %x% top.cBeta1), top.c22 + (mod1 %x% top.cBeta2), 0, 0),
cbind(0, 0, top.c11, 0),
cbind(0, 0, top.c21 + (mod2 %x% top.cBeta1), top.c22 + (mod2 %x% top.cBeta2)))
),
mxAlgebra(name = "chE",
rbind(
cbind( top.e11, 0, 0, 0),
cbind(top.e21 + (mod1 %x% top.eBeta1), top.e22 + (mod1 %x% top.eBeta2), 0, 0),
cbind( 0, 0, top.e11, 0),
cbind( 0, 0, top.e21 + (mod2 %x% top.eBeta1), top.e22 + (mod2 %x% top.eBeta2)))
),
mxAlgebra(name = "Amz", chA %&% top.PsAmz),
mxAlgebra(name = "C", chC %&% top.PsC),
mxAlgebra(name = "E", chE %*% t(chE)),
mxAlgebra(name = "expCovMZ", Amz + C + E),
mxData(mzData, type = "raw"),
mxExpectationNormal("expCovMZ", means = "top.expMean", dimnames = selVars), mxFitFunctionML()
),
mxModel("DZ",
umxMatrix("mod1", "Full", nrow = 1, ncol = 1, free = FALSE, labels = paste0("data.", selDefs[1])), # "data.defmod1"
umxMatrix("mod2", "Full", nrow = 1, ncol = 1, free = FALSE, labels = paste0("data.", selDefs[2])), # "data.defmod2"
mxAlgebra(name = "chA",
rbind(
cbind(top.a11, 0, 0, 0),
cbind(top.a21 + (mod1 %x% top.aBeta1), top.a22 + (mod1 %x% top.aBeta2), 0, 0),
cbind(0, 0, top.a11, 0),
cbind(0, 0, top.a21 + (mod2 %x% top.aBeta1), top.a22 + (mod2 %x% top.aBeta2)))
),
mxAlgebra(name = "chC",
rbind(
cbind(top.c11, 0, 0, 0),
cbind(top.c21 + (mod1 %x% top.cBeta1), top.c22 + (mod1 %x% top.cBeta2), 0, 0),
cbind(0, 0, top.c11, 0),
cbind(0, 0, top.c21 + (mod2 %x% top.cBeta1), top.c22 + (mod2 %x% top.cBeta2)))
),
mxAlgebra(name = "chE",
rbind(
cbind( top.e11, 0, 0, 0),
cbind(top.e21 + (mod1 %x% top.eBeta1), top.e22 + (mod1 %x% top.eBeta2), 0, 0),
cbind( 0, 0, top.e11, 0),
cbind( 0, 0, top.e21 + (mod2 %x% top.eBeta1), top.e22 + (mod2 %x% top.eBeta2)))
),
mxAlgebra(name = "Adz", chA %&% top.PsAdz ),
mxAlgebra(name = "C", chC %&% top.PsC ),
mxAlgebra(name = "E", chE %*% t(chE) ),
mxAlgebra(name = "expCovDZ", Adz + C + E ),
mxData(dzData, type = "raw"),
mxExpectationNormal("expCovDZ", means = "top.expMean", dimnames = selVars), mxFitFunctionML()
),
mxFitFunctionMultigroup(c("MZ", "DZ"))
)
if (!is.na(lboundACE)) {
model = omxSetParameters(model, labels = c("a11_r1c1", "c11_r1c1", "e11_r1c1","a22_r1c1", "c22_r1c1", "e22_r1c1"), lbound = lboundACE)
}
if (!is.na(lboundM)) {
model = omxSetParameters(model, labels = c("a22_r1c1", "c22_r1c1", "e22_r1c1"), lbound = lboundM)
}
model = as(model, "MxModelGxEbiv")
model = xmu_safe_run_summary(model, autoRun = autoRun, tryHard = tryHard)
invisible(model)
}
#' Summarize a bivariate GxE twin model
#'
#' `umxSummaryGxEbiv` summarizes a bivariate moderation model, as returned by [umxGxEbiv()].
#'
#' @aliases umxSummary.MxModelGxEbiv
#' @param model A fitted [umxGxEbiv()] model to summarize
#' @param std Whether to show the standardized model (not implemented! TRUE)
#' @param CIs Confidence intervals (FALSE)
#' @param xlab label for the x-axis of plot
#' @param location default = "topleft"
#' @param comparison mxCompare model with this model if offered up (default = NULL).
#' @param reduce Whether to run and tabulate a complete model reduction...(Defaults to FALSE)
#' @param separateGraphs Std and raw plots in separate graphs? (default = FALSE)
#' @param report markdown or html (html opens in browser)
#' @param digits round to how many digits (default = 2)
#' @param file The name of the dot file to write: NA = none; "name" = use the name of the model
#' @param returnStd Whether to return the standardized form of the model (default = FALSE)
#' @param ... Optional additional parameters
#' @return - optional [mxModel()]
#' @family Twin Modeling Functions
#' @export
#' @seealso - \code{\link{umxGxEbiv}()}, [plot()], [umxSummary()] work for IP, CP, GxE, and ACE models.
#' @references - <https://github.com/tbates/umx>, <https://tbates.github.io>
#' @md
#' @examples
#' data(twinData)
#' df = umx_scale_wide_twin_data(twinData, varsToScale = c("ht", "wt"), sep = "")
#' mzData = subset(df, zygosity %in% c("MZFF", "MZMM"))
#' dzData = subset(df, zygosity %in% c("DZFF", "DZMM", "DZOS"))
#'
#' \dontrun{
#' m1 = umxGxEbiv(selDVs = "wt", selDefs = "ht",
#' dzData = dzData, mzData = mzData, sep = "", dropMissingDef = TRUE)
#' # Plot Moderation
#' umxSummary(m1)
#' umxSummary(m1, location = "topright")
#' umxSummary(m1, separateGraphs = FALSE)
#' }
umxSummaryGxEbiv <- function(model = NULL, digits = 2, xlab = NA, location = "topleft", separateGraphs = FALSE, file = getOption("umx_auto_plot"), comparison = NULL, std = NULL, reduce = FALSE, CIs = NULL, report = c("markdown", "html"), returnStd = NULL,...) {
message("umxSummaryGxEbiv is beta code: You likely might also want to run summary(model) ")
report = match.arg(report)
umx_has_been_run(model, stop = TRUE)
if(any(!is.null(c(returnStd, std, CIs) ))){
message("For GxE, returnStd, extended, std, comparison or CIs are not yet implemented...")
}
xmu_show_fit_or_comparison(model, comparison = comparison, digits = digits)
selDVs = xmu_twin_get_var_names(model, source = "observed")
nVar = length(selDVs)
# TODO umxSummaryGxEbiv these already exist if a_std exists..
# TODO replace all this with umx_standardizeGxEbiv
# Calculate standardized variance components
# chA <- model$MZ.chA$values # Path coefficients
# message("chA")
# umx_print(chA)
tmp = cbind(model$top.aBeta1$values, model$top.cBeta1$values, model$top.eBeta1$values)
tmp = data.frame(tmp)
names(tmp) = c("aBeta1", "cBeta1", "eBeta1")
message("a, c, and e Betas on defVar")
umx_print(tmp, digits = digits)
tmp = cbind(model$top.aBeta2$values, model$top.cBeta2$values, model$top.eBeta2$values)
tmp = data.frame(tmp)
message("a. c, and e Betas on DV")
names(tmp) = c("aBeta2", "cBeta2", "eBeta2")
umx_print(tmp, digits = digits)
# message("Amz")
# print(model$MZ$Amz$result) # chA %*% top.PsAdz %*% t(chA)), # variance component A
# message("Adz")
# umx_print(model$DZ$Adz$result) # chA %*% top.PsAdz %*% t(chA)), # variance component A
# message("C from DZ model")
# umx_print(model$DZ$C$result) # chA %*% top.PsAdz %*% t(chA)), # variance component A
# message("E from DZ model")
# umx_print(model$DZ$E$result) # chA %*% top.PsAdz %*% t(chA)), # variance component A
tmp = rbind(
cbind(model$top$a11$values, NA , model$top$c11$values, NA , model$top$e11$values, NA),
cbind(model$top$a21$values, model$top$a22$values, model$top$c21$values, model$top$c22$values, model$top$e21$values, model$top$e22$values)
)
tmp = data.frame(tmp)
names(tmp) <-c("a1", "a2", "c1", "c2","e1", "e2")
umx_print(tmp, digits = 3)
# umxPlotGxEbiv(model, xlab = xlab, location = location, separateGraphs = separateGraphs)
if(reduce){
stop("Sorry, umxReduce not yet implemented for umxGxEbiv")
# TODO umxReduce not yet implemented for umxGxEbiv!
umxReduce(model = model, report = report)
}
}
#' @export
umxSummary.MxModelGxEbiv <- umxSummaryGxEbiv
#' Plot the results of a GxE univariate test for moderation of ACE components.
#'
#' Plot GxE results (univariate environmental moderation of ACE components).
#' Options include plotting the raw and standardized graphs separately, or in a combined panel.
#' You can also set the label for the x axis (xlab), and choose the location of the legend.
#'
#' @aliases plot.MxModelGxEbiv
#' @param x A fitted [umxGxEbiv()] model to plot
#' @param xlab String to use for the x label (default = NA, which will use the variable name)
#' @param location Where to plot the legend (default = "topleft")
#' see ?legend for alternatives like bottomright
#' @param separateGraphs (default = FALSE)
#' @param ... Optional additional parameters
#' @return None
#' @family Plotting functions
#' @export
#' @seealso - [plot()], [umxSummary()] work for IP, CP, GxE, SAT, and ACE models.
#' @seealso - [umxGxEbiv()]
#' @references - <https://tbates.github.io>
#' @md
#' @examples
#' require(umx)
#' data(twinData)
#' \dontrun{
#' selDVs = "wt"; selDefs = "ht"
#' df = umx_scale_wide_twin_data(twinData, varsToScale = c("ht", "wt"), suffix = "")
#' mzData = subset(df, zygosity %in% c("MZFF", "MZMM"))
#' dzData = subset(df, zygosity %in% c("DZFF", "DZMM", "DZOS"))
#'
#' m1 = umxGxEbiv(selDVs = selDVs, selDefs = selDefs,
#' dzData = dzData, mzData = mzData, sep = "", dropMissingDef = TRUE)
#' # Plot Moderation
#' plot(m1)
#' umxPlotGxEbiv(m1, xlab = "wt", separateGraphs = TRUE, location = "topleft")
#' }
umxPlotGxEbiv <- function(x, xlab = NA, location = "topleft", separateGraphs = FALSE, ...) {
message("umxGxEbiv plot is in early alpha: it isn't feature complete and has bugs as it's just the umxGxE code. Let me know how you'd like this model displayed")
if(class(x)[[1]] != "MxModelGxEbiv"){
stop("The first parameter of umxPlotGxE must be a umxGxEbiv model, you gave me a ", class(x))
}
model = x # to remind us that x has to be a umxGxEbiv model
# get unique values of moderator
mzData = model$MZ$data$observed
dzData = model$DZ$data$observed
selDefs = names(mzData)[3:4]
if(is.na(xlab)){
xlab = selDefs[1]
}
mzdef1 = as.vector(mzData[, selDefs[1]])
mzdef2 = as.vector(mzData[, selDefs[2]])
dzdef1 = as.vector(dzData[, selDefs[1]])
dzdef2 = as.vector(dzData[, selDefs[2]])
allValuesOfDefVar= c(mzdef1, mzdef2, dzdef1, dzdef2)
defVarValues = sort(unique(allValuesOfDefVar))
a11 = model$top.a11$values
a21 = model$top.a21$values
a22 = model$top.a22$values
Ba1 = model$top.aBeta1$values
Ba2 = model$top.aBeta2$values
c11 = model$top.c11$values
c21 = model$top.c21$values
c22 = model$top.c22$values
Bc1 = model$top.cBeta1$values
Bc2 = model$top.cBeta2$values
e11 = model$top.e11$values
e21 = model$top.e21$values
e22 = model$top.e22$values
Be1 = model$top.eBeta1$values
Be2 = model$top.eBeta2$values
Va = (c(a21 + a22) + (defVarValues * c(Ba1 + Ba2)))^2
Vc = (c(c21 + c22) + (defVarValues * c(Bc1 + Bc2)))^2
Ve = (c(e21 + e22) + (defVarValues * c(Be1 + Be2)))^2
Vt = Va + Vc + Ve
out = as.matrix(cbind(Va, Vc, Ve, Vt))
outStd = as.matrix(cbind(Va/Vt, Vc/Vt, Ve/Vt))
# message("Va, Vc, Ve, Vt")
# umx_print(out , digits = 2)
# message("std: Va, Vc, Ve")
# umx_print(outStd, digits = 2)
tmp= data.frame(rbind(
cbind(a11, NA , c11, NA, e11, NA, Ba1, Bc1, Be1),
cbind(a21, a22, c21, c22, e21, e22, Ba2, Bc2, Be2))
)
names(tmp) = c("a1", "a2", "c1", "c2", "e1", "e2", "a_betas", "c_betas", "e_betas")
umx_print(tmp, digits=2)
if(separateGraphs){
print("Outputting two graphs")
}else{
graphics::par(mfrow = c(1, 2)) # one row, two columns for raw and std variance
# par(mfrow = c(2, 1)) # two rows, one column for raw and std variance
}
acergb = c("red", "green", "blue", "black")
graphics::matplot(x = defVarValues, y = out, type = "l", lty = 1:4, col = acergb, xlab = xlab, ylab = "Variance", main= "Raw Moderation Effects")
graphics::legend(location, legend = c("genetic", "shared", "unique", "total"), lty = 1:4, col = acergb)
# legend(location, legend= c("Va", "Vc", "Ve", "Vt"), lty = 1:4, col = acergb)
graphics::matplot(defVarValues, outStd, type = "l", lty = 1:4, col = acergb, ylim = 0:1, xlab = xlab, ylab = "Standardized Variance", main= "Standardized Moderation Effects")
# legend(location, legend= c("Va", "Vc", "Ve"), lty = 1:4, col = acergb)
graphics::legend(location, legend = c("genetic", "shared", "unique"), lty = 1:4, col = acergb)
graphics::par(mfrow = c(1, 1)) # back to black
}
#' @export
plot.MxModelGxEbiv <- umxPlotGxEbiv
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.