R資料分析:樣本量計算的底層邏輯與實操,pwr包

Codewar發表於2022-03-03

樣本量問題真的是好多人的老大難,是很多同學科研入門第一個攔路虎,今天給本科同學改大創標書又遇到這個問題,我想想不止是本科生對這個問題不會,很多同學從上研究生到最後脫離科研估計也沒能把這個問題弄得很明白,那麼希望大夥兒在看了這篇文章能夠更加深入地理解樣本量計算的邏輯,也能對大家的科研設計中的樣本量設計部分有所啟發。

樣本量計算的邏輯

還記得我們最開始接觸統計推斷的時候,大家都知道一個詞叫做原假設,原假設一般來講都是“陰性的”,我們統計推斷要做的事情便是推翻原假設從而得出有“統計學意義的結果”,怎麼去推翻?就是在一次隨機事件中小概率的我們認為不可能發生的符合備擇假設的事件發生了,“不可能發生的事情竟然在現實中隨便一次抽樣研究中竟然就發生了”,就反證出原假設站不住腳,所以我們的結果就有“統計學意義”就得到了所謂的陽性結果了。

上面的這段話和樣本量又如何扯上關係呢?繼續往下看:

舉個例子,比如現在兩個研究團隊都正在做一個性別比例的研究,兩個團隊都想通過自己的研究去證明世界上男性數量和女性數量不相等。(例子純屬虛構方便大家理解)

那麼原假設就是世界上男性佔比等於女性佔比,現在兩個研究團隊都期望能用他們的抽樣樣本去推翻原假設,得到男性數量不等於女性數量的結論。

A團隊為了省事將樣本量定為了4個,那麼A團隊做一次研究最極端(可能性最小)的情況是4個全是同一個性別,那麼此時這種最極端的情況發生的概率也為2*0.5^5=0.0625,其實也還不算我們通常認為的小概率(0.05),所以在4個樣本量的情況下,A團隊無論如何進行精妙的研究設計,並且A團隊把這個設計精妙的研究無論重複多少次,都不可能得到陽性的研究結果,即使真實世界的男女比例真的不一樣,A團隊也100%發現不了陽性結果,就是因為樣本量太小導致的。

大家好好體會下上面的例子,然後再看B團隊的做法:

B團隊就很雞賊了,將樣本量定為5個,此時最極端的情況是5個全是同一個性別,那麼此時最極端的情況發生的概率為2*0.5^6=0.03125,哇,有點意思了,就是這個小概率事件理論上是可以發生的,而且,B團隊只要有耐心,理論上重複100次該研究設計,大約會有3次都可以拒絕原假設,所以B團隊如果在一次實驗中剛好選的樣本中5個都是同一性別,那麼B團隊就成功得到了陽性結果,但是仔細想想,如果真實世界中真的男女比例不一樣,B團隊因為樣本量只有5個,所以僅僅只能寄希望5個樣本全部是同性才能得到想要的結果,而且這種情況100次才大概能出現三次,是不是這個發現陽性結果的能力有點小呀?

就是即使真實世界中男女佔比不一樣,B團隊也1-3%=97%的概率發現不了。

大家認真閱讀上面我虛構的例子,你就能體會到,樣本量估計的極端重要性以及幾個可以影響樣本量大小的指標:

分別是:

小概率事件的定義α,也就是顯著性水平,我們的樣本統計量達到了這個顯著性水平就可以拒絕原假設(不管原假設是不是真)從而得到陽性結果,那麼這就存在一種可能性就是原假設是真的話因為達到了顯著性水平,你也把人家給拒絕了,這個就叫做一類錯誤;就是在前面的例子中,其實真實世界中男女比例是相等的,但是你就偏偏真的抽了5個全是同性的樣本得到了男女比例不同的錯誤結論,如果這個時候你將α調低一點,這樣的錯誤就會減少,比如調到0.3以下,那麼在例子中B團隊就不會拒絕原假設了,它需要更多的樣本量才有可能拒絕原假設,同樣的你也可以將α調高一點,這樣A團隊在理論上也可以得到男女數量不等的陽性結果。

我們知道了α是小概率事件的概率,相應地1-α為confidence level,就是就是假設做了無數次研究,每次研究都有相應的結局推斷,那麼1-α%的研究得到推斷我們有信心認為一定是正確反應了真實情況的,或者說1-α%的研究樣本估計值的某個區間一定是包含總體的真實值的。

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

統計效能1-β,也就前面例子中的能力,就是如果原假設真的為假,我可以正確拒絕原假設的能力,比如對於團隊B來講5個樣本量的時候,理論上如果原假設為假,就是男女數量真的是不同的,5個樣本量是可以得到拒絕原假設的結果的,但是隻有在5人全為同性別的時候才能得到(只有2中情況,5個樣本總的情況有2^5種,所以統計效能只有2/2^5,非常小),就是對於團隊B來講雖然理論上有能力得到陽性結果,但這個能力非常小,這個就叫做統計效能不足,我們一般期望統計效能能達到0.8,所以說團隊B還是應該繼續增加樣本量。

統計效能1-β是如果原假設真的為假,我可以正確拒絕原假設的概率,那麼自然你就理解了β就是原假設真的為假,卻沒能拒絕原假設的概率,叫做二類錯誤。

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

當然α和β不是本文的重點,但是希望看了上面的文字大家可以明白α和β為什麼能影響樣本量,我們繼續思考,如果現在我們的研究變成了“是不是世界上的男(女)性比女(男)性多100%?”,此時是不是我的備擇假設就成了“世界上性別構成比之差大於33.3%”,此時B團隊5個樣本理論上也不再能得到陽性的結果了。就是這個33.3%這個數也能影響樣本量,這個叫做效應值d:

The size of the effect that is worthwhile to detect (d)

到這兒,推斷統計的樣本計算的邏輯基本上就給大家寫完了(還有一個是總體的變異SD or σ,也會影響樣本量,但在本例中不適用),本文的重點是教會大家如何計算統計推斷的樣本量和統計描述的樣本量。

樣本量計算的不同情形

我們先看統計推斷的樣本量估計,分兩種情況,一種是以連續變數為結局,一種是以分類變數為結局,當結局是連續變數時,並且只有兩組比較時(常見的RCT研究樣本量計算都可以參照這個公式),有計算公式如下:

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

Gogtay N. J. (2010). Principles of sample size calculation. Indian journal of ophthalmology, 58(6), 517–518. https://doi.org/10.4103/0301-4738.71692

比如在下面的例子中應該是要85個樣本(84.6)每組的:

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

R資料分析:樣本量計算的底層邏輯與實操,pwr包

其實是84.6,算下來應該是85

上面的例子中給了標準差,但是這兒要特別注意d的計算,d的計算方法是均值差除以合併標準差:

This effect size is equal to the difference between the means at the endpoint, divided by the pooled standard deviation

但是合併標準差我們一般是不知道的,這個時候可以去找文獻,去預試驗,去大資料中挖掘或者直接極差/4估計一下,都是可以。

我們再看一個例子:

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

上面的例子是用PASS軟體計算出來的樣本量。依然是利用了了α,β,d和σ,得到的結果是85.

當結局是分類變數時,我們的計算公式如下:

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

參考文獻如下:

Wang, H. and Chow, S.-C. 2007. Sample Size Calculation for Comparing Proportions. Wiley Encyclopedia of Clinical Trials.

依然是給大家一個例子:

一項研究欲探討視訊觀察治療(VOT)對結核病治療是否優於直接監督治療(DOT),隨機化後2個月內是否完成80%或更多的預定治療觀察作為主要指標。預計VOT指標完成率為75%,DOT指標完成率為60%,檢驗效能為90%,檢驗水準為0.05,雙側檢驗,組間比例為1:1。試用PASS軟體估算所需最小樣本量。

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

上面的例子中,用PASS軟體得到的樣本量應該為200,如下圖:

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

可以注意到,這個時候計算樣本量用到了α,β,和d,就並沒有σ,因為分類結局本身就不存在σ嘛。

樣本量計算的R語言實操

剛剛給大家寫了樣本量的理論知識,以及在RCT中結局是連續變數和分類變數時樣本量計算的例項。這部分教大家在R語言中進行相應的計算,R中有很多包可以幫助我們計算樣本量,今天介紹pwr包。

比如我們的結局是連續變數時,連續變數部分給大家寫的前面兩個例子的樣本量計算程式碼都如下:

pwr.t.test(n = NULL,
                sig.level = 0.05, 
                type = "two.sample", 
                alternative = "greater", 
                power = 0.90, 
                d = 0.5)

執行上面的程式碼即可得到樣本量n為85和文獻中和PASS軟體的結果都是一樣的:

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

當然啦,如果大家的實驗是非劣效性試驗或者優效性實驗,只需要將alternaive引數改改就可以了。如果大家的實驗兩組不等,比如只定了對照組的人數,想要確定實驗組多少人也是可以用pwr包的pwr.t2n.test計算的。

再看結局是分類變數時的樣本量計算方法,這個時候要用到的函式是pwr.2p.test。比如上一部分例子中樣本量PASS算出來是200,同樣的我們用R語言寫出如下程式碼:

pwr.2p.test(h = ES.h(p1 = 0.75, p2 = 0.60),
           n = NULL,
           sig.level = 0.05,
           power = 0.90,
           alternative = "two.sided")

執行後我們可以看到,結果也是202(結果的些微不同是因為一個軟體用的是normal approximation method,一個用的arcsine transformation,都不算錯):

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

當然啦,如果大家的實驗是非劣效性試驗或者優效性實驗,只需要將alternaive引數改改就可以了。如果大家的實驗兩組不等,比如只定了對照組的人數,想要確定實驗組多少人也是可以用pwr包的pwr.2p2n.test計算的。

上面的整個邏輯是統計推斷的樣本量計算,但是整體看,我們樣本量推斷是有自己的流程的,如下:

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

第一步就是弄明白你的研究變數,比如你的研究就是一個現況調查,這個時候你關心的是樣本量是不是能夠正確地估計出總體引數,比如估計出總體均值,總體率的估計等等;比如你的研究是隨機對照實驗,你關心的是如果幹預真的有效,你現有的樣本量是不是有很大的把握識別出來有統計學意義的結果。

所以計算樣本量之前你首先要弄明白你的研究要搞啥,不同的研究型別決定了樣本量計算方法的不同:

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

其實做結果方程,做因子分析,做預測模型考慮的樣本量東西又不一樣了,反正大家要知道研究設計不同資料型別不同考慮的樣本量又不一樣了,但本文不展開。

R資料分析:樣本量計算的底層邏輯與實操,pwr包

 

本文只給大家寫了RCT的樣本量計算。

因為篇幅所限,描述統計的樣本量,觀察性研究的樣本量計算就下回再寫。

最後再給大家貼兩篇文獻,感興趣的同學可以去瞅瞅:

Wang, X. and Ji, X., 2020. Sample size estimation in clinical research: from randomized controlled trials to observational studies. Chest, 158(1), pp.S12-S20.

Wang, X. and Ji, X., 2020. Sample size formulas for different study designs: supplement document for sample size estimation in clinical research.

小結

今天給大家寫了統計推斷樣本量計算的原理和pwr包如何計算樣本量,還有個叫Openepi的軟體也可以計算樣本量,PASS也可以計算樣本量,而且這些都是免費的,大家可以擇一使用。感謝大家耐心看完,自己的文章都寫的很細,重要程式碼都在原文中,希望大家都可以自己做一做,請轉發本文到朋友圈後私信回覆“資料連結”獲取所有資料和本人收集的學習資料。如果對您有用請先記得收藏,再點贊分享。

也歡迎大家的意見和建議,大家想了解什麼統計方法都可以在文章下留言,說不定我看見了就會給你寫教程哦,有疑問歡迎私信,有合作意向請莫要猶豫直接滴滴我。

相關文章