R語言從入門到精通：Day10
- 2019 年 10 月 7 日
- 筆記
是時候
關注
我們一波了
到目前為止，R語言的數據操作和基礎繪圖部分已經講解完畢，換句話說，大家應該已經能將數據導入R中，並運用各種函數處理數據使其成為可用的格式，然後將數據用各種基礎圖形展示。完成前面這些步驟之後，我們接下來要探索數據中變數的分布以及各組變數之間的關係。
在課題或者項目中，你往往會遇到這樣的問題：參與本次實驗的病人的年齡的分布如何（均值、、標準差、中位數等）？實驗中不同組病人的生存時間有沒有差異？病人性別對實驗結果有無影響？接下來的幾次教程內容就是為了解決這些問題，我們會逐步學習R語言中的一些統計方法，希望大家在學習新內容的同事，也可以回顧一下自己之前學過的統計學課程。本次教程將主要關注R語言中生成基本的描述性統計量和推斷統計量的R函數。
寫在開篇的話，本篇教程內容較多，請務必靜下心來學習。
溫馨提示
1、本節內容重點內容較多，
務必緊跟紅色標記。
2、測試數據及程式碼
見文末客服小姐姐二維碼。
1、連續型變數的統計描述
生成描述性統計量的R函數中，連續型變數和類別型變數的統計方法有所不同，首先介紹連續型變數的統計函數（以R中自帶的mtcars數據集為例），summary()是R中基礎安裝的獲取描述性統計量的函數。函數summary()提供了最小值、最大值、四分位數和數值型變數的均值，以及因子向量和邏輯型向量的頻數統計。其中典型函數還有mean()、sd()、var()、min()、max()、median()、length()、range()和quantile()。同時，函數fivenum()可返回圖基五數總括(Tukey`s five-number summary，即最小值、下四分位數、中位數、上四分位數和最大值)。（上述函數的使用比較基礎，就不一一舉例了。）不過，R基礎安裝中沒有提供偏度和峰度的計算函數，下面是一個自定義計算偏度和峰度的函數實例。
圖1，偏度和峰度的示例。
圖1中，函數mystats()是自定義的函數（用於計算圖中所示的五個描述性統計量），函數sapply()和函數apply()使用類似，在之前的教程中介紹過。（具體程式碼見後台。）
當然，除了基礎安裝，還有很多獲取描述性統計量的方法。比如：Hmisc包中的describe()函數、pastecs包中的stat.desc()函數和psych包也擁有的describe()的函數。這裡我們給出Hmisc包的例子：（兩個包中的函數名稱重複時，在函數名前面加上包名稱即可，如Hmisc::describe()）。
圖2:Hmisc包中describe()示例
2、"分組"連續型變數的統計描述
上面介紹了獲取整體數據的描述性統計量的方法，更多時候我們需要將數據分組後分別計算各組的描述性統計量，函數by()或者aggregate()可以解決這個難題。下面是函數by()的一個例子，以變數am為分類標準，分別計算兩組的描述性統計量。其中函數dstats()是在函數mystats()基礎上定義的。
圖3:by()示例
當然也有很多函數包提供了分組計算描述性統計量的方法。比如：doBy包中summaryBy()函數、psych包中的describeBy()函數。下面是psych包中的describeBy()的示例：
圖4:describeBy()示例
描述性統計量的計算是很基礎的分析步驟，R中用於獲取描述性統計量的方法很多，大家可以根據自己的需要或者喜好選擇，或者你還可以自己寫一個函數出來！
3、分類變數的統計描述
對於連續型變數，我們可以計算均值、標準差等，那麼對於類別型變數該怎麼辦呢？頻數表和列聯表可以解決這個問題。（示例數據來自vcd包中的Arthritis數據集。）創建頻數表和列聯表的幾種重要方法如下表：
表1: 用於創建和處理列聯表的函數
具體的示例程式碼可以直接找客服胖雨小姐姐要（文末二維碼），就不在這裡一一展示了。使用gmodels包中的CrossTable()函數也是創建二維列聯表的一種方法，示例如下圖5.
圖5:函數CrossTable示例
函數CrossTable()有很多選項，可以做許多事情:計算(行、列、單元格)的百分比;指定小數位數;進行卡方、Fisher和McNemar獨立性檢驗;計算期望和(皮爾遜、標準化、調整的標準化)殘差;將缺失值作為一種有效值;進行行和列標題的標註;生成SAS或SPSS風格的輸出。參閱help(CrossTable)以了解詳情。
當有兩個以上的類別變數時，就需要生成多維列聯表，table() 和 xtabs() 都 可 以 基 於 三 個 或 更 多 的 類 別 型 變 量 生 成 多 維 列 聯 表 。表1中其它函數也都可以依次推廣到多維的情形（考慮篇幅有限，程式碼見文末客服二維碼）。
4、連續型變數的相關性檢驗
連續型變數中的相關性可以用相關係數來描述，相關係數的符號(±)表明關係的方向(正相 關或負相關)，其值的大小表示關係的強弱程度(完全不相關時為0，完全相關時為1)。（示例數據來自於R基礎安裝中的state.x77數據集）。R可以計算多種相關係數，包括Pearson相關係數、Spearman相關係數、Kendall相關係數、偏相關係數、多分格(polychoric)相關係數和多系列(polyserial)相關係數。cor()函數可以計算Pearson、Spearman、Kendall這三種相關係數，而cov()函數可用來計算協方差。這兩個函數中的use參數用來指定處理缺失數據的方式，method參數用來指定相關係數的類型。
圖6:函數cov()和cor()舉例
從上面的例子可以看出對角線上的相關係數都是1（變數和自身的相關性顯然為1），收入（income）和高中畢業率（HS Grad）相關性很強（0.61993232）。
而偏相關是指在控制一個或多個定量變數時，另外兩個定量變數之間的相互關係。你可以使用 ggm包中的pcor()函數計算偏相關係數。函數pcor()的參數為一個數值向量，前兩個數值表示要計算相關係數的變數下標，其餘的數值為條件變數(即要排除影響的變數)的下標，參數S為變數的協方差陣。
圖7，偏相關係數計算。
最後，polycor包中的hetcor()函數可以計算一種混合的相關矩陣，其中包括數值型變數的Pearson積差相關係數、數值型變數和有序變數之間的多系列相關係數、有序變數之間的多分格相關係數以及二分變數之間的四分相關係數。多系列、多分格和四分相關係數都假設有序變數或二分變數由潛在的正態分布導出。請參考此程式包所附文檔以了解更多。
在計算好相關係數以後，如何對它們進行統計顯著性檢驗呢? 函數cor.test()可以對單個的Pearson、Spearman和Kendall相關係數進行檢驗。其中的參數x和參數y為要檢驗相關性的變數，參數alternative則用來指定進行雙側檢驗或單側檢驗(取值為"two.side"、"less"或"greater")，而參數method用以指定要計算的相關類型("pearson"、" kendall"或"spearman" )。當研究的假設為總體的相關係數小於0時，使用alternative="less"。在研究的假設為總體的相關係數大於0時，應使用alternative="greater"。在默認情況下，假設為alternative="two.side"(總體相關係數不等於0)。cor.test()每次只能檢驗一種相關關係。但幸運的是，psych包中提供的corr.test()函數可以一次做更多事情，並且用法類似。psych包中的pcor.test()函數可以用於偏相關性係數的顯著性檢驗。另外，psych包中的r.test()函數提供了多種實用的顯著性檢驗方法。
圖8，corr.test()示例
5、分類變數的相關性檢驗
列聯表可以告訴你組成表格的各種變數組合的頻數或比例，不過你可能還會對列聯表中的變數是否相關或獨立感興趣。R提供了多種檢驗類別型變數獨立性的方法，接下來給大家介紹的三種檢驗分別為卡方獨立性檢驗、 Fisher精確檢驗和Cochran-Mantel-Haenszel檢驗。
圖6是用chisq.test()對示例數據做的卡方檢驗示例，說明了治療效果和性別是否獨立。但是下面的warning message是怎麼回事呢？因為在表中一個有一個小於5的值， 這可能會使卡方近似無效。
圖9:卡方檢驗示例。
可以使用fisher.test()函數進行Fisher精確檢驗來解決卡方檢驗無效的問題。
圖10，fisher.test()示例。
那麼這裡治療效果和性別是否獨立呢？mantelhaen.test()函數可用來進行Cochran-Mantel-Haenszel卡方檢驗，其原假設是，兩個名義變數在第三個變數的每一層中都是條件獨立的。用法和之前的兩個函數完全類似。
圖11:mantelhaen.test()示例。
從上面的獨立性檢驗結果可以看出我們關注的變數之間並不獨立，那自然可以考慮檢查變數之間的相關性。對於類別型變數，vcd包中提供了函數assocstats()用來計算二維列聯表的phi係數、列聯繫數和Cramer』s V係數（由於用法與前面三個函數太類似，不再贅述）。
6、連續型變數的比較檢驗
變數之間的關係除了獨立性、相關性之外，還可以進行比較，對於符合正態分布的連續型變數組間比較，我們一般採用t檢驗（示例數據為MASS包中的UScrime數據集）。T檢驗的函數為t.test()，有兩種調用格式：
- t.test(y~x, data) 其中的y是一個數值型變數，x是一個二分變數。
- t.test(y1, y2) 其中的y1和y2為數值型向量(即各組的結果變數)。
從下面的圖12，大家能不能看出group 0和group 1比較的結論呢？
圖12，t檢驗示例
函數t.test()也可以利用參數進行有方向的檢驗，不妨查看一下幫助文檔。
上面的例子是對於兩組獨立樣本的t檢驗，如果是非獨立樣本，將函數t.test()的參數paired設置為TRUE即可。如果是多於兩組的比較，需要用到方差分析，我們下一次再討論這部分內容。
當變數明顯不符合t檢驗或者方差分析的條件（比如非正態分布或者呈現有序關係），我們可以用非參數檢驗。若兩組數據獨立，可以使用Wilcoxon秩和檢驗。函數wilcox.test()就派上用場了。用非參數檢驗重複一下前面圖12中的比較。該函數調用方式與t.test()類似。這次的結論是否和圖12的結論一致呢？
圖13，非參數檢驗示例
有沒有多於兩組的非參數檢驗方法呢，答案是肯定的，如果各組獨立，則Kruskal-Wallis檢驗（函數kruskal.test()）將是一種實用的方法。如果各組不獨立(如重複測量設計或隨機區組設計)，那麼Friedman檢驗（函數friedman.test()）會更合適。（示例數據來自於R基礎安裝中的state.x77數據集。）
圖14，多組間的非參數檢驗。
如果要回答多組間的差異的細節，這裡有一個來自於網路的文件wmc.txt，包含了函數wmc(),可以解決這個難題（已經下載好放在後台）。結果如下。
圖15，函數wmc()示例
小結
這次的課程內容可以說是目前整個《R語言從入門到精通》系列課程中內容最多的一篇，而且涉及統計，理解上難度也比較大。希望大家能夠認真看，認真練，不要辜負老師的一番苦心啊~ 還是那句話：編程是不會把電腦編壞的，不要把電腦當作嬌花嫩草，多上手多練習才能記憶深刻。還在堅持學習的各位，要加油哦～
線下課程熱烈報名中，點擊查看詳情
本期乾貨
！R語言統計入門程式碼大全 ！
文章詳情：「科研貓」公眾號
科研貓原創系列，未經許可嚴禁轉載，版權事宜由上海辰明律師事務所提供法務支援。
