R-无序的定类数据分析：列联表、热力图、和弦图、桑基图和统计检验

2022-05-08 本文已影响0人睿睿爱汪汪也爱喵喵

今天我们将通过一个例子来说明如何分析两个定类变量。

文章背景：我们想研究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)