频数表与列联表

本次训练关注类别型变量的频数表和列联表，以及相应的独立性检验、相关性的度量、图形化展示结果的方法，本次训练会用到vcd包和gmodels包中的函数，首次使用需要用install.packages进行安装。

数据集

本次训练的数据集来自vcd包中的Arthritis数据集，这份数据来自Kock & Edward（1988），表示了一项风湿性关节炎新疗法的双盲临床实验的结果，以下查看数据集的一些数据信息。 主要治疗手段（安慰剂治疗、用药治疗），性别（男性、女性）和改善情况（无改善、一定程度改善、显著改善）

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

常用的生成频数表和列联表的函数

函数 描述
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) 创建一个紧凑的“平铺”式列联表

1、一维列联表

使用table()函数生成简单的频数统计表 table函数默认是忽略缺失值（NA）的，要在频数统计中将NA视为一个有效的类别，可以设定参数useNA=“ifany”

(mytable <- with(Arthritis,table(Improved,useNA = "ifany")))
## 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转化为百分比

prop.table(mytable)*100
## Improved
##     None     Some   Marked 
## 50.00000 16.66667 33.33333

2、二维列联表

通过两个类别变量构建二维列联表 table函数的调用格式为table(A,B)，A是行变量，B是列变量

我们常用xtabs函数来构建二维列联表，调用格式如下： xtabs(* ~ A + B , data = mydata *)

(mytable2 <- xtabs(~Treatment+Improved,data = Arthritis))
##          Improved
## Treatment None Some Marked
##   Placebo   29    7      7
##   Treated   13    7     21

使用margin.table()和prop.table()函数分别生成边际频数和比例，计算方式如下：

#以行变量为维度进行计算行和行比例,"1"代表table()或xtabs()语句中的第一个变量
margin.table(mytable2,1)
## Treatment
## Placebo Treated 
##      43      41
prop.table(mytable2,1)
##          Improved
## Treatment      None      Some    Marked
##   Placebo 0.6744186 0.1627907 0.1627907
##   Treated 0.3170732 0.1707317 0.5121951
#以列变量为维度进行计算列和列比例,"2"代表table()或xtabs()语句中的第二个变量，多变量同理
margin.table(mytable2,2)
## Improved
##   None   Some Marked 
##     42     14     28
prop.table(mytable2,2)
##          Improved
## Treatment      None      Some    Marked
##   Placebo 0.6904762 0.5000000 0.2500000
##   Treated 0.3095238 0.5000000 0.7500000

使用addmargins()函数为上述列联表添加边际和：

#直接添加频数表的边际和
addmargins(mytable2)
##          Improved
## Treatment None Some Marked Sum
##   Placebo   29    7      7  43
##   Treated   13    7     21  41
##   Sum       42   14     28  84
#与prop.table组合只用，添加列联表比例的边际和
addmargins(prop.table(mytable2))
##          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(prop.table(mytable2,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
#按照第二个维度计算列比例，并添加每列的边际和
addmargins(prop.table(mytable2,2),1)
##          Improved
## Treatment      None      Some    Marked
##   Placebo 0.6904762 0.5000000 0.2500000
##   Treated 0.3095238 0.5000000 0.7500000
##   Sum     1.0000000 1.0000000 1.0000000

使用gmodels包中的CrossTable函数创建二维列联表（第三种方式） CrossTable函数的功能很强大，可以一次制定要计算的项目，如行、列、单元格的百分比，制定小数位数，进行卡方、Fisher和McNemar独立性检验；计算期望和残差；将缺失值作为一种有效值；进行行和列标题的标注；生成SAS和SPSS风格的输出。 help(“CrossTable”)

CrossTable(x, y, digits=3, max.width = 5, expected=FALSE, prop.r=TRUE, prop.c=TRUE, prop.t=TRUE, prop.chisq=TRUE, chisq = FALSE, fisher=FALSE, mcnemar=FALSE, resid=FALSE, sresid=FALSE, asresid=FALSE, missing.include=FALSE, format=c(“SAS”,“SPSS”), dnn = NULL, …)

因此，我们可以用过设定chisq和expected来要求函数进行卡方检验两个变量的独立性以及求出期望值。

library(gmodels)
CrossTable(Arthritis$Treatment,Arthritis$Improved,chisq = TRUE,expected=TRUE)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |              Expected N |
## | Chi-square contribution |
## |           N / Row Total |
## |           N / Col Total |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  84 
## 
##  
##                     | Arthritis$Improved 
## Arthritis$Treatment |      None |      Some |    Marked | Row Total | 
## --------------------|-----------|-----------|-----------|-----------|
##             Placebo |        29 |         7 |         7 |        43 | 
##                     |    21.500 |     7.167 |    14.333 |           | 
##                     |     2.616 |     0.004 |     3.752 |           | 
##                     |     0.674 |     0.163 |     0.163 |     0.512 | 
##                     |     0.690 |     0.500 |     0.250 |           | 
##                     |     0.345 |     0.083 |     0.083 |           | 
## --------------------|-----------|-----------|-----------|-----------|
##             Treated |        13 |         7 |        21 |        41 | 
##                     |    20.500 |     6.833 |    13.667 |           | 
##                     |     2.744 |     0.004 |     3.935 |           | 
##                     |     0.317 |     0.171 |     0.512 |     0.488 | 
##                     |     0.310 |     0.500 |     0.750 |           | 
##                     |     0.155 |     0.083 |     0.250 |           | 
## --------------------|-----------|-----------|-----------|-----------|
##        Column Total |        42 |        14 |        28 |        84 | 
##                     |     0.500 |     0.167 |     0.333 |           | 
## --------------------|-----------|-----------|-----------|-----------|
## 
##  
## Statistics for All Table Factors
## 
## 
## Pearson's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 =  13.05502     d.f. =  2     p =  0.001462643 
## 
## 
##

通过结果我们可以看到卡方检验的自由度为2，卡方检验统计值为13.06，p值为0.0015，显著＜0.05。因此可以拒绝原假设（H0：Treatment治疗手段与Improved改善结果无关），得出结论治疗手段与改善结果有着密切的关系。（假设检验会在下次训练中详细操作，这里作为一个方法的延展）

3、多维列联表

table()函数和xtabs()函数均可构建三维以上的多维列联表，同理上述使用的margin.table，prop.table，addmargins函数都可以用于三维以上的多维列联表。而ftable可以以一种紧凑的方式输出多维列联表。如下：

(mytable3 <- 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(mytable3)
##                  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函数作用于多维列联表

margin.table(mytable3,1)
## Treatment
## Placebo Treated 
##      43      41
margin.table(mytable3,2)
## Sex
## Female   Male 
##     59     25
margin.table(mytable3,3)
## Improved
##   None   Some Marked 
##     42     14     28
margin.table(mytable3,c(1,3))#个人也理解为降维计算处理
##          Improved
## Treatment None Some Marked
##   Placebo   29    7      7
##   Treated   13    7     21

将prop.table作用于多维列联表

#ftable 让结果比较直观
ftable(prop.table(mytable3,c(1,2)))#按照治疗手段和性别的维度进行统计
##                  Improved       None       Some     Marked
## Treatment Sex                                             
## Placebo   Female          0.59375000 0.21875000 0.18750000
##           Male            0.90909091 0.00000000 0.09090909
## Treated   Female          0.22222222 0.18518519 0.59259259
##           Male            0.50000000 0.14285714 0.35714286
ftable(addmargins(prop.table(mytable3,c(1,2)),3))#按照第三个维度进行求和
##                  Improved       None       Some     Marked        Sum
## Treatment Sex                                                        
## Placebo   Female          0.59375000 0.21875000 0.18750000 1.00000000
##           Male            0.90909091 0.00000000 0.09090909 1.00000000
## Treated   Female          0.22222222 0.18518519 0.59259259 1.00000000
##           Male            0.50000000 0.14285714 0.35714286 1.00000000
ftable(addmargins(prop.table(mytable3,c(1,2)),3))*100#给出百分比
##                  Improved       None       Some     Marked        Sum
## Treatment Sex                                                        
## Placebo   Female           59.375000  21.875000  18.750000 100.000000
##           Male             90.909091   0.000000   9.090909 100.000000
## Treated   Female           22.222222  18.518519  59.259259 100.000000
##           Male             50.000000  14.285714  35.714286 100.000000

到此，就是列联表操作的全部训练内容，尽管我们能看出看出变量间的差异，但我们无法直观看出各种变量间是否相关或独立。如果不独立，那相关性的强度如何？这些问题需要我们进一步探讨的。

独立性检验

R提供了多种检验类别型变量的独立性方法，此处训练将继承上述的结果进行独立性检验，这里训练三种检验方法：

1、卡方独立检验

R中可以非常容易的使用chisq.test()函数对二维列联表的行变量和列变量进行卡方独立性检验。 我们检验患者接受治疗和改善水平之间的独立性 方式如下：

library(vcd)
#mytable2为上述二维列联表训练用xtabs函数创建
(csq <- chisq.test(mytable2))
## 
##  Pearson's Chi-squared test
## 
## data:  mytable2
## X-squared = 13.055, df = 2, p-value = 0.001463
#检验的结果可以通过csq获取其他你关心的结果
csq$expected
##          Improved
## Treatment None     Some   Marked
##   Placebo 21.5 7.166667 14.33333
##   Treated 20.5 6.833333 13.66667

这里的结果与前面我们用CrossTable函数创建二维列联表时，使用参数chisq=TRUE所显示的检验结果一致，R中有多种方法实现卡方检验。这里是最简单的一种。 结果表明患者接受治疗和改善水平看上去存在着某种关系（p<0.01）。 下面是解释过程：

  1. 原假设H0：患者采用的治疗手段与改善的结果之间没有关系；备择假设H1：患者采用的治疗手段与改善的结果之间存在某种关系
  2. 选择1%的显著水平（常用5%，显著水平越低，检验效果越好，得出结论的力度越大）
  3. p值<0.01,则可以拒绝原假设H0，从而接受备择假设H1
  4. 得出结论，患者采用的治疗手段与改善的结果之间存在某种关系

我们再检验改善效果与性别之间的独立性

#为了方便调用前面的结果，这里利用了margin.table函数可以降维的特性，构建了性别和改善结果的二维列联表（或者称为边际表）
chisq.test(margin.table(mytable3,c(2,3)))
## Warning in chisq.test(margin.table(mytable3, c(2, 3))): Chi-squared
## approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  margin.table(mytable3, c(2, 3))
## X-squared = 4.8407, df = 2, p-value = 0.08889

这里的结果表明患者性别和改善情况之间不存在关系或者说相互独立（p>0.05）

2、Fisher精确检验

可以使用fisher.test()函数进行Fisher精确检验。 Fisher精确检验的原假设是H0：边界固定的列联表中行和列是相互独立的。使用方法如下

fisher.test(mytable2)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  mytable2
## p-value = 0.001393
## alternative hypothesis: two.sided

结果表明拒绝原假设

3、Cochran-Mantel-Haenszel检验

可以使用mantelhean.test()函数进行Cochran-Mantel-Haenszel卡方检验，其原假设是H0：两个名义变量在第三个变量的每一层中都是条件独立的。 例如：我们要检验治疗手段和改善情况在性别的每一水平下是否独立。

mytable4 <- xtabs(~Treatment+Improved+Sex,data = Arthritis)
mantelhaen.test(mytable4)
## 
##  Cochran-Mantel-Haenszel test
## 
## data:  mytable4
## Cochran-Mantel-Haenszel M^2 = 14.632, df = 2, p-value = 0.0006647

结果表明拒绝原假设

相关性的度量

以上的结果均表明我们所检验的变量有某种相关性，我们可以用assocstats()函数来计算二维列联表的phi系数、列联系数和Cramer’s V系数

library(vcd)
#治疗手段和改善情况的二维列联表
assocstats(mytable2)
##                     X^2 df  P(> X^2)
## Likelihood Ratio 13.530  2 0.0011536
## Pearson          13.055  2 0.0014626
## 
## Phi-Coefficient   : NA 
## Contingency Coeff.: 0.367 
## Cramer's V        : 0.394

可视化

往后统一训练