#' Extends Mendelian randomization with the twin design to test evidence of causality
#'
#' @description
#' Testing causal claims is often difficult due to an inability to conduct experimental randomization of traits and
#' situations to people. When twins are available, even when measured on a single occasion, the pattern of cross-twin
#' cross-trait correlations can (given distinguishable modes of inheritance for the two traits) falsify causal hypotheses.
#'
#' `umxMRDoC` implements a 2-group model to form latent variables for each of two traits, and allows testing whether
#' trait 1 causes trait 2, vice-versa, or even reciprocal causation. This is robust to several types of confounding
#' due to the instrumental variable approach included in the model.
#'
#' This function applies both the MRDoC model and the MRDoC2 model depending on how many PRSs are passed as arguments.
#'
#' @param name The name of the model (defaults to either "MRDoC" or "MRDoC2).
#' @param pheno Phenotypes of interest, order matters ("exposure", "outcome")
#' @param prss Polygenic score(s). If a single one is passed MRDoC is run, MRDoC2 otherwise.
#' @param dzData The DZ dataframe
#' @param mzData The MZ dataframe
#' @param sep The separator in twin variable names, default = "_T", e.g. "dep_T1".
#' @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).
#' @param summary Optionally show a summary.
#' @param data = NULL If building the MZ and DZ datasets internally from a complete data set.
#' @param zyg = "zygosity" (for the data= method of using this function)
#' @return - [mxModel()] of subclass MxModelMRDoC
#' @export
#' @family Twin Modeling Functions
#' @seealso - [umxDoC()]
#' @references - Minica CC, Dolan CV, Boomsma DI, et al. (2018) Extending Causality Tests with Genetic Instruments: An Integration of Mendelian Randomization with the Classical Twin Design. Behavior Genetics 48(4): 337–349. \doi{10.1007/s10519-018-9904-4}
#' * McGue, M., Osler, M., & Christensen, K. (2010). Causal Inference and Observational Research: The Utility of Twins. *Perspectives on Psychological Science*, **5**, 546-556. \doi{10.1177/1745691610383511}
#' * Castro-de-Araujo LFS, Singh M, Zhou Y, et al. (2022) MR-DoC2: Bidirectional Causal Modeling with Instrumental Variables and Data from Relatives. Behavior Genetics. \doi{10.1007/s10519-022-10122-x}
#' @md
#' @examples
#' \dontrun{
#'
#' # ================
#' # = 1. Load Data =
#' # ================
#' data(docData)
#' mzData = subset(docData, zygosity %in% c("MZFF", "MZMM"))
#' dzData = subset(docData, zygosity %in% c("DZFF", "DZMM"))
#'
#' # ============================
#' # = 2. Make a MRDoC2 model =
#' # ============================
#' out = umxMRDoC(mzData = mzData, dzData = dzData,
#' pheno = c("varA1", "varA2"), prss = c("varB1", "varA3") )
#'}
umxMRDoC <- function(pheno, prss, mzData = NULL, dzData = NULL, data = NULL, zyg = NULL, sep = "_T", summary = !umx_set_silent(silent = TRUE), name = NULL, autoRun = getOption("umx_auto_run"), tryHard = c("no", "yes", "ordinal", "search"), optimizer = NULL) {
sketch = FALSE
if (missing(data) && missing(mzData) && missing(dzData)) sketch = TRUE
tryHard = match.arg(tryHard)
# options(mxByrow = TRUE)
umx_set_silent(TRUE)
# Managing data
if (!is.null(data)) {
if ("tbl" %in% class(data)) {
data = as.data.frame(data)
}
mzData = data[data[, zyg] %in% ifelse(is.null(mzData), "MZ", mzData), ]
dzData = data[data[, zyg] %in% ifelse(is.null(dzData), "DZ", dzData), ]
} else {
if ("tbl" %in% class(mzData)) {
mzData = as.data.frame(mzData)
dzData = as.data.frame(dzData)
}
}
vnames = tvars(c(pheno, prss), sep = sep)
if (sketch==F){
xmu_twin_check(
selDVs = c(pheno, prss),
sep = sep, dzData = dzData, mzData = mzData, enforceSep = TRUE,
nSib = 2, optimizer = optimizer
)
}
if (!missing(mzData)) {mzData = xmu_make_mxData(mzData, manifests = vnames)}
if (!missing(dzData)) {dzData = xmu_make_mxData(dzData, manifests = vnames)}
x = y = 1
if (length(prss) == 1) {
if (is.null(name)) { name = "MRDoC"}
top = mxModel("top",
umxMatrix("BE", type = "Full",nrow=3, ncol = 3, byrow = TRUE,
labels = c(NA, "g2", "b1",
"g1", NA, "b2",
NA, NA, NA),
free = c(FALSE, FALSE, TRUE,
TRUE, FALSE, TRUE,
FALSE, FALSE, FALSE)),
umxMatrix('I', type='Iden', nrow= 3,ncol= 3),
umxMatrix('F', type='Full', nrow=5, ncol=6, free=FALSE, byrow = TRUE,
values=c(1,0,0,0,0,0,
0,1,0,0,0,0,
0,0,1,0,0,0,
0,0,0,1,0,0,
0,0,0,0,1,0)),
umxMatrix('LY', type='Full',nrow=3, ncol = 3, free = FALSE, values = diag(3), labels = NA),
umxMatrix('A', type='Symm', nrow=3,ncol = 3,byrow = TRUE,
labels=c("ab2","abraas",NA,
"abraas","as2",NA,
NA,NA,"x2"),
free=c(TRUE,TRUE,FALSE,
TRUE,TRUE,FALSE,
FALSE,FALSE,TRUE)),
umxMatrix('C', type='Symm',nrow=3, ncol = 3,byrow = TRUE,
labels =c("cb2","cbrccs",NA,
"cbrccs","cs2",NA,
NA,NA,NA),
free=c(TRUE,TRUE,FALSE,
TRUE,TRUE,FALSE,
FALSE,FALSE,FALSE)),
umxMatrix('E', type='Symm', nrow=3, ncol = 3,byrow = TRUE,
labels =c("eb2",NA,NA,
NA,"es2",NA,
NA,NA,NA),
free= c(TRUE,FALSE,FALSE,
FALSE,TRUE,FALSE,
FALSE,FALSE,FALSE)),
umxMatrix('dzmu', type='Full', nrow=1, ncol=6, free=TRUE, values= 0, byrow = TRUE,
labels=c('muPh1','muPh2','muPS1','muPh1','muPh2','muPS1')
),
mxAlgebra('mzmu', expression = dzmu%*%t(F)),
mxAlgebra('A_' , expression = solve(I - BE)%&%A),
mxAlgebra('C_' , expression = solve(I - BE)%&%C),
mxAlgebra('E_' , expression = solve(I - BE)%&%E),
mxAlgebra('SPh' , expression = A_ + C_ + E_),
mxAlgebra('Smz_', expression = rbind(
cbind(SPh, A_+C_),
cbind(A_+C_, SPh))),
mxAlgebra('Sdz', expression=rbind(
cbind(SPh, .5%x%A_+C_),
cbind(.5%x%A_+C_, SPh))),
mxAlgebra('Smz', expression= F%&%Smz_)
)
MZ = mxModel("MZ",
mxExpectationNormal(covariance = "top.Smz",means = "top.mzmu", vnames[1:5]),
mxFitFunctionML()
)
DZ = mxModel("DZ",
mxExpectationNormal(covariance = "top.Sdz",means = "top.dzmu", vnames),
mxFitFunctionML()
)
if (sketch == F){
MZ = mxModel(MZ, mzData)
DZ = mxModel(DZ, dzData)
}
model = mxModel(name, top, MZ, DZ, mxFitFunctionMultigroup(c("MZ","DZ") ))
} else if (length(prss) == 2) {
if (is.null(name)) { name = "MRDoC2"}
top = mxModel("top",
umxMatrix("BE", type = "Full",nrow=4, ncol = 4,byrow = TRUE,
labels = c(NA, "g2", "b1", "b4",
"g1", NA, "b2", "b3",
NA, NA, NA, NA,
NA, NA, NA, NA),
free = c(FALSE, TRUE, TRUE, FALSE,
TRUE, FALSE, FALSE, TRUE,
FALSE, FALSE, FALSE, FALSE,
FALSE, FALSE, FALSE, FALSE)),
umxMatrix('I', type='Iden', nrow= 4,ncol= 4 ),
umxMatrix('F', type='Full', nrow=6, ncol=8, free=FALSE,byrow = TRUE,
values=c(1,0,0,0,0,0,0,0,
0,1,0,0,0,0,0,0,
0,0,1,0,0,0,0,0,
0,0,0,1,0,0,0,0,
0,0,0,0,1,0,0,0,
0,0,0,0,0,1,0,0)),
umxMatrix('LY', type='Full',nrow=4, ncol = 4, free = FALSE, values = diag(4), labels = NA),
umxMatrix('A', type='Symm', nrow=4,ncol = 4,byrow = TRUE,
labels=c("ab2","abraas",NA,NA,
"abraas","as2",NA,NA,
NA,NA,"x2" ,"xrfy",
NA,NA,"xrfy","y2"),
free=c(TRUE,TRUE,FALSE,FALSE,
TRUE,TRUE,FALSE,FALSE,
FALSE,FALSE,TRUE,TRUE,
FALSE,FALSE,TRUE,TRUE)),
umxMatrix('C', type='Symm',nrow=4, ncol = 4,byrow = TRUE,
labels =c("cb2","cbrccs",NA,NA,
"cbrccs","cs2" ,NA,NA,
NA,NA,NA,NA,
NA,NA,NA,NA),
free=c(TRUE,TRUE,FALSE,FALSE,
TRUE,TRUE,FALSE,FALSE,
FALSE,FALSE,FALSE,FALSE,
FALSE,FALSE,FALSE,FALSE)),
umxMatrix('E', type='Symm', nrow=4, ncol = 4,byrow = TRUE,
labels =c("eb2","ebrees",NA,NA,
"ebrees"," es2" ,NA,NA,
NA,NA,NA,NA,
NA,NA,NA,NA),
free= c(TRUE,TRUE,FALSE,FALSE,
TRUE,TRUE,FALSE,FALSE,
FALSE,FALSE,FALSE,FALSE,
FALSE,FALSE,FALSE,FALSE)),
umxMatrix('dzmu', type='Full', nrow=1, ncol=8, free=TRUE, byrow = TRUE,
values = 0, labels=c('muPh1','muPh2','muPS1','muPS2','muPh1','muPh2','muPS1','muPS2')
),
mxAlgebra('mzmu', expression = dzmu %*% t(F)),
mxAlgebra('A_' , expression = solve(I - BE) %&% A),
mxAlgebra('C_' , expression = solve(I - BE) %&% C),
mxAlgebra('E_' , expression = solve(I - BE) %&% E),
mxAlgebra('SPh' , expression= A_ + C_ + E_),
mxAlgebra('Smz_', expression=rbind(
cbind(SPh,A_+C_),
cbind(A_+C_,SPh))),
mxAlgebra('Sdz', expression=rbind(
cbind(SPh,.5%x%A_+C_),
cbind(.5%x%A_+C_,SPh))),
mxAlgebra('Smz', expression= F%&%Smz_)
)
MZ = mxModel("MZ", mzData,
mxExpectationNormal(covariance = "top.Smz",means = "top.mzmu", vnames[1:6]),
mxFitFunctionML()
)
DZ = mxModel("DZ", dzData,
mxExpectationNormal(covariance = "top.Sdz",means = "top.dzmu", vnames),
mxFitFunctionML()
)
if (sketch == F){
MZ = mxModel(MZ, mzData)
DZ = mxModel(DZ, dzData)
}
model = mxModel(name, top, MZ, DZ, mxFitFunctionMultigroup(c("MZ","DZ") ))
} else {
stop("Only 1 or 2 PRSs are supported")
}
model = as(model, "MxModelMRDoC") # set class so that S3s dispatch e.g. plot()
if (sketch == F){
model = mxAutoStart(model)
model = xmu_safe_run_summary(model, autoRun = TRUE, summary = summary, tryHard = tryHard)
}
return(model)
}
#' Present the results of a Mendelian Randomization Direction of Causation Model in a table
#'
#' Summarizes a MR Direction of Causation model, as returned by [umxMRDoC()]
#'
#' @aliases umxSummary.MxModelMRDoC
#' @param model A fitted [umxDoC()] model to summarize
#' @param digits Round to how many digits (default = 2)
#' @param std Whether to show the standardized model (TRUE) (ignored: used extended = TRUE to get unstandardized)
#' @param CIs Confidence intervals (default FALSE)
#' @param comparison Run mxCompare on a comparison model (default NULL)
#' @param RMSEA_CI Optionally compute CI on RMSEA.
#' @param report Print tables to the console (as 'markdown'), or open in browser ('html')
#' @param file The name of the dot file to write: NA = none; "name" = use the name of the model
#' @param ... Optional additional parameters
#' @return - nothing
#' @export
#' @family Summary functions
#' @seealso - \code{\link{umxDoC}()}, [plot()], [umxSummary()] work for DoC models.
#' @md
umxSummaryMRDoC <- function(model, digits = 2, std = TRUE, CIs = FALSE, comparison = NULL, RMSEA_CI = FALSE, report = c("markdown", "html"), file = getOption("umx_auto_plot"), ...) {
if (CIs){ model = umxCI(model, run = "if necessary")}
model_summary = summary(model, refModels = mxRefModels(model, run = TRUE))
if (CIs == TRUE) {
tmp = cbind(model_summary$parameters, model_summary$CI)
# print(tmp)
tmp = tmp[, c(1,12, 6, 11,13)]
} else {
tmp = model_summary$parameters
tmp[, c(1, 5:8)]
}#
umx_print(tmp, digits = digits, report =report, caption = paste0("Parameter loadings for model ", omxQuotes(model$name)), na.print = "", zero.print = "0", justify = "none")
with(model_summary, {
if(!is.finite(TLI)){
TLI_OK = "OK"
} else {
if(TLI > .95) {
TLI_OK = "OK"
} else {
TLI_OK = "bad"
}
}
if(!is.finite(RMSEA)) {
RMSEA_OK = "OK"
} else {
if(RMSEA < .06){
RMSEA_OK = "OK"
} else {
RMSEA_OK = "bad"
}
}
if(RMSEA_CI){
RMSEA_CI = RMSEA(model_summary)$txt
} else {
RMSEA_CI = paste0("RMSEA = ", round(RMSEA, 3))
}
Chi = model_summary$fit
ChiDoF = model_summary$degreesOfFreedom
fitMsg = paste0("\nModel Fit: \u03C7\u00B2(", ChiDoF, ") = ", round(Chi, 2), # was A7
# "Chi2(", ChiDoF, ") = ", round(Chi, 2), # was A7
", p " , umx_APA_pval(p, .001, 3, addComparison = TRUE),
"; CFI = " , round(CFI, 3),
"; TLI = " , round(TLI, 3),
"; ", RMSEA_CI
)
message(fitMsg)
if(TLI_OK != "OK"){ message("TLI is worse than desired (>.95)") }
if(RMSEA_OK != "OK"){ message("RMSEA is worse than desired (<.06)")}
})
}
#' @export
umxSummary.MxModelMRDoC <- umxSummaryMRDoC
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.