Apan’s Notes

《初學 R 語言的 60 分鐘》投影片

2015 年 3 月 14 日於蓮華池研究中心所舉行的台灣生態研究網年會中，我針對 R 語言的初學者所設計的課程投影片。

初學R語言的60分鐘 from Chen-Pan Liao

統計檢驗簡介：一個二項式分配的簡單例子

※ 引述《mrlee112233 (小史)》之銘言：
: 不好意思
: 我想問一下顯著水準
: 有沒有簡單一點的解釋
: 我上網查過資料
: 但因爲我不是相關科系
: 所以我跟本看不懂在寫什麼
: 虛無假設.type l error .type ll error..等
: 跟本不懂=   =
: 我文獻報告裏
: 他假設alpha=0.05 然後用二相分佈去做計算
: 所以想問有沒有比較淺顯易懂的解釋
: 謝謝

以下則是我的回應，也歡迎從 PTT 上我原本的回應觀看相同的內容。

我會建議你從教科書或聽課來學習這一連串的概念，但我也可以體會非相關科系的朋友不容易了解這些概念，所以寫一個例子給你參考，並省略一些太艱澀的話語。

假如有一個袋子，裡有無數顆球，其中球不是黑的就是白的。你可以獨立抽出 10 顆球，並記錄每顆球是什麼顏色。在抽出球之前，你可以設立一個虛無假說：「袋中的黑球和白球比例是 1:1」。當然，你並不知道這是不是對的，但就先這麼假設吧。

假如這個虛無假說是真實的，又因為你會抽出 10 顆球，所以你可以預先算出你抽出 0 顆黑球到 10 個顆黑球的機率。例如，抽出 0 顆黑球的機率是 0.0009766，抽出 1 顆黑球的機率是 0.009766，抽出 5 顆黑球的機率是 0.2461，…… 記得，這在還沒有抽球之前，就可以算得的。

接下來，你可以設立一個和虛無假說相對的假說，叫對立假說。我們就說這個對立假說是「袋中的黑球和白球比例不是 1:1」。這裡會跳出單尾和雙尾檢驗的概念，但我就不多說了，反正我們就做雙尾檢驗吧。另外，我們也要設定顯著水準，常見為 0.05。

這時候，你可以抽出 10 顆球了。如果你抽出 5 顆黑球和 5 顆白球，那你可能會相信虛無假說是對的。但如果你抽出 0 顆黑球和 10 顆白球，或是 10 顆黑球和 0 顆白球，那你可能會大大地懷疑虛無假說，而相信對立假說才是對的。問題來了：到底要多麼地違背虛無假說，你才相信對立假說？這就是顯著水準 = 0.05 的作用。

例如，如果你真的抽出 1 顆黑球和 9 顆白球好了。發生這種情況的機率是 0.009766（假如虛無假說是對的），另外，還有三種情況和 1 黑 9 白一樣程度或更違反虛無假說，分別是 0 黑 10 白（機率是 0.0009766）、9 黑 1 白（機率是 0.009766）、10 黑 0 白（機率是 0.0009766）。這四種情況的機率總共是 0.02148，稱為 p-value。

這時候，因為這個 p-value = 0.02148 比你設定的顯著水準小，所以你可以下一個結論：在 0.05 的顯著水準下，虛無假說不被接受。

當然，你也可能猜錯了，因為即使虛無假說是真的，你還是有個很小的機率抽到 1 黑 9 白或更偏激的結果。這種猜錯的情況就稱為型一錯誤。不過因為通常我們會設定一個滿小的顯著水準，所以型一錯誤不甚容易發生。

換一個情況，假如你抽出的是 2 黑 8 白呢？這時候，下列的所有情況的機率和就是 p-value：0 黑 10 白（機率是 0.0009766）、1 黑 9 白（機率是 0.009766）、2 黑 8 白（機率是 0.0439）、8 黑 2 白（機率是 0.0439）、9 黑 1 白（機率是 0.009766）、10 黑 0 白（機率是 0.0009766）。加起來的機率是 p-value = 0.1094，比顯著水準大了。這時候，你可以下一個結論：沒有證據指出虛無假說錯了，在 0.05 的顯著水準下。

再一次地，你也可能猜錯了，因為，說不定虛無假說並不正確（例如真實情況是 30% 黑球 70% 白球）。這種猜錯的情形就叫型二錯誤。要儘量避免型二錯誤是可能的（關鍵字：檢定力），但就不多說了。

整個故事其實不複雜，好啦，我寫得讓它變得有點複雜了。在虛無假說為真的條件下，取得目前及更偏離虛無假說的結果之機率叫 p-value，並拿它與顯著水準相比較。

因此，顯著水準被當做一個門檻、一個標準，來決定我們拒絕或不拒絕虛無假說。

我原本以為可以用最簡單的例子說明這一連串的概念，結果還是用了這麼多字。算是騙 p 幣好了。

An R function: Exact one/paired sample test for mean(s)

Here I provide a exact one sample or paired sample test for mean(s), which is the exact version of permutation test to compare a single sample against a mean or to compare paired samples’ differences against a mean.

The algorithm in this function is based on Bryan F. J. Manly. 1997. Randomization, bootstrap and Monte Carlo methods in biology. 2nd edition. pp. 91-97. That is, all possible permutations (i.e. all possible cases of exchange between $x_i$ and $\mu$ in one-sample test or all possible cases of exchange between $x_{1i}$ and $x_{2i}$) in paired-sample test are processed to calculate a exact p-value.

R code

Arguments

The usage of this function exactOneOrPairedSampleTest() is very similar to the R built-in function t.test(). There are four arguments:

x1: a numeric vector specifying $x_{1i}$
x2 = NULL: a numeric vector specifying $x_{2i}$ (in paired-sample case)
mu = 0: a number specifying true mean or true difference $\mu$
alternative = c("t","g","l")a single character specifying $H_0$: true mean of $x_1$ or mean of $x_{1i} - x_{2i}$ is equal, greater or less than $\mu$, respectively.

Example 1: one-sample test

Let $x_i = \{43,67,64,64,51,53,53,26,36,48,34,48,6 \}$, $i=1 \ldots 13$:

> x <- c(43,67,64,64,51,53,53,26,36,48,34,48,6)

Now we can call the above function exactOneOrPairedSampleTest() to test $H_0$: true mean of $x_i = 56$ against $H_A$: true mean of $x_i \neq 56$:

> test1 <- exactOneOrPairedSampleTest(x, alternative="t", mu=56)
> test1

 Exact one sample test

Alternative hypothesis: true mean is not equal to 56 

mean(x) - mu = -10.3846153846154
Number of total permutation = 8192 
Number of rejected permutation = 364 
P-value = 0.04443359

The results show that:

$\sum_{i=1}^{13} x_i / 13 - \mu = 45.61538 - 56 = -10.38462$;
Totally $8192 = 2^{13}$ permutations are processed;
364 permutations do not support $H_0$;
$P = 364/2^{13} = 0.0444$.

We can also call the function hist.exactOneOrPairedSampleTest():

> hist(test1)

As you can see, the vertical lines $x=\pm 10.38462$ shows the critical boundary of rejecting $H_0$. Note that it is a two-tail test.

Finally, you may call the returned object:

> str(test1)
List of 14
 $ x1          : num [1:13] 43 67 64 64 51 53 53 26 36 48 ...
 $ x2          : NULL
 $ x1.name     : chr "x"
 $ x2.name     : chr "NULL"
 $ n           : int 13
 $ mu          : num 56
 $ test.0      : num -10.4
 $ is.onesample: logi TRUE
 $ alternative : chr "t"
 $ test.perm   : num [1:8192] -10.38 -8.38 -12.08 -10.08 -11.62 ...
 $ DF          :'data.frame': 13 obs. of  3 variables:
  ..$ x1  : num [1:13] 43 67 64 64 51 53 53 26 36 48 ...
  ..$ mu  : num [1:13] 56 56 56 56 56 56 56 56 56 56 ...
  ..$ diff: num [1:13] -13 11 8 8 -5 -3 -3 -30 -20 -8 ...
 $ N           : num 8192
 $ p.value     : num 0.0444
 $ rejected.N  : int 364
 - attr(*, "class")= chr "exactOneOrPairedSampleTest"

to get more details of the results if you need them.

Example 2: paired-sample test

Let $x_{1i} = \{92, 0,72,80,57,76,81,67,50,77,90\}$ and $x_{2i} = \{43,67,64,64,51,53,53,26,36,48,34\}$, where each $i$ is paired and $i=1 \ldots 11$:

> x1 <- c(92, 0,72,80,57,76,81,67,50,77,90)
> x2 <- c(43,67,64,64,51,53,53,26,36,48,34)

Now we can call the function exactOneOrPairedSampleTest() to test $H_0$: true mean of $x_{1i} - x_{2i} \leq 10$ against $H_A$: true mean of $x_{1i} - x_{2i} > 10$:

> test2 <- exactOneOrPairedSampleTest(x1, x2, alternative="g", mu=10)
> # equivalent to test2 <- exactOneOrPairedSampleTest(x1 - x2, alternative="g", mu=10)
> test2
 Exact paired sample test

Alternative hypothesis: means of x1 - x2 is greater than 10
mean((x1)-(x2)) - mu = 8.45454545454546
Number of total permutation = 2048 
Number of rejected permutation = 445 
P-value = 0.2172852