孟德爾隨機化(Mendelian Randomization) 統計功效（power）和樣本量計算
孟德爾隨機化(Mendelian Randomization) 統計功效（power）和樣本量計算
1 統計功效（power）概念
統計功效（power）指的是在原假設為假的情況下，接受備擇假設的概率。
用通俗的話說就是，P<0.05時，結果顯著（接受備擇假設）; 在此結論下，我們有多大的把握堅信結果的顯著性，此時需要用到power來表示這種“把握”。
統計功效（power）的計算公式為 1-β。
說到β，要提一下假設檢驗中的一型錯誤和二型錯誤。
一型錯誤，用 α 表示，全稱 Type-I error；
二型錯誤，用 β 表示，全稱 type-II error；
有個比較經典的圖表示 Type-I error 和 type-II error：
(圖片來源忘了，侵刪)
因此，Power越大，犯第二型錯誤的概率越小，我們就更有把握認為結果是顯著的。
下面分別從網頁版和程式碼版講一下怎麼計算power和樣本量，網頁版和程式碼版均可完成分析，任選其一。
2 網頁版計算孟德爾隨機化power和樣本量
網頁版的見地址：https://shiny.cnsgenomics.com/mRnd/
2.1 網頁版計算孟德爾隨機化power
計算power需要用到7個輸入引數，分別為sample size, α, βyx, βOLS, R2xz, σ2(x), σ2(y)。 見下圖：
第一個引數sample size，指的是研究的樣本量大小；在這裡假定樣本量是1000；
第二個引數是 α，指的是一型錯誤(Type-I error)，預設0.05；
第三個引數是βyx，指的是暴露變數和結局變數之間 真實 的相關係數。如何理解 真實 呢，以大胸和不愛運動為例，在校正了性別和年齡等一系列可能會影響大胸和不愛運動的變數後得到的迴歸係數，稱為暴露變數(不愛運動)和結局變數（大胸）真實的相關係數；
第四個引數是βOLS，指的是暴露變數(不愛運動)和結局變數（大胸）之間 觀察到 的相關係數，跟βyx的區別在於，這裡不校正協變數；
第五個引數是R2xz，指的是工具變數（一般指SNP）對暴露變數(不愛運動)的解釋度；
第六個引數是σ2(x)，指的是暴露變數(不愛運動)的方差；
第七個引數是σ2(y)，指的是結局變數（大胸）的方差；
有了這7個引數以後，我們就可以計算power了。 power結果如下所示：
2.2 網頁版計算孟德爾隨機化樣本量
這個步驟同計算power的步驟，唯一不同的是，這個步驟是通過給定power，計算該power下需要的樣本量；
在這裡，我們給定的power是0.8，其他的引數同上面的步驟，得到的樣本量如下所示：
3 程式碼版計算孟德爾隨機化power和樣本量
該程式碼出自網站：https://github.com/kn3in/mRnd
3.1 程式碼版計算孟德爾隨機化power
在Rscript中執行results函式（以下程式碼完全照搬，不要修改任何引數）：
results <- function(N, alpha, byx, bOLS, R2xz, varx, vary, epower) {
threschi <- qchisq(1 - alpha, 1) # threshold chi(1) scale
f.value <- 1 + N * R2xz / (1 - R2xz) #R2xz, Proportion of variance explained for the association between the SNP or allele score (Z) and the exposure variable (X)
con <- (bOLS - byx) * varx # covariance due to YX confounding
vey <- vary - byx * varx * (2 * bOLS - byx)
if (vey < 0) {
data.frame(Error = "Error: Invalid input. The provided parameters result in a negative estimate for variance of the error term in the two-stage least squares model.")
} else {
if (is.na(epower)) {
b2sls <- byx + con / (N * R2xz)
v2sls <- vey / (N * R2xz * varx)
NCP <- b2sls^2 / v2sls
# 2-sided test
power <- 1 - pchisq(threschi, 1, NCP)
data.frame(Parameter = c("Power", "NCP", "F-statistic"), Value = c(power, NCP, f.value), Description = c("", "Non-Centrality-Parameter", "The strength of the instrument"))
} else {
# Calculation of sample size given power
z1 <- qnorm(1 - alpha / 2)
z2 <- qnorm(epower)
Z <- (z1 + z2)^2
# Solve quadratic equation in N
a <- (byx * R2xz)^2
b <- R2xz * (2 * byx * con - Z * vey / varx)
c <- con^2
N1 <- ceiling((-b + sqrt(b^2 - 4 * a * c)) / (2 * a)) #ceiling返回對應數字的'天花板'值，就是不小於該數字的最小整數
data.frame(Parameter = "Sample Size", Value = N1)
}
}
}
隨後執行以下如下命令：
results(N=1000,alpha=0.05, byx=1.3, bOLS=1.41, R2xz=0.01, varx=1, vary=116.6, epower=NA)
各個引數代表的意義如下所示：
alpha=0.05 #Type-I error rate
N=1000 # Sample size
byx=1.3 #the regression coefficients for the association between exposure (X) and outcome (Y) variables (adjusted for confounders).
R2xz=0.01 # genetic instrument that explains R2xz=0.01 of variation in exposure (X)
bOLS=1.41 # the regression coefficients for the association between exposure (X) and outcome (Y) variables (no confounder-adjustment)
varx=1 # Variance of the exposure variable (X)
vary=116.6 #Variance of the outcome variable (Y)
得到的結果如下所示：
3.2 程式碼版計算孟德爾隨機化樣本量
該步驟與前面一致，執行results函式後，再執行如下命令：
results(N=NA,alpha=0.05, byx=1.3, bOLS=1.41, R2xz=0.01, varx=1, vary=116.6, epower=0.8)
各個引數代表的意義如下所示：
alpha=0.05 #Type-I error rate
epower=0.8 # 1-(type-II error rate)
byx=1.3 #the regression coefficients for the association between exposure (X) and outcome (Y) variables (adjusted for confounders).
R2xz=0.01 # genetic instrument that explains R2xz=0.01 of variation in exposure (X)
bOLS=1.41 # the regression coefficients for the association between exposure (X) and outcome (Y) variables (no confounder-adjustment)
varx=1 # Variance of the exposure variable (X)
vary=116.6 #Variance of the outcome variable (Y)
得到的結果如下所示：
原文出處：Brion M J A, Shakhbazov K, Visscher P M. Calculating statistical power in Mendelian randomization studies[J]. International journal of epidemiology, 2013, 42(5): 1497-1501.
此推文感謝彭師姐推薦~