有何AI与医学：孟德尔随机化准备从大家问题多地方以及简单介绍写起来，然后补全整个分析流程！欢迎转发关注哦！
之前的文章介绍了非线性孟德随机化的基本原理，以及多基因风险得分的计算（GRS），今天来介绍使用nlmr-R包完成非线孟德尔随机化分析。
nlmr包实践
1#library(devtools)
2#install_github("jrs95/nlmr")
3library(nlmr)
4### IV (g), exposure (x) & outcome (y)
5epsx <- rexp(10000)#生成10000个均值为1的指数分布随机数
6u <- runif(10000, 0, 1)#最小0最大1
7g <- rbinom(10000, 2, 0.3)#0 1 2对应三个基因型
8epsy <- rnorm(10000)#均值为0方差为1
9ag <- 0.25
10x <- 1 + ag * g + u + epsx#这里可以找最行的GRS算
11y <- 0.15 * x^2 + 0.8 * u + epsy#这里直接是我们对应的结果,疾病就二分类就行，family对应binomial
12
13### Covariates (c)
14c1 = rnorm(10000)#模拟连续协变量，例如年龄
15c2 = rnorm(10000)#同上
16c3 = rbinom(10000, 2, 0.33)#模拟离散协变量，例如吸烟
17c = data.frame(c1 = c1, c2 = c2, c3 = as.factor(c3))
18
19### Analyses
20#通过使用回归将分数多项式模型拟合到局部平均因果效应来进行工具变量分析，光滑曲线。
21#y结局x暴露g基因型c协变量举证，family方法，q分段数，d拟合次数次1或者2默认1不重要，ci置信区间
22fp = fracpoly_mr(y, x, g, c, family ="gaussian", q = 10, d ="both", ci ="model_se", fig = TRUE)
23summary(fp)
24#分段节点式
25plm = piecewise_mr(y, x, g, c, family ="gaussian", q = 10, nboot = 50, fig = TRUE)
26summary(plm)
markdown结果
1#library(devtools)
2#install_github("jrs95/nlmr")
3library(nlmr)
4### IV (g), exposure (x) & outcome (y)
5epsx <- rexp(10000)
6u <- runif(10000, 0, 1)
7g <- rbinom(10000, 2, 0.3)
8epsy <- rnorm(10000)
9ag <- 0.25
10x <- 1 + ag * g + u + epsx
11y <- 0.15 * x^2 + 0.8 * u + epsy
12
13#
## Covariates (c)
14c1 = rnorm(10000)
15c2 = rnorm(10000)
16c3 = rbinom(10000, 2, 0.33)
17c = data.frame(c1 = c1, c2 = c2, c3 = as.factor(c3))
18
19#
## Analyses
20fp = fracpoly_mr(y, x, g, c, family = "gaussian", q = 10, d = "both", ci = "model_se", fig = TRUE)
21summary(fp)
22## Call: fracpoly_mr
23##
24## Number of individuals: 10000; Quantiles: 10; 95%CI: Model based SEs
25##
26## Powers: 2
27##
28## Coefficients:
29## Estimate Std. Error 95%CI Lower 95%CI Upper p.value
30## 2 0.144987 0.012498 0.120491 0.1695 < 2.2e-16 ***
31## ---
32## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
33##
34## Non-linearity tests
35## Fractional polynomial degree p-value: 0.924
36## Fractional polynomial non-linearity p-value: 0.000535
37## Quadratic p-value: 0.000532
38## Cochran Q p-value: 0.103
39##
40## Heterogeneity tests
41## Cochran Q p-value: 0.937
42## Trend p-value: 0.215
43
44plm = piecewise_mr(y, x, g, c, family = "gaussian", q = 10, nboot = 50, fig = TRUE)
45summary(plm)
46## Call: piecewise_mr
47##
48## Number of individuals: 10000; Quantiles: 10; Number of bootstrap replications: 50
49##
50## LACE:
51## Estimate Std. Error 95%CI Lower 95%CI Upper p.value
52## 1 0.370222 0.194433 -0.010866 0.7513 0.0568951 .
53## 2 0.324060 0.194220 -0.056611 0.7047 0.0952126 .
54## 3 0.728189 0.196023 0.343984 1.1124 0.0002034 ***
55## 4 0.674931 0.191170 0.300237 1.0496 0.0004147 ***
56## 5 0.767218 0.207506 0.360507 1.1739 0.0002179 ***
57## 6 0.766719 0.191298 0.391775 1.1417 6.124e-05 ***
58## 7 0.673263 0.190727 0.299438 1.0471 0.0004156 ***
59## 8 0.888327 0.199124 0.498044 1.2786 8.151e-06 ***
60## 9 1.096948 0.196101 0.712590 1.4813 2.222e-08 ***
61## 10 1.389435 0.416601 0.572897 2.2060 0.0008525 ***
62## ---
63## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
64##
65## Non-linearity tests
66## Quadratic p-value: 0.000532
67## Cochran Q p-value: 0.103
68##
69## Heterogeneity tests
70## Cochran Q p-value: 0.937
71## Trend p-value: 0.215
