Standard errors of MR GxE

Author

Gibran Hemani

Published

December 1, 2023

Background

Coding the MR-GxE model to have standard errors obtained through parametric bootstrap

library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
est <- function(b_gx, se_gx, b_gy, se_gy, nboot=1000) {
    npop <- length(b_gx)
    stopifnot(length(se_gx) == npop)
    stopifnot(length(b_gy) == npop)
    stopifnot(length(se_gy) == npop)
    
    mod <- summary(lm(b_gy ~ b_gx))

    # standard errors
    o <- lapply(1:nboot, \(i) {
        bgxb <- rnorm(npop, b_gx, se_gx)
        bgyb <- rnorm(npop, b_gy, se_gy)
        modb <- summary(lm(bgyb ~ bgxb))$coef
        tibble(boot=i, a=modb[1,1], b=modb[2,1])
    }) %>% bind_rows()

    res <- tibble(
        a = mod$coef[1,1],
        b = mod$coef[2,1],
        a_se = sd(o$a),
        b_se = sd(o$b),
        a_pval = pnorm(abs(a) / a_se, lower.tail=FALSE),
        b_pval = pnorm(abs(b) / b_se, lower.tail=FALSE),
        a_mean = mean(o$a),
        b_mean = mean(o$b)
    ) %>% as.list()
    return(res)
}

o <- est(
    b_gx=c(0, 1),
    se_gx = c(0.2, 0.2),
    b_gy=c(1, 1),
    se_gy = c(0.2, 0.2)
)
o
$a
[1] 1

$b
[1] 2.507433e-16

$a_se
[1] 0.227826

$b_se
[1] 0.3302168

$a_pval
[1] 5.685413e-06

$b_pval
[1] 0.5

$a_mean
[1] 1.004683

$b_mean
[1] -0.00958491
sessionInfo()
R version 4.3.2 (2023-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

time zone: Europe/London
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] dplyr_1.1.4

loaded via a namespace (and not attached):
 [1] digest_0.6.33     utf8_1.2.4        R6_2.5.1          fastmap_1.1.1    
 [5] tidyselect_1.2.0  xfun_0.41         magrittr_2.0.3    glue_1.6.2       
 [9] tibble_3.2.1      knitr_1.45        pkgconfig_2.0.3   htmltools_0.5.7  
[13] rmarkdown_2.25    generics_0.1.3    lifecycle_1.0.4   cli_3.6.1        
[17] fansi_1.0.5       vctrs_0.6.4       compiler_4.3.2    tools_4.3.2      
[21] pillar_1.9.0      evaluate_0.23     yaml_2.3.7        rlang_1.1.2      
[25] jsonlite_1.8.7    htmlwidgets_1.6.3