1 統(tǒng)計功效（power）概念
統(tǒng)計功效（power）指的是在原假設(shè)為假的情況下濒旦，接受備擇假設(shè)的概率再登。
用通俗的話說就是锉矢，P<0.05時，結(jié)果顯著（接受備擇假設(shè)）; 在此結(jié)論下灯节，我們有多大的把握堅信結(jié)果的顯著性绵估，此時需要用到power來表示這種“把握”。
統(tǒng)計功效（power）的計算公式為 1-β形入。
說到β唯笙，要提一下假設(shè)檢驗中的一型錯誤和二型錯誤。
一型錯誤，用 α 表示苞慢，全稱 Type-I error英妓；
二型錯誤，用 β 表示辑畦，全稱 type-II error腿倚；
有個比較經(jīng)典的圖表示 Type-I error 和 type-II error：
(圖片來源忘了敷燎，侵刪)
因此硬贯，Power越大，犯第二型錯誤的概率越小鸵赖，我們就更有把握認為結(jié)果是顯著的拄衰。
下面分別從網(wǎng)頁版和代碼版講一下怎么計算power和樣本量肾砂，網(wǎng)頁版和代碼版均可完成分析，任選其一包吝。
2 網(wǎng)頁版計算孟德爾隨機化power和樣本量
網(wǎng)頁版的見地址：https://shiny.cnsgenomics.com/mRnd/
2.1 網(wǎng)頁版計算孟德爾隨機化power
計算power需要用到7個輸入?yún)?shù)诗越，分別為sample size, α, βyx, βOLS, R2xz, σ2(x), σ2(y)息堂。 見下圖：
第一個參數(shù)sample size床未，指的是研究的樣本量大小斋扰；在這里假定樣本量是1000啃洋；
第二個參數(shù)是 α宏娄，指的是一型錯誤(Type-I error)，默認0.05粮宛；
第三個參數(shù)是βyx十饥，指的是暴露變量和結(jié)局變量之間 真實 的相關(guān)系數(shù)逗堵。如何理解 真實 呢，以大胸和不愛運動為例汁咏，在校正了性別和年齡等一系列可能會影響大胸和不愛運動的變量后得到的回歸系數(shù)作媚，稱為暴露變量(不愛運動)和結(jié)局變量（大胸）真實的相關(guān)系數(shù)纸泡；
第四個參數(shù)是βOLS，指的是暴露變量(不愛運動)和結(jié)局變量（大胸）之間 觀察到 的相關(guān)系數(shù)蚤假，跟βyx的區(qū)別在于磷仰，這里不校正協(xié)變量境蔼；
第五個參數(shù)是R2xz伺通，指的是工具變量（一般指SNP）對暴露變量(不愛運動)的解釋度罐监；
第六個參數(shù)是σ2(x)瞒爬，指的是暴露變量(不愛運動)的方差疮鲫；
第七個參數(shù)是σ2(y)弦叶，指的是結(jié)局變量（大胸）的方差伤哺；
有了這7個參數(shù)以后，我們就可以計算power了绢彤。 power結(jié)果如下所示：
2.2 網(wǎng)頁版計算孟德爾隨機化樣本量
這個步驟同計算power的步驟茫舶，唯一不同的是刹淌，這個步驟是通過給定power，計算該power下需要的樣本量有勾；
在這里疹启，我們給定的power是0.8，其他的參數(shù)同上面的步驟蔼卡，得到的樣本量如下所示：
3 代碼版計算孟德爾隨機化power和樣本量
該代碼出自網(wǎng)站：https://github.com/kn3in/mRnd
3.1 代碼版計算孟德爾隨機化power
在Rscript中運行results函數(shù)（以下代碼完全照搬喊崖，不要修改任何參數(shù)）：
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返回對應(yīng)數(shù)字的'天花板'值，就是不小于該數(shù)字的最小整數(shù)
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)
各個參數(shù)代表的意義如下所示：
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)
得到的結(jié)果如下所示：
3.2 代碼版計算孟德爾隨機化樣本量
該步驟與前面一致荤懂，運行results函數(shù)后，再運行如下命令：
results(N=NA,alpha=0.05, byx=1.3, bOLS=1.41, R2xz=0.01, varx=1, vary=116.6, epower=0.8)
各個參數(shù)代表的意義如下所示：
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)
得到的結(jié)果如下所示：
原文出處：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.
此推文感謝彭師姐推薦~
