孟德爾隨機化(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和樣本量
網頁版的見地址：//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和樣本量
該代碼出自網站：//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.
此推文感謝彭師姐推薦~
