# PRS vs IVW

statistics
Mendelian randomization
Author

Gibran Hemani

Published

May 8, 2022

How does PRS compare to IVW fixed effects analysis

## Simulation study

``````library(simulateGP)
geno1 <- make_geno(10000, 500, 0.5)
b <- choose_effects(500, 0.3)
x1 <- make_phen(b, geno1)
y1 <- make_phen(0.4, x1)

geno2 <- make_geno(1000, 500, 0.5)
x2 <- make_phen(b, geno2)
y2 <- make_phen(0.4, x2)

bhat <- gwas(x1, geno1)
b_unweighted <- sign(b)``````

### Standard unweighted PRS analysis

``````prs_unweighted <- geno2 %*% b_unweighted
summary(lm(x2 ~ prs_unweighted))``````
``````
Call:
lm(formula = x2 ~ prs_unweighted)

Residuals:
Min      1Q  Median      3Q     Max
-3.2356 -0.5547 -0.0330  0.6087  3.1932

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)    -0.321168   0.035444  -9.061   <2e-16 ***
prs_unweighted  0.027292   0.001791  15.239   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9011 on 998 degrees of freedom
Multiple R-squared:  0.1888,    Adjusted R-squared:  0.1879
F-statistic: 232.2 on 1 and 998 DF,  p-value: < 2.2e-16``````

### Meta analysing per-SNP PRS scores

``library(meta)``
``````Loading 'meta' package (version 6.0-0).
Type 'help(meta)' for a brief overview.
Readers of 'Meta-Analysis with R (Use R!)' should install
older version of 'meta' package: https://tinyurl.com/dt4y5drs``````
``````o <- sapply(1:ncol(geno2), function(i)
{
prs_unweighted <- geno2[,i] * b_unweighted[i]
summary(lm(x2 ~ prs_unweighted))\$coef[2,1:2]
})
metafor::rma(yi=o[1,], sei=o[2,], method="EE")``````
``````
Equal-Effects Model (k = 500)

I^2 (total heterogeneity / total variability):   17.58%
H^2 (total variability / sampling variability):  1.21

Test for Heterogeneity:
Q(df = 499) = 605.4622, p-val = 0.0007

Model Results:

estimate      se     zval    pval   ci.lb   ci.ub
0.0277  0.0020  13.8615  <.0001  0.0238  0.0316  ***

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1``````

### Standard errors with number of SNPs

``````ses_unweighted <- sapply(1:ncol(geno2), function(i)
{
prs_unweighted <- geno2[,1:i, drop=FALSE] %*% b_unweighted[1:i]
summary(lm(x2 ~ prs_unweighted))\$coef[2,2]
})
plot(ses_unweighted)``````