|
| 1 | +#' MR-GRIP: a modified MR-Egger model with the Genotype Recoding Invariant Property |
| 2 | +#' |
| 3 | +#' This implements the modified MR-Egger model with the Genotype Recoding Invariant Property (MR-GRIP) due to Dudbridge and Bowden et al. (2025). |
| 4 | +#' It is well known that the results of MR-Egger are sensitive to which alleles are designated as the effect alleles. |
| 5 | +#' A pragmatic convention is to orient all SNPs to have positive effects on the exposure, which has some advantages in interpretation but also brings some philosophical limitations. |
| 6 | +#' The MR-GRIP model is a modification to the MR-Egger model in which each term is multiplied by the genotype-phenotype associations. |
| 7 | +#' This makes each term in the model invariant to allele coding. |
| 8 | +#' |
| 9 | +#' @param b_exp Vector of genetic effects on exposure. |
| 10 | +#' @param b_out Vector of genetic effects on outcome. |
| 11 | +#' @param se_exp Standard errors of genetic effects on exposure. |
| 12 | +#' @param se_out Standard errors of genetic effects on outcome. |
| 13 | +#' @param parameters List of parameters. |
| 14 | +#' |
| 15 | +#' @export |
| 16 | +#' @return List of with the following elements: |
| 17 | +#' \describe{ |
| 18 | +#' \item{b}{MR estimate} |
| 19 | +#' \item{se}{Standard error of MR estimate} |
| 20 | +#' \item{pval}{p-value of MR estimate} |
| 21 | +#' \item{Q, Q_df, Q_pval}{Heterogeneity stats} |
| 22 | +#' \item{b.wi}{MR estimate adjusting for weak instruments} |
| 23 | +#' \item{se.wi}{Standard error adjusting for weak instruments} |
| 24 | +#' \item{pval.wi}{p-value adjusting for weak instruments} |
| 25 | +#' \item{mod}{Summary of regression} |
| 26 | +#' \item{dat}{Original data used for MR-GRIP} |
| 27 | +#' } |
| 28 | +mr_grip <- function(b_exp, b_out, se_exp, se_out, parameters) { |
| 29 | + if (length(b_exp) != length(b_out)) stop("The lengths of b_exp and b_out are not equal.") |
| 30 | + if (length(se_exp) != length(se_out)) stop("The lengths of se_exp and se_out are not equal.") |
| 31 | + if (length(b_exp) != length(se_out)) stop("The lengths of b_exp and se_out are not equal.") |
| 32 | + |
| 33 | + nulllist <- list( |
| 34 | + b = NA, |
| 35 | + se = NA, |
| 36 | + pval = NA, |
| 37 | + nsnp = NA, |
| 38 | + Q = NA, |
| 39 | + Q_df = NA, |
| 40 | + Q_pval = NA, |
| 41 | + mod = NA, |
| 42 | + smod = NA, |
| 43 | + dat = NA |
| 44 | + ) |
| 45 | + if ( |
| 46 | + sum(!is.na(b_exp) & !is.na(b_out) & !is.na(se_exp) & !is.na(se_out)) < 3 |
| 47 | + ) { |
| 48 | + return(nulllist) |
| 49 | + } |
| 50 | + |
| 51 | + dat <- data.frame( |
| 52 | + b_out = b_out, |
| 53 | + b_exp = b_exp, |
| 54 | + se_exp = se_exp, |
| 55 | + se_out = se_out |
| 56 | + ) |
| 57 | + grip_out <- b_out * b_exp |
| 58 | + grip_exp <- b_exp^2 |
| 59 | + # GRIP regression. Includes intercept. Weights designed to replicate IVW under no intercept. |
| 60 | + mod <- stats::lm(grip_out ~ grip_exp, weights = 1 / (grip_exp * se_out^2)) |
| 61 | + smod <- summary(mod) |
| 62 | + b <- stats::coefficients(smod)[2, 1] |
| 63 | + se <- stats::coefficients(smod)[2, 2] |
| 64 | + b.adj <- NA |
| 65 | + se.adj <- NA |
| 66 | + pval.adj <- NA |
| 67 | + pval <- 2 * stats::pt(abs(b / se), length(b_exp) - 2L, lower.tail = FALSE) |
| 68 | + return(list( |
| 69 | + b = b, |
| 70 | + se = se, |
| 71 | + pval = pval, |
| 72 | + b.adj = b.adj, |
| 73 | + se.adj = se.adj, |
| 74 | + pval.adj = pval.adj, |
| 75 | + nsnp = length(b_exp), |
| 76 | + mod = smod, |
| 77 | + dat = dat |
| 78 | + )) |
| 79 | +} |
0 commit comments