An R function: OLS/Robust scaled mass index

R code

R 繪圖：柱狀圖的低階繪圖範例

R 繪圖時的柱狀圖如何加上誤差線？如何使用線條填色？如何把圖示放在外部？本文提供一個例子。

以下是產生上圖的原始碼。

「在 95% 信心水準，抽樣誤差於正負 3.1% 以內」到底是什麼意思？

如果你對統計學有點興趣，可以往下看看這篇文。不知不覺寫得有點長。

故事：「某民調訪問 1000 名受訪者，有 40% 比例的人支持；在 95% 信心水準，抽樣誤差於正負 3.1% 以內。」現實中，民調常常出現相似的句子。我相信很多人看不懂這句話到底是什麼意思，或是誤解了這句話。所以現在來聊聊這到底是什麼鬼好了。

我們先把上面這個故事改成另一種情境好了。

例子：「你有一個大袋子。袋子裡有無數顆球。球只可能是黑色或白色。你希望知道一件事：到底黑球在大袋子中的個數比例是多少。現在，你抽出 1000 顆球，然後你發現有 400 顆是黑球（40%）。在 95% 信心水準，抽樣誤差於正負 3.1% 以內，也就是黑球比例的區間估計是 40% − 3.1% 到 40% + 3.1%。」

如果你可以看出這個例子其實和民調結果的故事是同一件事，那就可以往下看了。

在例子中，每一次調查就是抽出 n 顆球（例子中的 n = 1000）。袋子中黑球的比例，也就是我們希望知道的未知數，我們叫 p 好了。因為球實在太多了，我們不可能知道 p 的大小，只有抽出 n 顆球來估 p 的大小。

在統計學理，我們可以在還沒抽出球之前就預知「每次調查有多準」。什麼叫「有多準」呢？統計學的辦法是決定一個機率和設計一個區間，並估計 p 有多大的機率會座落在這個設計出來的區間之內。如果你每次調查抽出的球數 n 很大很大，那這個區間會很容易計算。

有點難懂？我試著換個方式說。例如，如果我可以做 100 次相同調查，我將會得到 100 個黑球比例和 100 個區間。我可以保證，這 100 次調查中，約有 95 個區間會包括了那個我們不可能知道的 p。這就是「在 95% 信心水準」的實際意義。

你也許會問，可不可能算 100 個區間有 99 個區間會包括 p 的調查？可以，但每個區間大小會變寬（由 ± 3.1% 變成 ± 4.1%），但結果也變得更值得信賴（因為區間將有 99% 的機率可以包括了真實的 p）。

如果你每次抽出更多顆球，那每個區間也都會變窄（因為抽出越多球就準嘛）。例如，每次抽出 5000 顆球，那這個區間大小會由 ± 3.1% 變成 ± 1.4%。

回到原本的例子。如果以相同的調查方法（一樣訪問 1000 名受訪者的情況下），並且重覆很多很多次（當然，實際上只進行了一次）。我們可以保證，那個我們想估計的支持度將座落在 95% 次所算出的區間之內。在故事中，這次調查結果的區間是 40% − 3.1% 到 40% + 3.1%。我們 *永遠不可能* 知道這次得到的區間會不會包括了真實的支持度 p，但它很可能會包括，因為這個調查方法就保證了進行 100 次約有 95 次是成立的。

有些人會把 40% − 3.1% 到 40% + 3.1% 理解成因為調查進行所造成的誤差，例如打錯資料、拒絕受訪、無效結果之類的。這都不是信賴區間的意義。

我並不想把這中間的數學式寫出來啦。只想說個概念。

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

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

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

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

最近回應了一篇在 PTT 統計板的問題：「[問題] 顯著水準的意思」。

※ 引述《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 

The results show that:

• $\sum_{i=1}^{11} (x_{1i}-x_{2i}) / 11 - 10 = 8.45455$;
• Totally $2048 = 2^{11}$ permutations are processed;
• 445 permutations do not support $H_0$;
• $P = 445/2^{11} = 0.2172852$.

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

> hist(test2)

As you can see, the vertical lines $x=8.45455$ shows the critical boundary of rejecting $H_0$. Note that this is a right-tail test.

Again, you may call the returned object:

> str(test2)
List of 14
 $ x1          : num [1:11] 92 0 72 80 57 76 81 67 50 77 ...
 $ x2          : num [1:11] 43 67 64 64 51 53 53 26 36 48 ...
 $ x1.name     : chr "x1"
 $ x2.name     : chr "x2"
 $ n           : int 11
 $ mu          : num 10
 $ test.0      : num 8.45
 $ is.onesample: logi FALSE
 $ alternative : chr "g"
 $ test.perm   : num [1:2048] 8.45 1.36 22.45 15.36 8.82 ...
 $ DF          :'data.frame': 11 obs. of  3 variables:
  ..$ x1  : num [1:11] 92 0 72 80 57 76 81 67 50 77 ...
  ..$ x2  : num [1:11] 43 67 64 64 51 53 53 26 36 48 ...
  ..$ diff: num [1:11] 39 -77 -2 6 -4 13 18 31 4 19 ...
 $ N           : num 2048
 $ p.value     : num 0.217
 $ rejected.N  : int 445
 - attr(*, "class")= chr "exactOneOrPairedSampleTest"

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

An R function for bootstrap confidence intervals

Here I introduce a R function made by myself to calculate several confidence intervals of mean/median including exact CI, basic bootstrap CI, percentile bootstrap CI & studentized bootstrap CI.

R code

Arguments

There are five arguments in this R function bootCI():

• x — a numeric vector specifying $x_i$
• alpha = 0.05 — a numeric specifying $\alpha$ value
• alternative = c("t", "l", "g") — a single character string specifying two-sided, less or greater tail
• B = 1999 — a integer specifying number of bootstrap samples
• quantileAlgorithm = 7 — a integer specifying the argument type passed to quantile()

Examples

We can define a numeric vector $x = [5\, 2\, 3\, 6\, 8\, 19\, 1.5]$ including $N=7$ items as num:

> num <- c(5, 2, 3, 6, 8, 19, 1.5)

After loading the above function, we can call the function to calculate CIs:

> bootCI(num)
Summary of x
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.500   2.500   5.000   6.357   7.000  19.000 
CIs of mu
                     2.5%     97.5%
$CI.exact       0.7778728 11.936413
$CI.basic       1.7125000  9.714286
$CI.percentile  3.0000000 11.001786
$CI.studentized 2.6785623 17.669412

The above results show that

\begin{aligned} \text{exact CI} & = ( \hat{\mu} + t_{\frac{\alpha}{2}} \cdot \hat{se}_{\hat{\mu}} \;,\; \hat{\mu} + t_{1-\frac{\alpha}{2}} \cdot \hat{se}_{\hat{\mu}} ) \\ \text{basic CI} & = ( 2\hat{\mu} - \hat{\mu}_{1-\frac{\alpha}{2}}^{*} \;,\; 2\hat{\mu} - \hat{\mu}_{\frac{\alpha}{2}}^{*} ) \\ \text{percentile CI} & = ( \hat{\mu}_{\frac{\alpha}{2}}^* \;,\; \hat{\mu}_{1-\frac{\alpha}{2}}^* ) \\ \text{studentized CI} & = ( \hat{\mu} + t^*_{\frac{\alpha}{2}} \cdot \hat{se}_{\hat{\mu}} \;,\; \hat{\mu} + t^*_{1-\frac{\alpha}{2}} \cdot \hat{se}_{\hat{\mu}} ) \\ \end{aligned}

where $\hat{\mu}$ is mean of $x$, $t_{\frac{\alpha}{2}}$ is lower $\alpha/2$ critical value for the $t$ distribution given $\text{df} = N-1$, $\hat{se}_{\hat{\mu}}$ is the standard error of the mean of $x$, $\hat{\mu}_{\frac{\alpha}{2}}^*$ is $\alpha/2$ percentile of the mean of bootstrapped $x$, and $t^*_{\frac{\alpha}{2}} = \frac{\hat{\mu^*}-\hat{\mu}}{\hat{se}^*_{\hat{\mu^*}}}$ is t-value based on bootstrapped $x$. See

• J. Carpenter and J. Bithell (2000). Bootstrap confidence intervals: when, which, what? A practical guide for medical statisticians. Statist. Med. 19:1141-1164
• Wikipedia contributors, 'Bootstrapping (statistics)', Wikipedia, The Free Encyclopedia, 17 December 2013, 21:28 UTC, <http://en.wikipedia.org/w/index.php?title=Bootstrapping_(statistics)&oldid=586549893> [accessed 31 December 2013]

for more details.

We can assign a different $\alpha$, direction or number of bootstrap samples:

> bootCI(num, alpha=0.01)
Summary of x
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.500   2.500   5.000   6.357   7.000  19.000 
CIs of mu
                      0.5%    99.5%
$CI.exact       -2.0962642 14.81055
$CI.basic       -0.1428571 10.28571
$CI.percentile   2.4285714 12.85714
$CI.studentized  1.1394952 25.22030

> bootCI(num, alternative="g")
Summary of x
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.500   2.500   5.000   6.357   7.000  19.000 
CIs of mu
                      5% 100%
$CI.exact       1.926445  Inf
$CI.basic       2.714286  Inf
$CI.percentile  3.285714  Inf
$CI.studentized 3.332740  Inf

> bootCI(num, B=99)
Summary of x
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.500   2.500   5.000   6.357   7.000  19.000 
CIs of mu
                     2.5%     97.5%
$CI.exact       0.7778728 11.936413
$CI.basic       1.6053571  9.330357
$CI.percentile  3.3839286 11.108929
$CI.studentized 2.5532298 15.532786

The returned list may be helpful for someone:

> myCI <- bootCI(num)
> myCI$CI.percentile
    2.5%    97.5% 
 3.00000 10.85714 
> str(myCI)
List of 8
 $ x             : num [1:7] 5 2 3 6 8 19 1.5
 $ alpha         : num 0.05
 $ alternative   : chr "t"
 $ B             : num 1999
 $ CI.exact      : Named num [1:2] 0.778 11.936
  ..- attr(*, "names")= chr [1:2] "2.5%" "97.5%"
 $ CI.basic      : Named num [1:2] 1.86 9.71
  ..- attr(*, "names")= chr [1:2] "2.5%" "97.5%"
 $ CI.percentile : Named num [1:2] 3 10.9
  ..- attr(*, "names")= chr [1:2] "2.5%" "97.5%"
 $ CI.studentized: Named num [1:2] 2.65 17.73
  ..- attr(*, "names")= chr [1:2] "2.5%" "97.5%"
 - attr(*, "class")= chr "bootCI"

Paired t-test vs general mixed model with R

Dataset

Consider the following dataset:

> x1 <- c(3,4,3,5,6,6,6,7,7,6,8,9,11)
> x2 <- c(3,4,5,5,8,6,8,7,9,8,8,9,14)
> cbind(x1,x2)
      x1 x2
 [1,]  3  3
 [2,]  4  4
 [3,]  3  5
 [4,]  5  5
 [5,]  6  8
 [6,]  6  6
 [7,]  6  8
 [8,]  7  7
 [9,]  7  9
[10,]  6  8
[11,]  8  8
[12,]  9  9
[13,] 11 14

where x1 and x2 are paired sample.

Paired t-test

The process of a typical paired t-test with R is following:

> t.test(x1, x2, paired=T)

 Paired t-test

data:  x1 and x2
t = -3.1225, df = 12, p-value = 0.008814
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.6977786 -0.3022214
sample estimates:
mean of the differences 
                     -1

The results show that 1) average of $x_1 - x_2$ equals to -1; 2) $t_{12} = -3.1225$, $p = 0.008814$.

Mixed model for paired sample

We can build a mixed model for this paired-sample example. Consider the following dataset:

> y <- c(x1,x2)
> trt <- gl(2,13,26)
> subject <- gl(13,1,26)
> cbind(y, trt, subject)
       y trt subject
 [1,]  3   1       1
 [2,]  4   1       2
 [3,]  3   1       3
 [4,]  5   1       4
 [5,]  6   1       5
 [6,]  6   1       6
 [7,]  6   1       7
 [8,]  7   1       8
 [9,]  7   1       9
[10,]  6   1      10
[11,]  8   1      11
[12,]  9   1      12
[13,] 11   1      13
[14,]  3   2       1
[15,]  4   2       2
[16,]  5   2       3
[17,]  5   2       4
[18,]  8   2       5
[19,]  6   2       6
[20,]  8   2       7
[21,]  7   2       8
[22,]  9   2       9
[23,]  8   2      10
[24,]  8   2      11
[25,]  9   2      12
[26,] 14   2      13

where y is dependent variable, trt is a two-level fixed factor, and subject is experimental units considered as blocks (random factor). A mixed model can fit this dataset:

\[ y_{ij} = \tau_i + b_{j} + \varepsilon_{ij} \]

where $i=\{1,2\}$ showing two-level in fixed factor $\tau$, $j=\{1,2,\ldots , 13\}$ showing 13 subjects as a random factor $b$.

We can find the solution of this mixed model by using lme4:lmer():

> require(lme4)
> mod <- lmer(y ~ trt + (1|subject) )
> summary(mod)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ trt + (1 | subject) 

REML criterion at convergence: 98.5708 

Random effects:
 Groups   Name        Variance Std.Dev.
 subject  (Intercept) 5.8590   2.4205  
 Residual             0.6667   0.8165  
Number of obs: 26, groups: subject, 13

Fixed effects:
            Estimate Std. Error t value
(Intercept)   6.2308     0.7085   8.794
trt2          1.0000     0.3203   3.123

Correlation of Fixed Effects:
     (Intr)
trt2 -0.226

As you can see, the estimate of trt2 equals 1 showing the average of $x_2 - x_1$; the statistic $t=3.123$ is equivalent to the t-value I get in paired t-test.

We can find the p-value by using car:Anova():

> require(car)
> Anova(mod, type=3, test.statistic="F")
Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)

Response: y
                F Df Df.res    Pr(>F)    
(Intercept) 77.34  1 13.288 6.636e-07 ***
trt          9.75  1 12.000  0.008814 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

and then we can find that $F_{1,12} = 9.75 = 3.123^2$ and $p=0.008814$. These outcomes are all equivalent to what I got in paired t-test.

Multiple comparisons for generalized linear model in R

An example

Defining y and x

> y <- c (rpois(10,1), rpois(10,2), rpois(10,3), rpois(10,5))
> y
 [1] 0 1 0 2 1 2 0 1 1 0 3 2 2 4 2 3 3 1 1 1 3 0 4 4 2 3 1 4 3 6 3 8 6 3 2 6 5 5 2 8
> x <- gl(4,10)
> x
 [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4
Levels: 1 2 3 4

Fitting a Poisson regression

> f <- glm(y~x, family=poisson)
> summary(f)

Call:
glm(formula = y ~ x, family = poisson)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.44949  -0.90724   0.04533   0.52699   1.52242  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.2231     0.3536  -0.631 0.527945    
x2            1.0116     0.4129   2.450 0.014277 *  
x3            1.3218     0.3979   3.322 0.000895 ***
x4            1.7918     0.3819   4.692 2.71e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 67.967  on 39  degrees of freedom
Residual deviance: 34.884  on 36  degrees of freedom
AIC: 142.6

Number of Fisher Scoring iterations: 5

Multiple comparisons

> # install.packages("multcomp")
> require(multcomp)
> f.mc <- glht(f, linfct = mcp(x = "Tukey"))
> summary(f.mc)

  Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = y ~ x, family = poisson)

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)    
2 - 1 == 0   1.0116     0.4129   2.450  0.06433 .  
3 - 1 == 0   1.3218     0.3979   3.322  0.00481 ** 
4 - 1 == 0   1.7918     0.3819   4.692  < 0.001 ***
3 - 2 == 0   0.3102     0.2807   1.105  0.67729    
4 - 2 == 0   0.7802     0.2575   3.030  0.01198 *  
4 - 3 == 0   0.4700     0.2327   2.019  0.17303    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

> confint(f.mc)

  Simultaneous Confidence Intervals

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = y ~ x, family = poisson)

Quantile = 2.5471
95% family-wise confidence level
 

Linear Hypotheses:
           Estimate lwr      upr     
2 - 1 == 0  1.01160 -0.04002  2.06322
3 - 1 == 0  1.32176  0.30822  2.33529
4 - 1 == 0  1.79176  0.81905  2.76447
3 - 2 == 0  0.31015 -0.40481  1.02512
4 - 2 == 0  0.78016  0.12436  1.43596
4 - 3 == 0  0.47000 -0.12281  1.06282

> coef(f.mc)
    2 - 1     3 - 1     4 - 1     3 - 2     4 - 2     4 - 3 
1.0116009 1.3217558 1.7917595 0.3101549 0.7801586 0.4700036 
> vcov(f.mc)
              2 - 1         3 - 1         4 - 1         3 - 2         4 - 2         4 - 3
2 - 1  1.704542e-01  1.249996e-01  1.249996e-01 -4.545455e-02 -4.545455e-02  9.714451e-17
3 - 1  1.249996e-01  1.583330e-01  1.249996e-01  3.333333e-02  6.938894e-17 -3.333333e-02
4 - 1  1.249996e-01  1.249996e-01  1.458330e-01 -2.775558e-17  2.083333e-02  2.083333e-02
3 - 2 -4.545455e-02  3.333333e-02 -2.775558e-17  7.878788e-02  4.545455e-02 -3.333333e-02
4 - 2 -4.545455e-02  6.938894e-17  2.083333e-02  4.545455e-02  6.628788e-02  2.083333e-02
4 - 3  9.714451e-17 -3.333333e-02  2.083333e-02 -3.333333e-02  2.083333e-02  5.416667e-02

Read more

• Simultaneous Inferences in R — Using package:multcomp by J. C. Wang
• Simultaneous Inference in General Parametric Models by T. Hothorn et al.

在 R 中配置 Generalized Logit Model

已知資料 data.csv 如下:

treat, r.a, r.b, r.c
   t1,  10,  20,   5
   t2,  11,  22,   5
   t3,  21,  40,  11
   t4,  21,   5,  10

其中 treat 為具有 4 個水準的固定因子; r.a, r.b 及 r.c 分別為 3 個類型的反應類別, 其下為頻率.

以下為 R code:

# 系統第一次運作請先安裝 package nnet
install.packages("nnet")

# 引入 package nnet
require(nnet)

# 讀取資料
dat <- read.csv("data.csv")

# 建立 generalized logit model
m <- multinom(
  cbind(r.a, r.b, r.c) ~ treat,          # 其中 r.a 寫在最前面即可作為比較基準
  data=dat,
  contrasts=list(treat=contr.treatment), # 其中採用 treatment contrast 方便解釋
  Hess=T                                 # 其中採用 Hessian Matrix 求標準誤
)

# 係數與標準誤
(s.m <- summary(m))

# odds ratio
exp(s.m$coefficients)

# chi-square
(chisq <- (s.m$coefficients / s.m$standard.error)^2)

# p-value
pchisq(chisq, df=1, lower.tail=F)

可得到以下結果:

> (s.m <- summary(m))
Call:
multinom(formula = cbind(r.a, r.b, r.c) ~ treat, data = dat, 
    contrasts = list(treat = contr.treatment), Hess = T)

Coefficients:
    (Intercept)        treat2      treat3      treat4
r.b   0.6931308  2.158365e-05 -0.04877868 <strong>-2.12818025</strong>
r.c  -0.6931750 -9.525513e-02  0.04656602 -0.04876199

Std. Errors:
    (Intercept)    treat2    treat3    treat4
r.b   0.3872968 0.5351291 0.4718241 0.6305649
r.c   0.5477240 0.7687048 0.6622145 0.6690451

Residual Deviance: 346.6049 
AIC: 362.6049 


> # odds ratio
> exp(s.m$coefficients)
    (Intercept)   treat2    treat3    treat4
r.b   1.9999671 1.000022 0.9523919 <strong>0.1190537</strong>
r.c   0.4999861 0.909141 1.0476672 0.9524078


> # chi-square
> (chisq <- (s.m$coefficients / s.m$standard.error)^2)
    (Intercept)       treat2      treat3       treat4
r.b    3.202894 1.626795e-09 0.010688081 <strong>11.390878935</strong>
r.c    1.601630 1.535527e-02 0.004944712  0.005311928


> # p-value
> pchisq(chisq, df=1, lower.tail=F)
    (Intercept)    treat2    treat3      treat4
r.b  0.07350811 0.9999678 0.9176589 <strong>0.000738056</strong>
r.c  0.20567234 0.9013815 0.9439400 0.941899243

可以發現, 唯一達顯著之係數為 r.b 對比 r.a 時在 treat4 對比 treat1 的情況, 其係數 $\beta = -2.12818025$, 勝算比 $= \exp(-2.12818025) = 0.1190537$, 檢定量 $\chi^2_1 = 11.390878935$, $p = 0.000738056$, 顯示 treat4 比 treat1 的 r.b / r.a 達顯著差異. 由實際資料 $\frac{5/21}{20/10} = 0.1190476$ 也和模型所預測的勝算比 0.1190537 相當接近.

G-test for 2-way contingency data in R

A R function applying G-test for 2-way contingency data in R is demonstrated. Yates’ correction and Williams’ correction are also included in the output.

R Source

Example

R code

my.2way.contingency.table <- matrix(
  c(
    5,0,3,12,6,
    4,2,7,23,11
  ),
  nrow = 2, ncol=5, byrow=T
)

g.test.result <- g.test.2way(my.2way.contingency.table)
g.test.result
str(g.test.result)

R output

> my.2way.contingency.table <- matrix(
+   c(
+     5,0,3,12,6,
+     4,2,7,23,11
+   ),
+   nrow = 2, ncol=5, byrow=T
+ )
> 
> g.test.result <- g.test.2way(my.2way.contingency.table)
> g.test.result
Observed value:
     [,1] [,2] [,3] [,4] [,5]
[1,]    5    0    3   12    6
[2,]    4    2    7   23   11
Expected value:
         [,1]      [,2]     [,3]     [,4]      [,5]
[1,] 3.205479 0.7123288 3.561644 12.46575  6.054795
[2,] 5.794521 1.2876712 6.438356 22.53425 10.945205
Observed value with Yates' correction:
     [,1] [,2] [,3] [,4] [,5]
[1,]  4.5  0.5  3.5 12.5  6.5
[2,]  4.5  1.5  6.5 22.5 10.5
Degree of freedom = 4 
G-test:
 G = 3.41127, p = 0.491497
G-test with Yates' correction:
 G = 1.28744, p = 0.863503
G-test with Williams' correction:
 G = 3.07349, q_min = 1.1099, p = 0.545603
> str(g.test.result)
List of 11
 $ data.observed      : num [1:2, 1:5] 5 4 0 2 3 7 12 23 6 11
 $ data.expected      : num [1:2, 1:5] 3.205 5.795 0.712 1.288 3.562 ...
 $ data.observed.Yates: num [1:2, 1:5] 4.5 4.5 0.5 1.5 3.5 6.5 12.5 22.5 6.5 10.5
 $ df                 : num 4
 $ q.min              : num 1.11
 $ g                  : num 3.41
 $ g.Yates            : num 1.29
 $ g.Williams         : num 3.07
 $ p                  : num 0.491
 $ p.Yates            : num 0.864
 $ p.Williams         : num 0.546
 - attr(*, "class")= chr "g.test.2way"

Convert taxa matrix to guild matrix with R

在 R 中如何把分類群資料小計成為功能群資料? 可利用以下介紹的 R function taxaToGuild() 將一分類群資料檔及一分類群對應功能群資料檔轉化出一個功能群資料框. 一般來說, 該 function 也提供一個資料框資料以系統化的方式進行欄小計.

taxaToGuild.R 原始檔

例子

假設有一蜘蛛以科為分類群的物種數資料，以 CSV 檔儲存後如下 (taxa.csv):

"Plot","Leptonetidae","Clubionidae","Others","Araneidae","Oonopidae","Ctenizidae"
1,2,0,0,1,1,2
2,1,2,3,1,0,0
3,1,0,1,1,1,1
4,2,3,2,1,0,1
5,0,1,1,2,1,4
6,2,0,0,0,1,3
7,1,0,1,2,2,1
8,0,0,1,2,3,1

其中第一列皆為變數名且第一欄為可辦視的樣點名。另有一科名對照功能群的資料檔以 CSV 檔儲存後如下 (guild.csv):

"taxa","Guild1","Guild2","Guild3"
"Oonopidae","Ground runner","runner","Ground"
"Araneidae","Orb weaver","weaver","Space"
"Clubionidae","Foliage runner","runner","Foliage"
"Ctenidae","Ground runner","runner","Ground"
"Ctenizidae","Burrow dweller","dweller","Burrow"
"Gnaphosidae","Ground runner","runner","Ground"
"Leptonetidae","Ground weaver","weaver","Ground"
"Linyphiidae","Space weaver","weaver","Space"

其中第一列皆為變數名且第一欄為可辦視的科名, 第 2 至 4 欄為三套不同的功能群系統. 在將 taxa.csv, guild.csv, taxaToGuild.r 三檔案置於相同某路徑後，於 R 環境中以 setwd("某路徑") 後進行以下動作.

source("taxaToGuild.r")
new.dataframe <- taxaToGuild(
  taxa.file.csv = "taxa.csv",
  taxa.file.csv.row.num = 1,
  guild.file.csv = "guild.csv",
  guild.file.csv.row.num = 1,
  guild.system.row.num = 1
)
new.dataframe

其中 taxa.file.csv = "taxa.csv" 表示引入 taxa.csv 為分類群資料, taxa.file.csv.row.num = 1 表示 taxa.csv 的第一欄為列名, guild.file.csv = "guild.csv" 表示引入 guild.csv 為功能群對照檔, guild.file.csv.row.num = 1 表示 guild.csv 的第一欄為列名, guild.system.row.num = 1 表示以 guild.csv 的第 1 組對照表 (實際上是此例的第 2 欄). 新產生的功能群資料框賦予變數 new.dataframe. 操作可得以下結果.

> source("taxaToGuild.r")
> new.dataframe <- taxaToGuild(
+   taxa.file.csv = "taxa.csv",
+   taxa.file.csv.row.num = 1,
+   guild.file.csv = "guild.csv",
+   guild.file.csv.row.num = 1,
+   guild.system.row.num = 1
+ )
The guild system is Guild1.
Warning: the following taxa is merged into guild NA due to no reference: Others.
> new.dataframe
  Ground weaver Foliage runner NA Orb weaver Ground runner Burrow dweller
1             2              0  0          1             1              2
2             1              2  3          1             0              0
3             1              0  1          1             1              1
4             2              3  2          1             0              1
5             0              1  1          2             1              4
6             2              0  0          0             1              3
7             1              0  1          2             2              1
8             0              0  1          2             3              1

值得注意的是, 因為 "Others" 無法藉由 "Guild1" 這項對應表成功對應至某個功能群, 所以被丟棄至 new.dataframe 中的 "NA" 欄. 在本例中若欲刪去 new.dataframe 中的 "NA" 欄, 可利用 new.dataframe[,-3].

One-dimensional interpolation with R

The function approx() and spline() provide one-dimensional linear and non-linear interpolation with R, respectively. Syntax of them are similar, for example,

x.old <- 1:10
y.old <- rnorm(10)
x.new <- seq(1, 10, 0.001)

# Linear interpolation
y.new.linear <- approx(x.old, y.old, xout = x.new)$y

# Forsythe, Malcolm and Moler's spline (default of spline)
y.new.fmm <- spline(x.old, y.old, xout = x.new, method = "fmm")$y

# Natural splines
y.new.natural <- 
  spline(x.old, y.old, xout = x.new, method = "natural")$y

# Periodic splines 
y.new.periodic <-
  spline(x.old, y.old, xout = x.new, method = "periodic")$y

# Plot
plot(
  y.new.fmm ~ x.new,
  type="l", col=2, xlab="x", ylab="y",
  main="One-dimensional Interpolation with R"
)
lines(y.new.periodic ~ x.new, col=3)
lines(y.new.natural ~ x.new, col=4)
lines(y.new.linear ~ x.new, col=5)
points(y.old ~ x.old)
legend(8, 2, c("fmm", "periodic", "natural", "linear", "old"),
  col = c(2:5,1), lty = c(1,1,1,1,0), pch = c(NA,NA,NA,NA,1), merge = T
)

Nested ANOVA permutation test in R

A R function for the permutation test for onw-way nested design ANOVA.

This work is inspired by Dr. David C. Howell, University of Vermont. In Dr. Howell’s webpage “Permutation Tests for Factorial ANOVA Designs” (http://www.uvm.edu/~dhowell/StatPages/More_Stuff/Permutation%20Anova/PermTestsAnova.html [access on May 28, 2012]), he described the algorithm of permutation test for a one-way nested design ANOVA by using R language. I followed this algorithm and made a R function to do this test.

Example

dat <- read.csv(textConnection("
 trt , unit , obs
 1 , 1 , 11
 1 , 1 ,  9
 1 , 1 ,  9
 1 , 2 ,  8
 1 , 2 ,  7
 1 , 2 ,  6
 1 , 3 ,  8
 1 , 3 , 10
 1 , 3 , 11
 2 , 4 , 11
 2 , 4 ,  8
 2 , 4 ,  7
 2 , 5 , 10
 2 , 5 , 14
 2 , 5 , 12
 2 , 6 ,  9
 2 , 6 , 10
 2 , 6 ,  8
"))

## traditional nested ANOVA
mod.2 <- aov(
 obs ~ factor(trt) + Error(factor(unit)),
 data = dat
)
summary(mod.2) # alternative

## permutation nested ANOVA
nestedPermutationAnova(dat$obs, dat$trt, dat$unit, 499)

Logistic regression with package rms in R

The rms package by Frank E Harrell Jr provides a useful function lrm to apply logistic regression. Moreover, this function also provides many useful statistics, and a gooeness-of-fit test (le Cessie and Houwelingen test) could be applied by using the function residuals(a lrm model, "gof"). Read the official document of package rms for more details.

淺談貝氏定理

本文將介紹貝氏定理 (Bayes’ theorem) 的推導與實例。

推導

條件機率的定義是, \[ \Pr(A|B)= \frac{\Pr(A \cap B)}{\Pr(B)}, \] 其中 $\Pr(A|B)$ 是指在事件 $B$ 發生的前題下, 事件 $A$ 發生的機率; $\Pr(A \cap B)$ 是指事件 $A$ 發生且事件 $B$ 發生的機率. 當然, $A \cap B$ 與 $B \cap A$ 是完全同義的, 因此, 我們也可以得到 \[ \Pr(B|A)=\frac{\Pr(B \cap A)}{\Pr(A)}=\frac{\Pr(A \cap B)}{\Pr(A)}. \]

有趣的事發生了. 若將上述二式中的 $\Pr(A \cap B)$ 單離出來, 可以發現 \[ \Pr(A \cap B) = \Pr(A|B) \times \Pr(B) = \Pr(B|A) \times \Pr(A), \] 可得 \[ \Pr(A|B) = \frac{\Pr(B|A) \times \Pr(A)}{\Pr(B)}, \] 被稱為貝氏定理. 不過在實用上來說, 常利用全機率定理改寫分母, 成為定理的另一個形式 \[ \Pr(A|B) = \frac{\Pr(B|A) \times \Pr(A)}{\Pr(B|A)\times \Pr(A) + \Pr(B|\lnot A)\times \Pr(\lnot A) }, \] 其中 $\lnot A$ 表示 “非 $A$”.

例題

Wikipedia 中描述了一個貝氏定理中著名的 “吸毒測試” 例題, 簡單翻譯與修改後如下. 某種檢測吸毒的方法具有 99% 的敏感度和 97% 的可靠度, 也就是有吸毒者受測呈陽性的機率為 99% 而非吸毒者受測呈陰性的機率是 97%. 調查人員已經知道某群體的總體吸毒率為 0.1%. 已知, 在該群體中某一位人員受檢後呈陽性, 則該人員吸毒的機率是多少?

令吸毒事件為 $D$, 測試結果呈陽性為 $P$. 由命題可知, $\Pr(P|D)=0.99$, $\Pr(\lnot P|\lnot D)=0.97$, $\Pr(D)=0.001$. 由貝氏定義可知, \[ \begin{aligned} \Pr(D|P) =& \frac{\Pr(P|D) \times \Pr(D)}{\Pr(P)} \\ =& \frac{\Pr(P|D) \times \Pr(D)}{\Pr(P|D) \times \Pr(D) + \Pr(P|\lnot D) \times \Pr(\lnot D)} \\ =& \frac{\Pr(P|D) \times \Pr(D)}{\Pr(P|D) \times \Pr(D) + [1-\Pr(\lnot P|\lnot D)] \times [1-\Pr(\lnot D)]} \\ =& \frac{0.99 \times 0.001}{0.99 \times 0.001 + (1-0.97) \times (1-0.001)} \\ \simeq& 0.032 \end{aligned} \] 為所求, 也就是該人員受測呈陽性但有吸毒的機率只有 3.2%, 即使該檢測有 99% 的敏感度. 這種結果是不是與你預期的差距很大?

利用貝氏定理, 統計學家已經將之推廣出一派新的統計學門, 稱為貝氏推論. 貝氏推論的核心概念是, 在已知的資料中, 如何從過去的知識推論母體. 在例題中, 過去的知識 $\Pr(D)$ 被稱為事前機率 (prior), $\Pr(D|P)$ 被稱為事後機率 (posterior), 而 $\Pr(P|D) / \Pr(P)$ 被稱為標準概似度 (standardised likelihood). 因此, 貝氏定理亦常被表示成 \[ \text{posterior} = \text{standardised likelihood} \times \text{prior}. \] 假如 posterior 不易取得, 可藉由對 prior 更多的了解, 使我們可以更有信心地推論 posterior, 是貝氏推論的重要武器.

Create a NA matrix in R

matrix(data=NA, nrow=3, ncol=5) can create/initiate a 3 × 5 NA matrix on the fly. Of course, matrix(data=0, nrow=3, ncol=5) may be useful in case.

Statistical computing by using R

最近把東海大學生命科學系大學部的生物統計學實驗講義中的 SAS code 全改寫成 R code。內容是如何使用 R 計算常見的統計檢定，包括例題及 R code。該文件以 Creative Commons Attribution-ShareAlike 3.0 Unported License 授權。

View Statistical computing by using R on Scribd

A function of Pearson/Deviance goodness-of-fit for (generalized) linear regression

After an analysis of (generalized) linear regression, a Pearson/Deviance goodness-of-fit test is useful to test if the model fitted reasonably. I wrote a simple function in R to do it.

Source Code

#    Copyright 2011 Chen-Pan Liao 
#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    any later version.
#    
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#    
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.

# Usage example:
#   model.name <- glm (...)
#   lm.fit.test(model.name)

lm.fit.test <- function (my.model) {
  df <- my.model$df.residual;
  chisq.pearson <- sum(resid(my.model , type="pearson")^2);
  chisq.deviance <- sum(resid(my.model , type="deviance")^2);
  p.pearson <- pchisq( chisq.pearson , df , lower.tail=F );
  p.deviance <- pchisq( chisq.deviance , df , lower.tail=F );
  ratio.pearson <- chisq.pearson / df;
  ratio.deviance <- chisq.deviance / df;
  cat(
    "Pearson chisq = " , chisq.pearson ,
    ", df = " , df ,
    ", p = " , p.pearson ,
    ", chisq / df = " , chisq.pearson/df ,
    ".\n" ,
    "Deviance chisq = " , chisq.deviance ,
    ", df = " , df ,
    ", p = " , p.deviance ,
    ", chisq / df = " , chisq.deviance/df ,
    ".\n"
  );
}

Usage and example

After defining a glm() model as a custom name, copy the source code into R interpreter and then input lm.fit.test(your.model.name) where your.model.name is the custom glm model name. See the fallowing example

y <- c(rpois(20,1) , rpois(20,2));
x <- gl(2 , 20);
my.model <- glm(y~x , family=poisson);
summary(my.model);
lm.fit.test(my.model);

Games-Howell post-hoc in R

Mr. Shigenobu AOKI displayed a quick R function for both of Tukey and Games-Howell post-hoc comparisons after an ANOVA. The Games-Howell post-hoc is famous for non-equal variances between treatments. Read more on the webpage about Games-Howell post-hoc in R.