今天我们将通过一个例子来说明如何分析两个定类变量。
文章背景：我们想研究CFPS2010和CFPS2012青少年对自身的职业期望。
如表1，我们将原始的职业期望编码整合成9类（职业编码的大类）和其他。由于我们想分析同一个人在跨轮次调查中职业期望的稳定性情况，故将分析对象定义为在CFPS2010和CFPS2012中都回答了自己对自己职业期望的受访者。如表2所示，进行重编码后的数据是宽数据，样本量是1920，数据集名字为expect。我们在进行后续分析时，要将其转换为绘图所需的其他形式。
表1 职业期望编码的重编码对应关系示意表
表2 CFPS2010和CFPS2012青少年职业期望的情况(数据集expect）
⭐分析方式1——列联表、频数与频率
在表3中，我们展示了2010与2012年青少年职业期望的交叉统计情况。同时该表内，也附上了频数（落在各类别中的数据个数）、⽐例（某⼀类别数据占全部数据的⽐值）、百分⽐（将对⽐的基数作为100⽽计算的⽐值，包括百分比、行百分比和列百分比）。
表3 CFPS2010与CFPS2012青少年职业期望的频数分析
⭐分析方式2——统计图表
分析前色彩讲解：预设渐变色，我们这里介绍2个色彩包。
1）专门生产系列颜色的RColorBrewer包，详见图1中的系列颜色。
library(RColorBrewer)
display.brewer.all()
图1 RColorBrewer包系列颜色展示
2）色盲友好的配色方案viridis包，详见图2中的系列颜色。
library(viridis)
?viridis()#可以看到更多对这组包色彩的说明
图2 viridis包系列颜色展示
接下来我们来画图吧~【注：图3-图6中的类目数字的含义：1）国家机关、党群组织、企业、事业单位负责人；2）医生；3）教师；4）专业技术人员（刨除教师和医生）；5）办事人员和有关人员；6）商业、服务业人员；7）农、林、牧、渔、水利业生产人员；8）生产、运输设备操作人员及有关人员；9）军人；10）其他。】
😀绘图1：相邻矩阵的热力图【绘图工具：ggplot2或专业的ComplexHeatmap包】
相邻矩阵是指代表N个节点之间关系的N*N矩阵，矩阵内的位置（i,j）表示第i和j个节点之间的关系。对本文的有向关系网络来说，相邻矩阵不具有对称性。相邻矩阵的的对角线表示节点与自己的关系情况。
相邻矩阵可以用热力图来表示，将节点之间连接的权重用颜色来表达。
Method 1：用ggplot2绘制
library(ggplot2)
library(reshape2)
expect_Heat<-as.data.frame(round(prop.table(table(expect$code_new10,expect$code_new12)),5)*100)
colnames(expect_Heat)<-c("from","to","value")
#expect_Heat（表4）就是我们生成的用ggplot2来画热力图的数据集。
表4 用ggplot2画热力图所需要的数据集
ggplot(expect_Heat,aes(x=to,y=from,fill=value,label=value))+
geom_tile(color="black")+
scale_fill_gradientn(colors=brewer.pal(9,"YlGnBu"))+
xlab('CFPS2012')+
ylab('CFPS2010')+
coord_equal()+
theme(axis.text.x=element_text(angle=90,hjust=1,colour='black'),
axis.text.y=element_text(angle=0,hjust=1,colour='black'))
图3 2010和2012年青少年对自己职业期望的相关性热力图 （用 ggplot2包画的）
#####安装ComplexHeatmap前需先调入devtools包，然后我们再通过install_github函数安装ComplexHeatmap包
library(devtools)
install_github("jokergoo/ComplexHeatmap")
library(ComplexHeatmap)
myColors=brewer.pal(9,"YlGnBu")[1:9]
expect_Heat_matrix<-as.matrix(round(prop.table(table(expect$code_new10,expect$code_new12)),5)*100)
#expect_Heat_matrix（图4）就是我们生成的用ComplexHeatmap包绘制热力图所需要的矩阵。
Heatmap(expect_Heat_matrix,myColors)
图 4 用ComplexHeatmap包绘图所需要的矩阵
图5 2010和2012年青少年对自己职业期望的相关性热力图（用ComplexHeatmap包画的）
😀绘图2：和弦图【绘图工具：circlize包的chordDiagram函数】
和弦图(chord Diagram)是一种显示矩阵中数据之间相互关系的数据可视化方法，主要用来展示多个对象之间的关系。
和弦图主要由节点和弦构成，节点数据沿圆周径向排列，连接圆上任意两点的带权重（有宽度）的弧线称之为弦，弦（两点之间的连线）就代表着两者之间的关联关系。不同节点和弦之间用颜色将数据加以分类，能够直观得对数据进行比较和区分，十分适合用来表示复杂数据之间的关联关系。
图6表现了2010至2012年青少年职业期望的演变特征，节点数为10表示要分析十种职业期望，数量通过弧线的颜色和宽度使其相互间的关系表达清晰：不同颜色的连线表示不同职业期望之间的关联程度，权重（弧线宽度）可发现不同期望之间选择上的明显变化。
library(circlize)
library(viridis)
expect_chord<-expect_heat
chordDiagram(expect_chord,col=viridis::turbo(10))
图6 2010和2012年青少年对自己职业期望的和弦图
😀绘图3：桑基图【绘图工具：ggalluvial包】
桑基图有利于展现分类数据之间的相关性，以流的形式呈现同一类别的元素数量。在图7中，边表示职业期望之间的流动情况，流量表示流动数据的具体数值，节点表示不同的职业期望分类。
通过图7，我们可以总结一下桑基图的特点：
1）起始总流量=结束总流量，即能量守恒。不能在中间过程中创造出流量，也不能有流量总量上的耗损。故，假如出现耗损或新增，可以设一个新的类别，如“其他”，存储这些不属于分析对象的类目。
2）分析流量的流动过程，不同线条代表不同的流量分流情况，边的宽度与流量成比例，边越宽，流量数值越大，
#install.packages('ggalluvial')
library(ggalluvial)
#去掉pid列，对expect数据集进行格式转换，转换后生成的用来画桑基图的数据集如表5，其中person列在绘图时代表流量，year列区分了流动中的状态，expectancy代表节点。
sapply(names(expect),function(x) length(unique(expect[,x])))
expect_sanky<-to_lodes_form(co_2[,1:ncol(expect)],axes=1:ncol(expect),id="person")
colnames(expect_sanky) <- c("person","year","expectancy")
expect_sanky<-mutate(expect_sanky,year=ifelse(year=="code_new10","2010","2012"))
表5 绘制桑基图需要的数据集展示
ggplot(expect_sanky,
aes(x = year, stratum = expectancy ,alluvium=person,
fill = expectancy, label = expectancy)) +
geom_flow() +
geom_stratum(alpha = .5) +
geom_text(stat = "stratum", size = 3) +
theme(legend.position = "none") +
ggtitle("job expectancy at two points in time")+
theme(plot.title = element_text(hjust = 0.5))+
theme(panel.grid.major=element_line(colour=NA),
panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA),
panel.grid.minor = element_blank())
图7 2010-2012青少年对自己职业期望的变化的桑基图
注释：这三种图都可以进行美化，我在绘制部分图形时，并没有调参数，具体的内容修饰见参考文件。
⭐分析方式3——统计检验
😊方法一：卡方独立性检验：研究两组分类变量的关系
library(gmodels)
CrossTable(expect$code_new10, expect$code_new12,expected = T,format = "SAS",fisher = T,prop.c = T,prop.t = F,prop.chisq = T)
图8 CrossTable检验的结果
也可以用stats包中的统计检验分别计算。
1）Pearson 卡方独立性检验
chisq.test(expect$code_new10, expect$code_new12, correct = TRUE,
p = rep(1/length(x), length(x)), rescale.p = FALSE,
simulate.p.value = FALSE, B = 2000)
有多于20%的期望频数小于5，最小的期望频数T=0.025。依据卡方检验的使用条件，我们谨慎起见，选择读取Fisher's Exact Test。
2）Fisher精确检验
fisher.test(co_2$code_new10,co_2$code_new12, workspace = 5000, hybrid = FALSE,
control = list(), or = 1, alternative = "two.sided",simulate.p.value=TRUE,
conf.int = TRUE, conf.level = 0.95)
图9 Pearson 卡方独立性检验和 Fisher精确检验的结果
由P值<0.01，拒绝原假设，认为2010年和2012年青少年的职业期望是独立的。
😊方法二：一致性检验
Kappa系数用于一致性检验，也可以用于衡量分类精度，kappa系数的计算是基于混淆矩阵的。Kappa检验是用于检验两个变量的测量结果是否一致，强调一致性。从研究目的上来说，一致性、差异性和独立性是不同的概念。
kappa系数进行解读：一般认为kappa>0.8，则有很好的一致性，kappa<0.4，则一致性较差；此外还应通过显著性检验，也就是p值<0.05。
library(irr)
kappa2(expect)
kappa2(expect, "squared") # predefined set of squared weights
图10 kappa检验的结果
根据图10，无论是原数据，还是考虑平方后的数据，我们均可以得到结论：在0.01的显著性水平下，2010年和2012年青少年职业期望的一致性较差。
#kappam.fleiss()：在2010、2012的基础上，在添加新的列——2014、2016等，该函数可以检验多个变量的测量结果是否一致。Fleiss' kappa系数：取值范围是【-1,1】，其系数越大，一致性越强，一般经验认为低于0.7是不可接受的尚需要改进，0.7可接受，0.9优秀。
😊方法三：配对分类数据的差异性检验
分析配对分类数据的差异性。此类数据最常见于实验研究，⽤不同的⽅法或不同的时间点检测同⼀批⼈，看两个⽅法的效果是否有差异。此时可使⽤配对卡⽅检验。注意，行名和列名具有相同的分类，并且计数表示成对响应。本质上，那些前后反应相同的个体并不影响对反应变化的评估。我们主要研究“不一致”的计数。Kappa检验会利用列联表的全部数据，而配对卡方检验只利用“不一致“数据。
当数据维度是2*2时，用配对分类数据的McNemar 检验；当数据维度超过2*2时，用McNemar–Bowker检验。
注：1）McNemar–Bowker检验的一个缺点：如果矩阵中的某些位置有0，它可能会失败；
2）如果“不一致”分类的计数很低，McNemar 检验可能不可靠。
mcnemar.test(expect$code_new10,expect$code_new12)
图11 McNemar–Bowker检验的结果
根据图11，McNemar–Bowker检验失败了。
if(!require(rcompanion)){install.packages("rcompanion")}
还可以利用Rcompanion包的nominalSymmetryTest函数对数据进行对称性检验，对本数据也是失效的。
⭐分析方式3——统计模型
待更新。
参考文件
1、《R语言数据可视化之美：专业图表绘制指南》（增强版），张杰，中国工信出版集团&电子工业出版社。该书的电子版见微信读书。
2、 circlize官方教程。网址：R 中的循环可视化 (jokergoo.github.io)。
3、ComplexHeatmap官方教程。网址：Heatmap Complete Reference (jokergoo.github.io)
4、定类数据分析与定序数据分析。二、定类与定序变量分析 - 百度文库 (baidu.com)
5、定类数据与卡方检验。网址：定类数据如何分析？卡方检验有什么使用场景？ - 百度文库 (baidu.com)
6、R Handbook: Tests for Paired Nominal Data (rcompanion.org)
网友评论