Werden wir Helden für einen Tag
Home | About | Archive現實的 R: descriptive statistics
根據非正式統計，本博講 R 或 Stat 的話 ((另外三個題目是小說、讀書及棋類遊戲)) ，讀者人數會銳減。就算本人講三腳貓的電腦編程，或者政治塗鴉，關注人數會較多。
但博客本來的精神，是忠於自己的偉大。對於一個下半生都應該會投入統計事業騙口苦飯吃的人來說，沒有比在博客分享統計電算之樂更正常的事。 ((其實是因為我這個人是沒有生活的，像個活死人，除了講這些東西之外，我的生命空虛得像個荒廢的防空洞。我的博客是沒有精彩如重遇當年暗戀的女孩或者飲飲食食之類的題才。))
統計的其中一個重要部份，是了解數據。 Descriptive statistics 是第一步。
通常論文第一個 table ，是一個像是這樣的 Table 。 ((取自 Gault et al. Effect of oxandrolone and timing of pubertal induction on final height in Turner’s syndrome: randomised, double blind, placebo controlled trial. BMJ 2011; 342:d1980))
上碩士課時，教授曾說過，在研究只在寫 Protocol 的階段，也已經在 Word 或 excel 設計好空的表格，格式就是像上面的 Table 。在研究完成之時，只消將數字填入。他說這方法可以令他清楚自己的研究在做甚麼，自己想要看到甚麼數據及統計結果，也可以防止自己被引誘去幹 data dredging 的勾當。
其實 SPSS 也有輸出表格的功能，只需將表格 Copy and paste 到 Word 就可以了。但這些功能有些人會濫用，我的 Thesis supervisor 就常常說有博士研究生將 SPSS 的 raw output 如那些 Correlation matrix 抄到 Word 交貨。但是看了半天，都不知道在幹甚麼。
我通常的方法，都是用 Word 打好 Table ，再從 R 產生結果，再將 Table 用人手填滿。這個方法並沒有甚麼問題的，只是有時覺得慢。
最近在找有沒有自動產生這些 Table 的方法，我並不介意 Table 產生後要我用人手改 Description ，但我的要求是 Table 上的統計結果要是自動產生。
其實一早知道有類似的東西。例如晨早知道內建有 by() 這個 function 。以 R 內建的費雪的 iris 數據作例子。
#< -
by(iris$Sepal.Length, iris$Species, FUN=mean)
但是回來的結果如此：
iris$Species: setosa
[1] 5.006
------------------------------------------------------------
iris$Species: versicolor
[1] 5.936
------------------------------------------------------------
iris$Species: virginica
[1] 6.588
回報的東西叫做 list ，並不像一個 table 。我有想過寫一個 function ，將這個 list 砌成 table 模樣，但通常想要做而又如此常用的東西，理應有人做了才是，只不過是我沒有去找而已。找了一回，找到了 doBy 這個套件。其中 summaryBy() 此一 function 就是我想要的東西。
讀過 help 後，可以測試。
#< -
summaryBy(Sepal.Length+Sepal.Width~Species, data=iris, FUN=mean, keep.names=TRUE)
回來的東西是：
Species Sepal.Length Sepal.Width
1 setosa 5.006 3.428
2 versicolor 5.936 2.770
3 virginica 6.588 2.974
Check 過是 dataframe 。但慣常這類 table 是以類別為直行，簡單而言，以上的這個 table 要 transpose 。
data.frame(t(summaryBy(Sepal.Length+Sepal.Width~Species, data=iris, FUN=mean, keep.names=TRUE)))
X1 X2 X3
Species setosa versicolor virginica
Sepal.Length 5.006 5.936 6.588
Sepal.Width 3.428 2.770 2.974
頗像樣。但通常這些 table 又會是以 mean (SD) 的方式表達數據。那麼只要改改 FUN 這個 parameter 便可以。
#< -
meansd <- function (x, sig.fig=3) {
return(paste(round(mean(x),sig.fig),"(",round(sd(x),sig.fig),")"))
}
data.frame(t(summaryBy(Sepal.Length+Sepal.Width~Species, data=iris, FUN=meansd, keep.names=TRUE)))
X1 X2 X3
Species setosa versicolor virginica
Sepal.Length 5.006 ( 0.352 ) 5.936 ( 0.516 ) 6.588 ( 0.636 )
Sepal.Width 3.428 ( 0.379 ) 2.77 ( 0.314 ) 2.974 ( 0.322 )
雖然不知為何有 X1, X2, X3 ，但那些都不太要緊，切去它亦不太難。根據同樣的方法，只要自行撰寫如 meansd() 之類的 function ，就可以產生如 median (IQR) 或 N (%) 之類的 summary statistics 。
既然可以生產這樣的 data frame ，代表可以用 write.csv() 之類的東西將它存為 csv ，在試算表開啟。
甚至可以用 odfWeave 來生產文書處理檔案，方便編輯。
例如在 OpenOffice.Org Writer 打入如此 codechunk
假設將此文件存為 table1.odt ，再用 R 的 odfWeave 處理。
library(odfWeave)
odfWeave("table1.odt","output_t1.odt")
出來的東西是這個樣。
