要如何了解 P 值(以T test為例)

參考:Biostatistics for Dummies, P41~48)

對我來說,P值最直觀的解釋是: 在現在的虛無假說為真實的情況下,出現目前的結果的機率。

知道兩個重點後對我來說P值變得很好理解:

P值的P是Probability,可能性的P
計算P值的前提是假設虛無假說(H₀)為真

因為我們目前看到的結果,是從母群體抽樣得到的結果。

　眼前的世界並不是真實,只是真實的一部分。

也因為如此,目前看到的差異, 有可能只是抽樣的誤差, 而P值在告訴我們眼前看到的,其實只是抽樣誤差的機率會是多少。

而P值的用處是,出現這樣的結果的機率在虛無假說為真的情況下實在太低了,只是抽樣所導致的誤差的機率實在太低了,所以我們拒絕了虛無假說。拒絕了虛無假說後我們去接受另一個假說(Alternative Hypothesis, H_A),這個假說通常是研究想要證明的結果。因此,事實上你看到的每一個P值背後都有一個虛無假說, 為了理解這些P值的意義, 請你去了解他們的虛無假說是什麼。

Q:以一個Two sample的T test為例, 他的虛無假說是什麼呢?

單用講的太空泛,我們看一個實際的例子。

下面的例子a是平均約為10的10個樣本, b是平均為20的十個樣本。 Two sample的T test的虛無假說是a 和b 的平均相同。而Alternative hypothesis是a 和b的平均不同。來看看T-test有沒有辦法拒絕這個虛無假說(顯著水準假設為0.05)。

set.seed(123)
a <- rnorm(10,10)
b <- rnorm(10,20)

t.test(a,b, var.equal = T)

## 
##  Two Sample t-test
## 
## data:  a and b
## t = -22.733, df = 18, p-value = 1.044e-14
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -11.070569  -9.197423
## sample estimates:
## mean of x mean of y 
##  10.07463  20.20862

這個是什麼表,好多數字好難理解。

這是我一開始的感覺,但沒關係,先從最熟悉的P-value開始。

計算的結果p-value為10^{-14}, 而我們設定的拒絕虛無假說的P值為0.05,此處的P值小於我們設定的標準,因此虛無假說(兩者平均相同)被拒絕了,而我們去接受Alternative hypothsis: a和b的平均不同。

論文裡面當然不會明述這些虛無假說,以及拒絕假說的過程,最後只留下一句a和b的平均有顯著差異。

為什麼P<0.05

P<0.05就代表我們看到的差異是真的嗎?

我們把標準降低一點, 只要P值低於0.5就拒絕虛無假說, 這樣發表不就輕輕鬆鬆了嗎!! 如果我們把標準提高一點(P < 0.000001),這樣不就能保證結果正確無誤了嗎?

想要理解要如何設定拒絕虛無假說的標準,要先理解什麼是type I 和type II error。

Type I error是明明虛無假說是真的,你卻拒絕她了(你接受了一個不存在的Alternative hypothesis)。

Type II error是明明虛無假說是錯的,你卻無法拒絕她。以剛剛的例子來說,真實的情況是兩者平均不同,你卻無法發現。

Q:如果我們設定P<0.5, 哪種情況容易發生?

Q:如果我們設定P<0.0000000001,哪種情況容易發生?

此處設定的P(\(\alpha\))也就是type I error的機率。

如果我們設定P<0.05可以拒絕虛無假說, 代表我們接受type I error的機率可以有5%。5%真的錯了也只能接受。

下面是一個例子,就算兩個數值平均完全相同, 也有0.05機率會P < 0.05。

# 測試兩組來自平均完全相同(0)母群體的樣本的差異(1000次)
p <- map_dbl(1:1000, 
             ~ t.test(rnorm(10,0), rnorm(10,0), var.equal = T)$p.value)

# p值的分布
qplot(p)

# P值<0.05的比例
mean(p<0.05)

## [1] 0.052

我們就是接受了這個5%的type I error.

Note: P<0.05這個設定是歷史的因素所促成的(約定成俗),我自己並不喜歡這個設定。與其說P<0.05,我比較偏向給出計算出來的P是多少。

Power analysis (檢定力分析)

沒有人喜歡做不出結果,檢定不出差異的實驗。

我們想要盡量減少Type II error的機率。

我們有辦法在試驗之前知道檢定出差異的機率(\(\beta\);檢定力)有多少嗎?

\(\beta\) = 1 - Pr(type II error), 我們希望盡量提高\(\beta\).

決定一個試驗的檢定力(\(\beta\))的因素有:

\(\alpha\) value: 你所設定拒絕虛無假說的P, 以P<0.05為例, \(\alpha\) 則於0.05。另外\(\alpha\) 也等於Pr(type I error)

Q: 提高\(\alpha\), \(\beta\)也會隨之??降低/提高
樣本數
存在的差異和其noise的比例(effect size/variance)

variance會影響P值和檢定力,儘管平均真的不同。等等我提到T test 的P值怎麼計算的時候也會提到,但我們可以看看一個例子.

set.seed(123)
a_noisy <- rnorm(n = 10, mean = 10, sd = 100)
b_noisy <- rnorm(n = 10, mean = 20, sd = 100)

t.test(a_noisy, b_noisy, var.equal = T)

## 
##  Two Sample t-test
## 
## data:  a_noisy and b_noisy
## t = -0.5249, df = 18, p-value = 0.6061
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -117.05694   70.25768
## sample estimates:
## mean of x mean of y 
##  17.46256  40.86220

明明差異相同,為什麼這次沒有辦法找到顯著差異? 因為這次的目標的噪音(variance)太強,我們沒有辦法準確地測量效果量(effect size)。

這樣的情況下,為了檢定出差異，我們會需要更多的樣本(1000):

set.seed(123)
a_noisy_alot <- rnorm(n = 1000, mean = 10, sd = 100)
b_noisy_alot <- rnorm(n = 1000, mean = 20, sd = 100)

t.test(a_noisy_alot, b_noisy_alot, var.equal = T)

## 
##  Two Sample t-test
## 
## data:  a_noisy_alot and b_noisy_alot
## t = -2.8229, df = 1998, p-value = 0.004806
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -21.410635  -3.856842
## sample estimates:
## mean of x mean of y 
##  11.61279  24.24653

這三者和\(\beta\)的關係很複雜,但是可以被計算的。

# 這是第一次的例子: beta = 100%
pwr::pwr.t.test(n = 10, sig.level = 0.05, d = 10)

## 
##      Two-sample t test power calculation 
## 
##               n = 10
##               d = 10
##       sig.level = 0.05
##           power = 1
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

# 這是第二次的例子: beta = 5%
pwr::pwr.t.test(n = 10, sig.level = 0.05, d = 0.1)

## 
##      Two-sample t test power calculation 
## 
##               n = 10
##               d = 0.1
##       sig.level = 0.05
##           power = 0.05516129
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

# 這是第三次的例子: beta = 60%
pwr::pwr.t.test(n = 1000, sig.level = 0.05, d = 0.1)

## 
##      Two-sample t test power calculation 
## 
##               n = 1000
##               d = 0.1
##       sig.level = 0.05
##           power = 0.6083667
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

# 樣本數和檢定力的關係(alpha 為0.05)
ggplot() + 
  xlim(c(1, 100)) + 
  geom_function(fun = ~ pwr::pwr.t.test(n = .x, 
                                        sig.level = 0.05, 
                                        d = 0.5)$power, 
                aes(color = "Effect size = 0.5")) +
  geom_function(fun = ~ pwr::pwr.t.test(n = .x, 
                                        sig.level = 0.05, 
                                        d = 1)$power, 
                aes(color = "Effect size = 1")) +
  geom_function(fun = ~ pwr::pwr.t.test(n = .x, 
                                        sig.level = 0.05,
                                        d = 0.3)$power, 
                aes(color = "Effect size = 0.3")) +
  labs(x = "Sample number", y = "beta")

# 效果量和檢定力的關係(alpha 為0.05)
ggplot() + 
  xlim(c(0.1, 1.5)) + 
  geom_function(fun = ~ pwr::pwr.t.test(n = 10, 
                                        sig.level = 0.05, 
                                        d = .x)$power, 
                aes(color = "n = 10")) +
  geom_function(fun = ~ pwr::pwr.t.test(n = 50, 
                                        sig.level = 0.05, 
                                        d = .x)$power, 
                aes(color = "n = 50")) +
  geom_function(fun = ~ pwr::pwr.t.test(n = 100, 
                                        sig.level = 0.05,
                                        d = .x)$power, 
                aes(color = "n = 100")) +
  labs(x = "Effect size", y = "beta")

通常P會約定成俗的使用0.05。

做臨床試驗之前我們會先需要知道效果量(需要一個小的先行研究),知道效果量以後我們才知道要收集多少樣本(決定N是多少),才能達到檢定力80%的標準(常見的標準)

P 值是怎麼計算的?

參考:

Biostatistics for Dummies, Chapter 12, P161~
Wiki:https://en.wikipedia.org/wiki/Student%27s_t-test

以two sample T test為例, P值是怎麼計算的?

單用講的實在抽象,雖然現在有統計軟體的幫忙,我們通常不需要知道P值是怎麼被計算出來的,但實際自己算算看,不但能留下更深刻的印象,也能幫助我們理解剛剛所呈現統計結果的摘要。

計算兩組之間的差異

set.seed(123)
a <- rnorm(10,10)
b <- rnorm(10,10.5)

(d <- mean(a) - mean(b))

## [1] -0.6339963

計算我們得到的平均(或其他推定值)的(不)精確度(Standard error; SE):

比較兩個數值的平均我們會需要計算pooled standard deviation; S_p = \(\sqrt{(SD_a ^ 2 + SD_b^2)/2}\), 再除以\(\sqrt{(n/2)}\)

(SP <- sqrt((sd(a)^2 + sd(b)^2)/2))

## [1] 0.99682

(SE <- SP/sqrt(10/2))

## [1] 0.4457915

計算t-statistic:差異/SE

(t <- d/SE)

## [1] -1.422181

計算自由度

two sample test 的自由度是n_a + n_b - 2。這邊減去2是因為為了估算平均,兩組sample各用掉了一個自由度

df <- length(a) + length(b) -2

利用t-分布計算P值

(p <- 2 * pt(q = t, df = df, lower.tail = T))

## [1] 0.1720714

Note:這邊要乘以二是因為我們做的是雙尾檢定, 沒有預先設定要檢定左側或右側,所以計算出來的P值要乘以二。

來比較看看手動計算的結果跟t.test給我們的結果吧!

p_of_t.test <- t.test(a,b, var.equal = T)$p.value
t_of_t.test <- t.test(a,b, var.equal = T)$statistic

all.equal(t, t_of_t.test, check.names = F)

## [1] TRUE

all.equal(p, p_of_t.test, check.names = F)

## [1] TRUE

如果我們做的是單尾檢定, P 會變成原本的一半。

one_tail_p <- t.test(a,b, var.equal = T, 
                     alternative = "less")$p.value
all.equal(p/2, one_tail_p, check.names = F)

## [1] TRUE

用t分布來求p很抽象, 我們再用圖說明一次。

# 雙尾的P(左加右的面積)
ggplot() + 
  geom_function(fun = ~ dt(x = .x, df = 18)) + 
  xlim(c(-3,3)) +
  geom_vline(aes(xintercept = c(-t, t))) +
  labs(x = "t-statistic", y = "density")

# 單尾的P
ggplot() + 
  geom_function(fun = ~ dt(x = .x, df = 18)) + 
  xlim(c(-3,3)) +
  geom_vline(aes(xintercept = c(t))) +
  labs(x = "t-statistic", y = "density")

# 自由度對t分布的影響
ggplot() + 
  geom_function(fun = dt, args = list(df = 1), aes(color = "df = 1")) + 
  geom_function(fun = dt, args = list(df = 10), aes(color = "df = 10")) + 
  geom_function(fun = dt, args = list(df = 50), aes(color = "df = 50")) + 
  geom_function(fun = dnorm, aes(linetype = "normal dist.")) + 
  xlim(c(-3,3)) + labs(x = "t-statistic", y = "density")

map(c(3,5,10,100), ~ ggplot() + 
  geom_function(fun = dt, args = list(df = .x)) + 
  scale_x_continuous(limits = c(-4,4), breaks = -4:4) +
  geom_vline(aes(xintercept = c(qt(p = 0.025, df = .x), 
                                - qt(p = 0.025, df = .x)))) +
  geom_text(aes(0, 0.1, 
                label = str_c("p = 0.05,\n|t|= ", 
                              - round(qt(p = 0.025, df = .x), 2)))) +
  labs(x = "t-statistic", y = "density", title = str_c("df = ", .x))) %>% 
  ggpubr::ggarrange(plotlist = .)

Superiority, equivalence and non-inferiority test

參考:

Biostatistics for Dummies, Chapter 16
Turan, F. N., & Senocak, M. (2007). Evaluating “superiority”, “equivalence” and “non-inferiority” in clinical trials. Annals of Saudi medicine, 27(4), 284–288. https://doi.org/10.5144/0256-4947.2007.284

Superiority test 是為了證明一個藥物的臨床效果比另一個好(或是比安慰劑好), 通常會使用two sided test, alpha value 設在0.05。(也就是我們剛剛在做的事情。)

但有時候我們想要證明的不是A>B, 而是A=B或是A>=B。我們能夠用拒絕虛無假說A=B來證明A>B或A<B,但我們要怎麼證明A=B?

Q:你覺得當拒絕A=B的虛無假說時,P>0.05能夠證明A=B嗎?為什麼不能?

這個時候就需要進行Equivalence 或是Noninferiority Test.。

我們會需要先建立一個可以接受的範圍，如果這個藥物的95% CI能在這個範圍之內，我們說她的效能跟之前的藥物相同或不比他差。

Equivalence

這邊的CI通常是90%,會比95%CI短一點

Noninferiority Test

可以是90或是95%CI

結尾: 其他推薦書單

Wickham H, Grolemund G. R for data science: import, tidy, transform, visualize, and model data. ” O’Reilly Media, Inc.”; 2016 Dec 12.
Faraway, Julian J. Linear models with R. ed 2, Chapman and Hall/CRC, 2015.
Faraway JJ. Extending the linear model with R: generalized linear, mixed effects and nonparametric regression models. Chapman and Hall/CRC; 2016 Mar 23.
Gareth J, Daniela W, Trevor H, Robert T. An introduction to statistical learning: with applications in R. Spinger; 2013.
Wickham, H. (2016). ggplot: Elegant Graphics for Data Analysis.
Crawley, M. J. (2012). The R book. John Wiley & Sons.

P.value, power analysis and superior test