7.2 頻數表和列聯表
> library(vcd)
> head(Arthritis)
ID Treatment Sex Age Improved
1 57 Treated Male 27 Some
2 46 Treated Male 29 None
3 77 Treated Male 30 None
4 17 Treated Male 32 Marked
5 36 Treated Male 46 Marked
6 23 Treated Male 58 Marked
7.2.1 生成頻數表
函 數 描 述
table(var1, var2, …, varN) 使用 N 個類別型變量(因子)創建一個 N 維列聯表
xtabs(formula, data) 根據一個公式和一個矩陣或數據框創建一個 N 維列聯表
prop.table(table, margins) 依margins定義的邊際列表將表中條目表示為分數形式
margin.table(table, margins) 依margins定義的邊際列表計算表中條目的和
addmargins(table, margins) 將概述邊margins(默認是求和結果)放入表中
ftable(table) 創建一個緊湊的“平鋪”式列聯表
-
一維列聯表
> mytable<-with(Arthritis,table(Improved))
> mytable
Improved
None Some Marked
42 14 28
可以用prop.table()將這些頻數轉化為比例值:
> prop.table(mytable)
Improved
None Some Marked
0.5000000 0.1666667 0.3333333
或使用prop.table()*100轉化為百分比:
2. 二維列聯表
對於二維列聯表,table()函數的使用格式為:mytale<-table(A,B)
其中的A是行變量,B是列變量。除此之外,xtabs()函數還可使用公式風格的輸入創建列聯表,
格式為:mytable<-xtabs(~A+B,data=mydata)
其中的mydata是一個矩陣或數據框。總的來說,要進行交叉分類的變量應出現在公式的右側(即~符號的右方),以+作為分隔符。若某個變量寫在公式的左側,則其為一個頻數向量(在數據已經被表格化時很有用)。
對於Arthritis數據,有:
> mytable<-xtabs(~Treatment+Improved,data=Arthritis)
> mytable
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
可以使用margin.table()和prop.table()函數分別生成邊際頻數和比例。行和與行比
例可以這樣計算:
> margin.table(mytable,1)
Treatment
Placebo Treated
43 41
> prop.table(mytable,1)
Improved
Treatment None Some Marked
Placebo 0.6744186 0.1627907 0.1627907
Treated 0.3170732 0.1707317 0.5121951
列和與列比例可以這樣計算:
> margin.table(mytable,2)
Improved
None Some Marked
42 14 28
> prop.table(mytable,2)
Improved
Treatment None Some Marked
Placebo 0.6904762 0.5000000 0.2500000
Treated 0.3095238 0.5000000 0.7500000
各單元格所占比例可用如下語句獲取:
> prop.table(mytable)
Improved
Treatment None Some Marked
Placebo 0.34523810 0.08333333 0.08333333
Treated 0.15476190 0.08333333 0.25000000
可以使用addmargins()函數為這些表格添加邊際和
> addmargins(mytable)
Improved
Treatment None Some Marked Sum
Placebo 29 7 7 43
Treated 13 7 21 41
Sum 42 14 28 84
> addmargins(prop.table(mytable))
Improved
Treatment None Some Marked Sum
Placebo 0.34523810 0.08333333 0.08333333 0.51190476
Treated 0.15476190 0.08333333 0.25000000 0.48809524
Sum 0.50000000 0.16666667 0.33333333 1.00000000
在使用addmargins()時,默認行為是為表中所有的變量創建邊際和
> addmargins(prop.table(mytable,1),2)#僅添加了各行的和
Improved
Treatment None Some Marked Sum
Placebo 0.6744186 0.1627907 0.1627907 1.0000000
Treated 0.3170732 0.1707317 0.5121951 1.0000000
注意 table()函數默認忽略缺失值(NA)。要在頻數統計中將NA視為一個有效的類別,請設定參數useNA="ifany"。.
使用gmodels包中的CrossTable()函數是創建二維列聯表的第三種方法。CrossTable()
函數仿照SAS中PROC FREQ或SPSS中CROSSTABS的形式生成二維列聯表
> CrossTable(Arthritis$Treatment,Arthritis$Improved)
CrossTable()函數有很多選項,可以做許多事情:計算(行、列、單元格)的百分比;指
定小數位數;進行卡方、Fisher和McNemar獨立性檢驗;計算期望和(皮爾遜、標准化、調整的
標准化)殘差;將缺失值作為一種有效值;進行行和列標題的標注;生成SAS或SPSS風格的輸出。
3.多維列聯表
table()和xtabs()都可以基於三個或更多的類別型變量生成多維列聯margin.table()、prop.table()和addmargins()函數可以自然地推廣到高於二維的情況。另外,ftable()函數可以以一種緊湊而吸引人的方式輸出多維列聯表
> mytable<-xtabs(~Treatment+Sex+Improved,data=Arthritis)
, , Improved = None
Sex
Treatment Female Male
Placebo 19 10
Treated 6 7
, , Improved = Some
Sex
Treatment Female Male
Placebo 7 0
Treated 5 2
, , Improved = Marked
Sex
Treatment Female Male
Placebo 6 1
Treated 16 5
> ftable(mytable)
Improved None Some Marked
Treatment Sex
Placebo Female 19 7 6
Male 10 0 1
Treated Female 6 5 16
Male 7 2 5
> margin.table(mytable,c(1,3))#治療情況(Treatment) × 改善情況(Improved)的邊際頻數
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
7.2.2獨立性檢驗
1. 卡方獨立性檢驗
可以使用chisq.test()函數對二維表的行變量和列變量進行卡方獨立性檢驗
> library(vcd)
> mytable<-xtabs(~Treatment+Improved,data=Arthritis)
> chisq.test(mytable)
Pearson's Chi-squared test
data: mytable
X-squared = 13.055, df = 2, p-value = 0.001463#治療情況和改善情況不獨立
2. Fisher精確檢驗
可以使用fisher.test()函數進行Fisher精確檢驗。Fisher精確檢驗的原假設是:邊界固定
的列聯表中行和列是相互獨立的。其調用格式為fisher.test(mytable),其中的mytable是
一個二維列聯表
> fisher.test(mytable)
Fisher's Exact Test for Count Data
data: mytable
p-value = 0.001393
alternative hypothesis: two.sided
這里的fisher.test()函數可以在任意行列數大於等於2的二維列聯表上使用,但不能用於2×2的列聯表。
3.Cochran-Mantel—Haenszel檢驗
mantelhaen.test()函數可用來進行Cochran—Mantel—Haenszel卡方檢驗,其原假設是,兩
個名義變量在第三個變量的每一層中都是條件獨立的。
> mantelhaen.test(mytable)
Cochran-Mantel-Haenszel test
data: mytable
Cochran-Mantel-Haenszel M^2 = 14.6323, df = 2,
p-value = 0.0006647
7.2.3 相關性的度量
如果可以拒絕原假設,那么你的興趣就會自然而然地轉向用以衡量相關性強弱的相關性度量。vcd包中的assocstats()函數可以用來計算二維列聯表的phi系數、列聯系數和Cramer’sV系數
> mytable<-xtabs(~Treatment+Improved,data=Arthritis)
> assocstats(mytable)
X^2 df P(> X^2)
Likelihood Ratio 13.530 2 0.0011536
Pearson 13.055 2 0.0014626
Phi-Coefficient : 0.394
Contingency Coeff.: 0.367
Cramer's V : 0.394
總體來說,較大的值意味着較強的相關性。vcd包也提供了一個kappa()函數,可以計算混
淆矩陣的Cohen’s kappa值以及加權的kappa值。(舉例來說,混淆矩陣可以表示兩位評判者對於一系列對象進行分類所得結果的一致程度。)
7.2.5將表轉換為扁平格式
通過table2flat將表轉換為扁平格式
> table2flat<-function(mytable){
+ df<-as.data.frame(mytable)
+ rows<-dim(df)[1]
+ cols<-dim(df)[2]
+ x<-NULL
+ for(i in 1:rows){
+ for(j in 1:df$Freq[i]){
+ row<-df[i,c(1:(cols-1))]
+ x<-rbind(x,row)
+ }
+ }
+ row.names(x)<-c(1:dim(x)[1])
+ return(x)
+ }
使用table2flat()函數轉換已發表的數據
> treatment<-rep(c("Placebo","Treated"),times=3)
> improved<-rep(c("None","Some","marked"),each=2)
> Freq<-c(29,13,7,17,7,21)
> mytable<-as.data.frame(cbind(treatment,improved,Freq))
> mydata<-table2flat(mytable)
> head(mydata)
treatment inmproved
1 Placebo None
2 Placebo None
3 Placebo None
4 Placebo None
5 Treated None
6 Placebo Some
