1 分析前的环境设置

1.1 下载安装套件包

本章学习者須先下载安装相關套件包，方能正常運行。

#下载安装套件包，如果之前有下载过，会跳过不下载；如果之前没有下载过，会下载安装套件包
if (!require("showtext")) install.packages("showtext")
if (!require("plyr")) install.packages("plyr")
if (!require("dplyr")) install.packages("dplyr")
if (!require("sjlabelled")) install.packages("sjlabelled")
if (!require("sjmisc")) install.packages("sjmisc")
if (!require("sjPlot")) install.packages("sjPlot")
if (!require("descr")) install.packages("descr")

1.2 环境设置

我们通常会将R的工作环境先进行一些设置，这些设置包括了 “设定当前工作目录”、“设定系統中文文字编码”、以及“设定绘图物件中的中文文字”等项

#设置工作目录
setwd("/Users/simonfair/Desktop/量化研究数据分析/R Ch7 百分比交叉分析")

#查看目前工作目录
getwd()

[1] "/Users/simonfair/Desktop/量化研究数据分析/R Ch7 百分比交叉分析"

#设置系統中文文字编码(简体中文)
Sys.setlocale(category = "LC_ALL", locale = "zh_CN.UTF-8")

[1] "zh_CN.UTF-8/zh_CN.UTF-8/zh_CN.UTF-8/C/zh_CN.UTF-8/zh_CN.UTF-8"

#使绘图物件中的中文文字能正确呈现
#调用加载"showtext"套件包
library(showtext) 
showtext_auto(enable = TRUE)

2 数据预备

2.1 导入外部数据

#导入SPSS数据档(*.sav)
library(sjlabelled)
TCS2015sc <- read_spss("TCS2015sc.sav")

Converting atomic to factors. Please wait...

#存成R数据档(*.rda)
save(TCS2015sc, file = "TCS2015sc.rda")
#导入R数据档(*.rda)
load("TCS2015sc.rda")

3 R * C百分比交叉分析(Crosstabulate Analysis)

3.1 挑选分析变量 dplyr::select ( )

使用 dplyr 套件包 select( ) 函数，从R数据档(TCS2015sc)中挑选分析时所需要的变量（此处挑选5个变量），形成一个新的分析数据档(mydata)。用 dim( )函数查看数据档中的观察个案数与变量数，用names( )函数 查看数据档中的变量名称

library(dplyr)
mydata <- select(TCS2015sc,
  A1,A8,V1.5,age_strata, W_Raking)
dim(mydata) #2002行(观察值), 5列(变量)

[1] 2002    5

names(mydata) #查看数据档中的变量名称

[1] "A1"         "A8"         "V1.5"       "age_strata" "W_Raking"

3.2 呈现分析变量的频数分布 sjmisc::frq(, weights= )

本例，使用 sjmisc 套件包 frq( ) 函数，呈现分析自变量(A8：教育程度)，与因变量(V1.5, 日子过的快不快乐)的加权后频数分布（加权变量：W_Raking）。

library(sjmisc)
frq(mydata$A8,weights = mydata$W_Raking) #加权

A8. 教育程度 (xw) <categorical> 
# total N=2002 valid N=2002 mean=10.90 sd=6.53

Value |           Label |   N | Raw % | Valid % | Cum. %
--------------------------------------------------------
    1 |              无 |  66 |  3.30 |    3.30 |   3.30
    2 |            自修 |   5 |  0.25 |    0.25 |  45.25
    3 |            小学 | 220 | 10.99 |   10.99 |  62.84
    4 |            初中 | 210 | 10.49 |   10.49 |  73.33
    5 |            初职 |   9 |  0.45 |    0.45 |  73.78
    6 |      高中普通科 | 112 |  5.59 |    5.59 |  79.37
    7 |      高中职业科 | 112 |  5.59 |    5.59 |  84.97
    8 | 高职(高商/高工) | 300 | 14.99 |   14.99 |  99.95
    9 |        士官学校 |   1 |  0.05 |    0.05 | 100.00
   10 |            五专 |  74 |  3.70 |    3.70 |   6.99
   11 |            二专 | 144 |  7.19 |    7.19 |  14.19
   12 |            三專 |  23 |  1.15 |    1.15 |  15.33
   13 |      军警专修班 |   0 |  0.00 |    0.00 |  15.33
   14 |      军警专科班 |   4 |  0.20 |    0.20 |  15.53
   15 |        空中行专 |   2 |  0.10 |    0.10 |  15.63
   16 |        空中大学 |   7 |  0.35 |    0.35 |  15.98
   17 |  军警官校或大学 |   9 |  0.45 |    0.45 |  16.43
   18 |   技术学院/科大 | 130 |  6.49 |    6.49 |  22.93
   19 |            大学 | 442 | 22.08 |   22.08 |  45.00
   20 |            硕士 | 121 |  6.04 |    6.04 |  51.30
   21 |            博士 |  11 |  0.55 |    0.55 |  51.85
   88 |            其他 |   0 |  0.00 |    0.00 |  99.95
 <NA> |            <NA> |   0 |  0.00 |    <NA> |   <NA>

frq(mydata$V1.5,weights = mydata$W_Raking) #加权

V1.5. 整体来说，你觉得目前的日子过得快不快乐？ (xw) <categorical> 
# total N=2002 valid N=2002 mean=5.50 sd=15.00

Value |      Label |    N | Raw % | Valid % | Cum. %
----------------------------------------------------
    1 | 非常不快樂 |   22 |  1.10 |    1.10 |   1.10
    2 |     不快乐 |  232 | 11.59 |   11.59 |  12.69
    3 |       快乐 | 1490 | 74.43 |   74.43 |  87.11
    4 |   非常快乐 |  202 | 10.09 |   10.09 |  97.20
   94 |   无法选择 |   56 |  2.80 |    2.80 | 100.00
 <NA> |       <NA> |    0 |  0.00 |    <NA> |   <NA>

3.3 进行变量的编码转换

3.3.1 将因变量重新编码:V1.5(过的快不快乐)(5分类) -> V1.5.cat3(3分类)

加载 sjmisc 套件包，使用 rec( )函数，将原始变量(V1.5)依照我们的编码规则(1:2=1; 3:4=2; 94=3)，重新编码(rec)为另一新的变量(V1.5.cat3)；

使用 rec( )函数的 val.labels( ) 与 value.labels( ) 参数，赋予变量标签和变量值标签；

使用 rec( ) 函数的 as.num=FALSE 参数，将新的变量(V1.5.ca3)设置为类别变量；

最后，使用 sjmjsc 套件包的 flat_table( )函数，检查编码转换是否正确，并使用 frq( ) 函数，呈现重新编码后的频数分布。

library(sjmisc) #加载sjmisc套件
mydata$V1.5.cat3 <- rec(mydata$V1.5, #使用rec函数重新编码
                      rec =  "1:2=1; 3:4=2; 94=3", #设置编码规则
                      var.label = "日子过的快不快乐三分类",
                      val.labels = c("不快乐","快乐","无法选择"),
                      as.num = FALSE #设置为类别变量
                      )

flat_table(mydata, V1.5, V1.5.cat3)#检查编码转换是否正确

           V1.5.cat3 不快乐 快乐 无法选择
V1.5                                     
非常不快樂               21    0        0
不快乐                  232    0        0
快乐                      0 1498        0
非常快乐                  0  195        0
无法选择                  0    0       56

frq(mydata$V1.5.cat3,
  weights = mydata$W_Raking)#做V1.5.cat3的频数分布(加权)

日子过的快不快乐三分类 (xw) <categorical> 
# total N=2002 valid N=2002 mean=1.90 sd=0.38

Value |    Label |    N | Raw % | Valid % | Cum. %
--------------------------------------------------
    1 |   不快乐 |  254 | 12.69 |   12.69 |  12.69
    2 |     快乐 | 1692 | 84.52 |   84.52 |  97.20
    3 | 无法选择 |   56 |  2.80 |    2.80 | 100.00
 <NA> |     <NA> |    0 |  0.00 |    <NA> |   <NA>

3.3.2 将自变量重新编码:A8(教育程度)(21分类) -> edu.cat5(5分类)，

使用同样的处理程序，将自变量 A8 (教育程度)(21分类) 重新编码为 edu.cat5 (5分类)。

library(sjmisc)#加载sjmisc套件
mydata$edu.cat5 <- rec(mydata$A8,  #使用rec函数重新编码
    rec =  "1:3=1; 4:5=2; 6:9=3; 10:15=4; 16:21=5",#设置编码规则
    var.label = "受教育程度五分类",
    val.labels = c("小学及以下",
                   "初中",
                   "高中职", 
                   "专科",
                   "本科及以上"),
    as.num = FALSE #设定为类别变量
    )

flat_table(mydata, A8, edu.cat5)#检查编码转换是否正确

                edu.cat5 小学及以下 初中 高中职 专科 本科及以上
A8                                                             
无                               63    0      0    0          0
五专                              5    0      0    0          0
二专                            235    0      0    0          0
三專                              0  225      0    0          0
军警专修班                        0    9      0    0          0
军警专科班                        0    0    113    0          0
空中行专                          0    0    114    0          0
空中大学                          0    0    312    0          0
军警官校或大学                    0    0      1    0          0
技术学院/科大                     0    0      0   75          0
大学                              0    0      0  152          0
自修                              0    0      0   23          0
硕士                              0    0      0    4          0
博士                              0    0      0    2          0
小学                              0    0      0    0          7
初中                              0    0      0    0          9
初职                              0    0      0    0        113
高中普通科                        0    0      0    0        406
高中职业科                        0    0      0    0        123
高职(高商/高工)                   0    0      0    0         11
其他                              0    0      0    0          0
士官学校                          0    0      0    0          0

frq(mydata$edu.cat5,
    weights = mydata$W_Raking) #做edu.cat5的频数分布

受教育程度五分类 (xw) <categorical> 
# total N=2002 valid N=2002 mean=3.44 sd=1.43

Value |      Label |   N | Raw % | Valid % | Cum. %
---------------------------------------------------
    1 | 小学及以下 | 291 | 14.54 |   14.54 |  14.54
    2 |       初中 | 219 | 10.94 |   10.94 |  25.47
    3 |     高中职 | 526 | 26.27 |   26.27 |  51.75
    4 |       专科 | 247 | 12.34 |   12.34 |  64.09
    5 | 本科及以上 | 719 | 35.91 |   35.91 | 100.00
 <NA> |       <NA> |   0 |  0.00 |    <NA> |   <NA>

3.4 呈现交叉分析频数百分比表

在进行交叉分析时，每个交叉细格通常可以呈现三种百分比： row%, col%, 以及 total%。

学习者最困扰的是，不知要选择哪一种百分比作为比较的标准。一般的建议是：呈现自变量所在的百分比。

如果我们想将自变量放在表格的左边，就要呈现 row 百分比。

3.4.1 呈现方式(1)：自变量在row, 因变量在column

就本例来说，edu5.cat5(教育水平五分类) 是我们所设想的自变量。在进行交叉分析时，我们想将其放在表格的左边，所以就要呈现row%。我们可以调用 sjPlot 套件包的 tab_xtab ( ) 函数，呈现交叉分析的 row% (show.row.prc = TRUE)。

另外，由于我们使用的是抽样调查资料，该资料中亦有提供加权变量(W_Raking)供研究者使用，所以我们记得要使用加权变量进行分析(weight.by=mydata$W_Raking)，这样产生出来的数据，才会是一個可以代表调查总体的正确结果。

如 Table 1

library(sjPlot) #加载sjPlot套件
#使用tab_xtab函数做交叉分析表
tab_xtab(mydata$edu.cat5, mydata$V1.5.cat3, 
        #row自变量edu.cat5,col因变量V1.5.cat3
        weight.by = mydata$W_Raking,#加权
        show.row.prc = TRUE,#呈现 row %
        show.legend = TRUE,#呈现cell数据说明
        show.summary = FALSE,
        #单纯呈现交叉表,不呈现卡方检验相关数据
        encoding = "utf8") #(for win)防止报表出现中文乱码

Table 1: 交叉分析表示例(1)-自变量在row

受教育程度五分类	日子过的快不快乐三分类			Total
受教育程度五分类	不快乐	快乐	无法选择	Total
小学及以下	37 12.8 %	249 85.9 %	4 1.4 %	290 100 %
初中	42 19.2 %	171 78.1 %	6 2.7 %	219 100 %
高中职	69 13.1 %	444 84.4 %	13 2.5 %	526 100 %
专科	30 12.1 %	211 85.1 %	7 2.8 %	248 100 %
本科及以上	76 10.6 %	618 85.8 %	26 3.6 %	720 100 %
Total	254 12.7 %	1693 84.5 %	56 2.8 %	2003 100 %

observed values
% within 受教育程度五分类

3.4.2 呈现方式(2)：自变量在column, 因变量在row

另外一种呈现方式是将自变量放在表格的上头，此时就要呈现 column百分比。以本例来说，edu5.cat5是我们所设想的自变量，在进行交叉分析时，我们想将其放在表格的上头，所以就要呈现 col% (show.col.prc = TRUE)，一样要记得将加权变量(W_Raking)放入分析语法中。

如 Table 2

library(sjPlot)#加载sjPlot套件
#使用tab_xtab函数做交叉分析表
tab_xtab(mydata$V1.5.cat3, mydata$edu.cat5, 
      #row因变量V1.5.cat3, col自变量edu.cat5
                 weight.by = mydata$W_Raking,#加权
                 show.col.prc = TRUE,#show col %
                 show.legend = TRUE,#呈现cell数据说明
                 show.summary = FALSE, 
                 #单纯呈现交叉表,不呈现卡方检验相关数据
                 encoding = "utf8") 
                  #(for win)防止报表出现中文乱码

Table 2: 交叉分析表示例(1)-自变量在column

日子过的快不快乐三分类	受教育程度五分类					Total
日子过的快不快乐三分类	小学及以下	初中	高中职	专科	本科及以上	Total
不快乐	37 12.8 %	42 19.2 %	69 13.1 %	30 12.1 %	76 10.6 %	254 12.7 %
快乐	249 85.9 %	171 78.1 %	444 84.4 %	211 85.1 %	618 85.8 %	1693 84.5 %
无法选择	4 1.4 %	6 2.7 %	13 2.5 %	7 2.8 %	26 3.6 %	56 2.8 %
Total	290 100 %	219 100 %	526 100 %	248 100 %	720 100 %	2003 100 %

observed values
% within 受教育程度五分类

3.5 呈现交叉分析频数百分比图

在进行交叉分析时，我们也可以使用 sjPlot 套件包的 plot_xtab( ) 函数，呈现交叉分析百分比的条形图(type=“bar”)以及折线图(type=“line”)，以可视化的方式，查看两变量间的关系。

3.5.1 呈现方式(1)：交叉分析频数百分比条形图 sjPlot::plot_xtab (, type=“bar”)

自变量(edu.cat5)在x轴(row), 因变量(V1.5.cat3)在y轴(col)。

如 Figure 1

library(sjPlot)#加载sjPlot套件
#使用plot_xtab函数绘图
plot_xtab(mydata$edu.cat5, mydata$V1.5.cat3,
     #row:x轴变量(自变量)edu.cat5，col:y轴变量(因变量)V1.5.cat3
          type = "bar",  #呈现条形图 
          margin = "row", #呈现row% (自变量所在的%)
          weight.by = mydata$W_Raking)#加权

3.5.2 呈现方式(2)：交叉分析频数百分比折线图 sjPlot::plot_xtab (, type=“line”)

自变量(edu.cat5)在x轴(row), 因变量(V1.5.cat3)在y轴(col)。

如 Figure 2

library(sjPlot)#加载sjPlot套件
#使用plot_xtab函数绘图
plot_xtab(mydata$edu.cat5, mydata$V1.5.cat3, 
 #row:x轴变量(自变量)edu.cat5，col:y轴变量(因变量)V1.5.cat3
          type = "line", #呈现折线图 
          margin = "row", #呈现row%(自变量所在的%)
          weight.by = mydata$W_Raking) #加权

3.5.3 呈现方式(3)：交叉分析频数百分比雷达图 ggradar::ggradar( )

也可以下载 ggradar套件包，使用 ggradar( )函数，将事先定义好的数据档，呈现交叉分析百分比雷达图。

如 Figure 3

#输入表格资料
happy <- c("Not Happy", "Happy", "Hard to Choose") 
edu1 <- c(0.128,0.859,0.014)
edu2 <- c(0.192,0.781,0.027)
edu3 <- c(0.131,0.844,0.025)
edu4 <- c(0.121,0.851,0.028)
edu5 <- c(0.106,0.858,0.036)

#组合成data.frame
happy_edu<-data.frame(happy,edu1,edu2,edu3,edu4,edu5)

print(happy_edu)

           happy  edu1  edu2  edu3  edu4  edu5
1      Not Happy 0.128 0.192 0.131 0.121 0.106
2          Happy 0.859 0.781 0.844 0.851 0.858
3 Hard to Choose 0.014 0.027 0.025 0.028 0.036

#下载ggradar套件包
devtools::install_github("ricardo-bion/ggradar", dependencies = TRUE)

Skipping install of 'ggradar' from a github remote, the SHA1 (f99517ae) has not changed since last install.
  Use `force = TRUE` to force installation

#呈现雷达图
library(ggradar)
ggradar(happy_edu)

Ignoring unknown labels:
• size : "14"

3.6 呈现卡方检验值,对应的P值，以及Cramer’s V

为了要检验在 样本交叉表中(例如：edu.cat5教育程度; V1.5.cat3快乐感)所观察到的百分比差距，能否推论到总体，达到统计上的显著差异，我们必须进一步使用 皮尔森卡方检验 (Pearson Chi-Square Test)，来做统计检验的工作。

我们可以使用 sjPlot套件包的 tab_xtab( ) 函数，加上 show.summary=TRUE这个参数，呈现卡方检验值,对应的P值，以及Cramer’V 相关系数(0 ~ 1)等这些统计数据。

如 Table 3 (a) 与 Table 3 (b)

library(sjPlot)#加载sjPlot套件
#使用tab_xtab函数做交叉分析表
#自变量在row
tab_xtab(mydata$edu.cat5,mydata$V1.5.cat3,  
        #row自变量edu.cat5,  col因变量V1.5.cat3
        weight.by = TCS2015sc$W_Raking,#加权
        show.row.prc = TRUE,#呈现 row %
        show.legend = TRUE,#呈现cell数据说明
        show.summary = TRUE,
        # 呈现卡方检验值,P值,以及Cramer'V相关系数
        encoding = "utf8") #(for win)防止报表出现中文乱码
#自变量在col
#使用tab_xtab函数做交叉分析表
tab_xtab(mydata$V1.5.cat3, mydata$edu.cat5, 
      #row因变量V1.5.cat3, col自变量edu.cat5,  
      weight.by = mydata$W_Raking,#加权
      show.col.prc = TRUE,#呈现 col %
      show.legend = TRUE,#呈现cell数据说明
      show.summary = TRUE,  
      #呈现卡方检验值,P值,以及Cramer'V相关系数
      encoding = "utf8") #(for win)防止报表出现中文乱码

Table 3: 交叉分析卡方检验表

(a) 交叉分析卡方检验表-自变量row

受教育程度五分类	日子过的快不快乐三分类			Total
受教育程度五分类	不快乐	快乐	无法选择	Total
小学及以下	37 12.8 %	249 85.9 %	4 1.4 %	290 100 %
初中	42 19.2 %	171 78.1 %	6 2.7 %	219 100 %
高中职	69 13.1 %	444 84.4 %	13 2.5 %	526 100 %
专科	30 12.1 %	211 85.1 %	7 2.8 %	248 100 %
本科及以上	76 10.6 %	618 85.8 %	26 3.6 %	720 100 %
Total	254 12.7 %	1693 84.5 %	56 2.8 %	2003 100 %
χ²=15.289 · df=8 · Cramer's V=0.062 · p=0.054

observed values
% within 受教育程度五分类

(b) 交叉分析卡方检验表-自变量column

日子过的快不快乐三分类	受教育程度五分类					Total
日子过的快不快乐三分类	小学及以下	初中	高中职	专科	本科及以上	Total
不快乐	37 12.8 %	42 19.2 %	69 13.1 %	30 12.1 %	76 10.6 %	254 12.7 %
快乐	249 85.9 %	171 78.1 %	444 84.4 %	211 85.1 %	618 85.8 %	1693 84.5 %
无法选择	4 1.4 %	6 2.7 %	13 2.5 %	7 2.8 %	26 3.6 %	56 2.8 %
Total	290 100 %	219 100 %	526 100 %	248 100 %	720 100 %	2003 100 %
χ²=15.289 · df=8 · Cramer's V=0.062 · p=0.054

observed values
% within 受教育程度五分类

统计报表显示，Pearson ${X^2}$＝15.289, df(degree of freedom,自由度)=8, 其所对应的P = 0.054 (=0.05)，

代表我们可以拒绝 H0 (虚无假设)，不排斥或接受H1(研究假设)。也就是说，如果将样本的结果推论回总体，则民众的 快乐感 的确会因为其 受教育程度 的不同，而存在 显著差异，此一差异已超过抽样误差的范围，可以推论到总体。

当两个变量都是属于 定类变量 时(nominal scale)，一般会使用 Cramer’s V (介于0~1之间)，来衡量两变量彼此间的相关强度。

判断相关系数的经验法则是：

Note

Cramer’s V：

0 - > 零相关

0.01 ~ 0.30 - > 弱相关

0.31 ~ 0.69 - > 中相关

0.70 ~ 0.99 ->强相关

1 - >完全相关

本例中，Cramer’s V = 0.062，表示民众的 快乐感 与其 受教育程度 之间的关连，呈现一种极弱的 弱相关 关系。

综合卡方检验与 Cramer’s V 的统计分析结果，我们可以这样说：对于民众的快乐感受而言，民众的教育水平是一个达到统计显着差异的解释(预测)变量；然而，由于两变量间的相关关系很弱，所以，民众的教育水平可能不是其快乐感受的最佳解释(预测)变量。

当卡方检验结果显示 因变量 (例如：V1.5.cat3民众的”快乐感”)的确会因为我们所设置的 自变量 (例如：edu.cat5民众的”受教育程度”)而存在显著差异时，我们进一步可以观察交叉方格的 调整后标准化残差，判断到底是哪一类别或哪一细格产生显着差异。

3.7 呈现交叉方格的调整后标准化残差

要观察卡方检验的显着差异是发生在那一方格，需要借重 调整后标准化残差 (Adjusted Standardized Residual)此一參數。如果該方格内 调整后标准化残差的絕對值 大於1.96，则该方格即为显着差异之所在（方格内的百方比显着高于或低于某一类目全体的比例）。反之，即表示无显着差异。

我们可以使用 descr 套件包的 crosstab( ) 函数，来呈现交叉表中每一方格的调整后标准化残差。

#呈现交叉分析变量每个選項的频数
table(mydata$edu.cat5)


  1   2   3   4   5 
303 234 540 256 669

table(mydata$V1.5.cat3)


   1    2    3 
 253 1693   56

#使用plyr套件的mapvalues()函数赋予变量值标签
library(plyr)
mydata$edu.cat5.label <- mapvalues(mydata$edu.cat5,
         from = c("1","2", "3", "4", "5"),
         to = c("1.小学及以下","2.初中", 
                "3.高中职","4.专科", "5.本科及以上"))
mydata$V1.5.cat3.label <- mapvalues(mydata$V1.5.cat3,
         from = c("1","2","3"),
         to = c("1.不快乐", "2.快乐", "3.无法选择"))

#加载descr套件
library(descr)
#使用crosstab函数产生交叉表
crosstab(mydata$edu.cat5.label, mydata$V1.5.cat3.label, 
         #row自变量edu.cat5   ＃col因变量V1.5.cat3
                mydata$W_Raking, #加权
                prop.r = TRUE, # 显示row百分比
                chisq = TRUE, # 顯示卡方值
                asresid = TRUE) #显示调整后的标准化残差值

   Cell Contents 
|-------------------------|
|                   Count | 
|             Row Percent | 
|           Adj Std Resid | 
|-------------------------|

==========================================================
                    日子过的快不快乐三分类
受教育程度五分类    1.不快乐   2.快乐   3.无法选择   Total
----------------------------------------------------------
1.小学及以下             37      249            4     290 
                       12.8%    85.9%         1.4%   14.5%
                      0.043    0.682       -1.582         
----------------------------------------------------------
2.初中                   42      171            6     219 
                       19.2%    78.1%         2.7%   10.9%
                      3.062   -2.792       -0.053         
----------------------------------------------------------
3.高中职                 69      444           13     526 
                       13.1%    84.4%         2.5%   26.3%
                      0.351   -0.083       -0.525         
----------------------------------------------------------
4.专科                   30      211            7     248 
                       12.1%    85.1%         2.8%   12.4%
                     -0.295    0.259        0.027         
----------------------------------------------------------
5.本科及以上             76      618           26     720 
                       10.6%    85.8%         3.6%   35.9%
                     -2.141    1.214        1.658         
----------------------------------------------------------
Total                   254     1693           56    2003 
==========================================================

Statistics for All Table Factors

Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 = 15.28888      d.f. = 8      p = 0.0538 

        Minimum expected frequency: 6.122816

从上表中，我们可以看到：

受教度程度为初中者，所属方格的 Adj Std Resid 值为 +3.062（绝对值大于1.96）, 代表初中教育程度者自觉 不快乐 的比例(19.2%），显着高于 (▲)全体自觉 不快乐 的比例。

另一方面，本科及以上 者所属方格的 Adj Std Resid 值为 -2.141（绝对值大于1.96），代表本科及以上教育程度者自觉 不快乐 的比例 (10.6%），显着低于 (▼) 全体自觉 不快乐 的比例。

同理可以得知，初中教育程度者自觉快乐的比例(78.1%），显着低于 (▼)全体自觉快乐的比例。

所以，通过皮尔森卡方检验以及调整后的标准化残差分析，我们可以得到以下的结论：

民众的 快乐感 会因为其 受教育程度 的不同，而存在统计上的显着差异（Pearson $X^2$ = 15.289, df = 8, P = 0.005）。
具体来说，初中教育程度者有较高的 不快乐 比例，较低的快乐比例。 本科及以上 教育程度者有较低的 不快乐 比例。

4 2 * 2 交叉分析

[2 * 2] 表格交叉分析，基本上和 [R * C 表格] 有着相同的分析步骤。

但[2 * 2]交叉分析在卡方值的选择上，主要是以 Yates’s Continuity correction 的卡方值，来代替 Pearson卡方值；而在相关系数上，则是选择使用 φ (读做Phi)(0 ~ 无限大)，来代替 Cramer’s V值，做为相关系数的测量统计值。此处以具体例子来展示[2 * 2]表格的交叉分析步骤。

4.1 进行编码转换

4.1.1 将因变量重新编码：V1.5(过的快不快乐)(5分类) - > V1.5.cat2(2分类)

我们要做的是将 V1.5 变量的原始五分类(1.非常不快乐, 2. 不快乐, 3.快乐, 4. 非常快乐, 94.无法选择)，转换为一个新的二分类变量 V1.5.cat2(1.快乐, 0.其他)

此处可以使用 sjmisc 套件包，将原始变量(V1.5)依照我们的编码规则(3,4=1; else=0)，重新编码(rec)为另一新的变量(V1.5.cat2)；给予变量名称(var.label)和变量标签名称(val.labels)；将新的变量(V1.5.ca2)设置为类别变量(as.num=FALSE)；检查编码转换是否正确(flat_table)；并呈现重新编码后的频数分布(frq)。

#加载sjmisc套件
library(sjmisc) 
mydata$V1.5.cat2 <- rec(mydata$V1.5,  #使用rec函数重新编码
               rec = "else=0; 3,4 = 1", #设置编码规则
               var.label = "日子过的快不快乐二分类",
               val.labels = c("其他","快乐"),
               as.num = FALSE) #设定为类别变量

flat_table(mydata, V1.5, V1.5.cat2)#检查编码转换是否正确

           V1.5.cat2 其他 快乐
V1.5                          
非常不快樂             21    0
不快乐                232    0
快乐                    0 1498
非常快乐                0  195
无法选择               56    0

frq(mydata$V1.5.cat2, 
    weights = mydata$W_Raking) #频数分布(加权)

日子过的快不快乐二分类 (xw) <categorical> 
# total N=2002 valid N=2002 mean=0.85 sd=0.36

Value | Label |    N | Raw % | Valid % | Cum. %
-----------------------------------------------
    0 |  其他 |  310 | 15.48 |   15.48 |  15.48
    1 |  快乐 | 1692 | 84.52 |   84.52 | 100.00
 <NA> |  <NA> |    0 |  0.00 |    <NA> |   <NA>

4.1.2 将自变量重新编码：A1(性别) -> sex(2分类)

仿照前述程序将原始变量(A1)依照我们的编码规则(1=1; 2=2)，重新编码(rec)为另一新的变量(sex)，并赋予变量及变量值标签。

#加载sjmisc套件
library(sjmisc) 
frq(mydata$A1)

性别 (x) <categorical> 
# total N=2002 valid N=2002 mean=1.54 sd=0.50

Value | Label |    N | Raw % | Valid % | Cum. %
-----------------------------------------------
    1 |    男 |  928 | 46.35 |   46.35 |  46.35
    2 |    女 | 1074 | 53.65 |   53.65 | 100.00
 <NA> |  <NA> |    0 |  0.00 |    <NA> |   <NA>

mydata$sex <- rec(mydata$A1,  #使用rec函数重新编码
               rec = "1=1; 2 = 2",  #设置编码规则
               var.label = "sex性别",
               val.labels = c("男","女"),
               as.num = FALSE) #设定为类别变量

flat_table(mydata, A1,sex)#检查编码转换是否正确

   sex   男   女
A1              
男      928    0
女        0 1074

frq(mydata$sex, 
    weights = mydata$W_Raking)#频数分布(加权)

sex性别 (xw) <categorical> 
# total N=2002 valid N=2002 mean=1.51 sd=0.50

Value | Label |    N | Raw % | Valid % | Cum. %
-----------------------------------------------
    1 |    男 |  988 | 49.35 |   49.35 |  49.35
    2 |    女 | 1014 | 50.65 |   50.65 | 100.00
 <NA> |  <NA> |    0 |  0.00 |    <NA> |   <NA>

我们可以使用两种方式来呈现[2 * 2]交叉表格： (1)使用 sjPlot 套件包的 tab_xtab( )函数呈现， (2)使用 descr 套件包的 crosstab( )函数呈现；藉以比较两种统计表格呈现方式的不同之处。

#调用sjPlot套件的tab_xtab()函数，做交叉分析(1)
#加载sjPlot套件
library(sjPlot)
#使用tab_xtab函数做交叉表
tab_xtab(mydata$sex, mydata$V1.5.cat2, 
        #row自变量sex,  col因变量V1.5.cat2
        weight.by = mydata$W_Raking,#加权
        show.row.prc = TRUE, #呈现 row % 
        show.legend = TRUE, #呈现 cell数据说明 
        show.summary = TRUE, #呈现卡方检验值
        encoding = "utf8") #(for win)防止报表出现中文乱码
#加载plyr套件
library(plyr)
#使用plyr套件的mapvalues函数赋予变量值标签
mydata$sex.label <- mapvalues(mydata$sex,
         from = c("1","2"),
         to = c("1.男", "2.女"))
mydata$V1.5.cat2.label <- mapvalues(mydata$V1.5.cat2,
         from = c("0","1"),
         to = c("0.其他", "1.快乐"))

#调用descr套件包的crosstab()函数，做交叉分析(2)
library(descr) #加载descr套件
#使用crosstab函数做交叉表
crosstab(mydata$sex.label, mydata$V1.5.cat2.label, 
                mydata$W_Raking, #加权
                prop.r = TRUE, #显示row %
                chisq = TRUE)# 显示卡方值

   Cell Contents 
|-------------------------|
|                   Count | 
|             Row Percent | 
|-------------------------|

==================================
           日子过的快不快乐二分类
sex性别    0.其他   1.快乐   Total
----------------------------------
1.男         170      818     988 
            17.2%    82.8%   49.4%
----------------------------------
2.女         139      874    1013 
            13.7%    86.3%   50.6%
----------------------------------
Total        309     1692    2001 
==================================

Statistics for All Table Factors

Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 = 4.651843      d.f. = 1      p = 0.031 

Pearson's Chi-squared test with Yates' continuity correction 
------------------------------------------------------------
Chi^2 = 4.388788      d.f. = 1      p = 0.0362 
        Minimum expected frequency: 152.5697

sex性别	日子过的快不快乐二分类	Total
其他	快乐
男	170 17.2 %	818 82.8 %	988 100 %
女	139 13.7 %	874 86.3 %	1013 100 %
Total	309 15.4 %	1692 84.6 %	2001 100 %
χ²=4.389 · df=1 · &phi=0.048 · p=0.036

使用 descr::crosstab( )所产生的 Pearson’s Chi-squared test 的值为 4.651843，Yates’s continuity correction Chi-squared test的值为 4.38,

此时，我们应以 Yates’s continuity correction Chi-squared test的值(4.38)为准，其所对应的P-value = 0.0362 < 0.05，代表民众日子过的快不快乐，会因为其性别不同而有差异。

另外，在 Table 4 使用 sjPlot::tab_xtab( ) 所产生的 $X^2$＝4.389，此值即是 Yates’s continuity correction 的卡方值。另外，亦呈现了φ(phi) 相关系数的大小(0.048)。

学习者可继续按照前述[R * C] 表格中所介绍的方式，呈现[2 * 2]交叉方格的调整后标准化残差，具体观察显着差异之所在。

特别要说的是，在 [2 * 2]交叉表格中，如果有任何一方格的期望值 < 5，卡方值的选用就需要使用 Fisher’s exact test (费雪精确检验)来代替。统计软件都会自动跑出 Fisher’s exact test，供研究者选用。

另一方面，[2 * 2] 交叉表格还可以用 发生比 (odds) 与 发生比的倍数变化 (odds ratio) 的概念来做分析。

发生比(odds)，其意为：“发生某一事件的概率(p)，相对于不发生该事件概率(1-p)的比例”，即 $\frac{p}{1 - p}$；

发生比的倍数变化 (odds ratio)，其意为 “发生比的比值”。

以 Table 4 为例，

我们可以得知，

男性自觉 快乐相对于其他 的比例＝4.81 (82.8% / 17.2%)，此 4.81 称为 发生比 (odds)；
女性自觉 快乐相对于其他 的比例＝6.30 (86.3% / 13.7%)，此 6.30 称为 发生比 (odds)；
全体自觉 快乐相对于其他 的比例＝5.49(84.6% / 15.4%)，此 5.49 称为 发生比 (odds)。

所以，

男性自觉 快乐相对于其他 的比例，是女性的 0.763倍 (4.81 / 6.30)；
男性自觉 快乐相对于其他 的比例，是全体的 0.876倍 (4.81 / 5.49)；

此 0.763, 0.876 称为 发生比的倍数变化 (odds ratio)。

女性自觉 快乐相对于其他 的比例，是男性的 1.31倍 (6.30 /4.81)；
女性自觉 快乐相对于其他 的比例，是全体的 1.15倍 (6.30 / 5.49)；

此 1.31, 1.15 亦称为 发生比的倍数变化(odds ratio)。

简单来说，就是 女性比男性更觉得自己过的比较快乐。

5 交叉分析卡方检验的前提及违反前提时之处理

进行交叉分析卡方检验时，有一个重要的前提是： 期望值小于5的方格数，不能超过全部方格数的20%；

如果超过，研究者需要将变量进行适度的类别合并，再进行卡方检验，这样得到的卡方值($X^2$)以及 P 值检验，才会是比较正确的。

什么是 期望值 呢？ 期望值 (expected value) 指的是如果兩变量无关，则交叉表格中每一方格内理论上应该要出现的期望次数。

具体来说，就是如果自变量在该方格内的反应比例与总体的反应比例并无二致的话，理论上该方格应该要出现的期望次数。

至于合并类别的原则，一般的经验是，自变量与因变量的交叉方格数目，最好不要超过 40(R * C <= 40)。实务上，自变量与因变量的类别数，各以6类以下为宜。

5.1 一个不佳的交叉分析卡方检验例子

此处以 age_strata (年龄分层,7分类)与 V1.5 (民众过的快不快乐, 5分类) 的[7 * 5] 交叉分析，做为一个不佳的卡方检验例子。

#使用plyr套件的mapvalues函数赋予变量值标签
#加载plyr套件
library(plyr)
mydata$age_strata.label <- mapvalues(mydata$age_strata,
         from = c("1","2", "3", "4", "5", "6", "7"),
         to = c("18-19", "20-29", "30-39", "40-49",
                "50-59","60-69", "70以上"))
mydata$V1.5.label <- mapvalues(mydata$V1.5,
         from = c("1","2", "3", "4", "94"),
         to = c("非常不快乐", "不快乐", "快乐", 
                "非常快乐", "无法选择"))

#交叉分析
#加载descr套件
library(descr)
#使用crosstab函数做交叉表
crosstab(mydata$age_strata.label, mydata$V1.5.label,
                mydata$W_Raking, #加权
                expected = TRUE, #呈現期望值
                prop.r = TRUE, #呈現row %
                chisq = TRUE)# 呈現卡方值

Warning in chisq.test(tab, correct = FALSE, ...): Chi-squared近似算法有可能不准

   Cell Contents 
|-------------------------|
|                   Count | 
|         Expected Values | 
|             Row Percent | 
|-------------------------|

=====================================================================
            V1.5. 整体来说，你觉得目前的日子过得快不快乐？
年龄分层    非常不快乐   不快乐    快乐   非常快乐   无法选择   Total
---------------------------------------------------------------------
18-19               0        5      41         18          2      66 
                  0.7      7.6    49.1        6.7        1.8         
                  0.0%     7.6%   62.1%      27.3%       3.0%    3.3%
---------------------------------------------------------------------
20-29               7       30     239         43         10     329 
                  3.6     38.1   244.8       33.3        9.2         
                  2.1%     9.1%   72.6%      13.1%       3.0%   16.4%
---------------------------------------------------------------------
30-39               5       56     298         35         10     404 
                  4.4     46.8   300.6       40.9       11.3         
                  1.2%    13.9%   73.8%       8.7%       2.5%   20.2%
---------------------------------------------------------------------
40-49               3       58     277         21         14     373 
                  4.1     43.2   277.5       37.8       10.4         
                  0.8%    15.5%   74.3%       5.6%       3.8%   18.6%
---------------------------------------------------------------------
50-59               5       42     275         38         12     372 
                  4.1     43.1   276.8       37.7       10.4         
                  1.3%    11.3%   73.9%      10.2%       3.2%   18.6%
---------------------------------------------------------------------
60-69               2       26     193         31          6     258 
                  2.8     29.9   192.0       26.1        7.2         
                  0.8%    10.1%   74.8%      12.0%       2.3%   12.9%
---------------------------------------------------------------------
70以上              0       15     168         17          2     202 
                  2.2     23.4   150.3       20.5        5.6         
                  0.0%     7.4%   83.2%       8.4%       1.0%   10.1%
---------------------------------------------------------------------
Total              22      232    1491        203         56    2004 
=====================================================================

Statistics for All Table Factors

Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 = 59.6325      d.f. = 24      p = 7.2e-05 

        Minimum expected frequency: 0.7245509 
Cells with Expected Frequency < 5: 8 of 35 (22.85714%)

上表显示，交叉表中，期望值小于 5 的方格数约有8个，占全部35方格数的22.85%，已违反卡方检验的分析前提，据此而计算出来的卡方值，有可能不准确。因此，研究者需要合并类别，重新进行交叉分析与卡方检验，方为正道。

5.2 合并类别后重新进行交叉分析卡方检验

5.2.1 进行编码转换

5.2.1.1 将因变量重新编码:V1.5(过的快不快乐)(5分类) - > V1.5.cat3(3分类)

#加载sjmisc套件
library(sjmisc)
mydata$V1.5.cat3 <- rec(mydata$V1.5,#使用rec函数重新编码
          rec = "1,2=1; 3,4=2; 94=3;", #设置编码规则
          var.label = "V1.5.cat3快不快乐三分类",
          val.labels = c("不快乐","快乐", "无法选择"),
          as.num = FALSE) #设定为类别变量

5.2.1.2 将自变量重新编码:age_strata(年龄分层)(7分类) - > age.cat5(5分类)

#加载sjmisc套件
library(sjmisc)
mydata$age.cat5 <- rec(mydata$age_strata, #rec函数重新编码
                rec = "1,2=1; 3=2; 4=3; 5=4; 6,7=5", 
                var.label = "age.cat5年龄五分类",
                val.labels = c("18-29","30-39", "40-49",
                               "50-59", "60以上"),
                as.num = FALSE) #设定为类别变量

5.2.2 呈现交叉分析及卡方检验(合并类别后)

经过对于变量类别进行合并后，统计报表显示，重新进行的卡方检验结果，并无任何一个方格的期望值小于或等于5。所以这是一个合格的卡方检验分析。

如 Table 5

library(plyr)#加载plyr套件
#使用plyr套件的mapvalues函数赋予变量值标签
mydata$age.cat5.label <- mapvalues(mydata$age.cat5,
         from = c("1","2", "3", "4","5"), #界定变量值
         to = c("18-29", "30-39", "40-49", "50-59",
                "60及以上"))#給予变量标签

mydata$V1.5.cat3.label <- mapvalues(mydata$V1.5.cat3,
      from = c("1","2", "3"),#界定变量值
      to = c("不快樂", "快乐", "無法选择"))#給予变量标签

#交叉分析及卡方检验
library(descr)#加载descr套件
crosstab(mydata$age.cat5.label, mydata$V1.5.cat3.label, 
         #自变量age.cat5.label  #因变量V1.5.cat3.label
                mydata$W_Raking, #加权
                expected = TRUE, #显示期望值
                prop.r = TRUE, #显示row %
                chisq = TRUE)# 显示卡方值

   Cell Contents 
|-------------------------|
|                   Count | 
|         Expected Values | 
|             Row Percent | 
|-------------------------|

=======================================================
                      V1.5.cat3快不快乐三分类
age.cat5年龄五分类    不快樂    快乐   無法选择   Total
-------------------------------------------------------
18-29                    41     341         12     394 
                       49.8   333.2       11.0         
                       10.4%   86.5%       3.0%   19.7%
-------------------------------------------------------
30-39                    61     333         10     404 
                       51.1   341.6       11.3         
                       15.1%   82.4%       2.5%   20.2%
-------------------------------------------------------
40-49                    61     298         14     373 
                       47.1   315.4       10.4         
                       16.4%   79.9%       3.8%   18.6%
-------------------------------------------------------
50-59                    47     313         12     372 
                       47.0   314.6       10.4         
                       12.6%   84.1%       3.2%   18.6%
-------------------------------------------------------
60及以上                 43     408          8     459 
                       58.0   388.2       12.8         
                        9.4%   88.9%       1.7%   22.9%
-------------------------------------------------------
Total                   253    1693         56    2002 
=======================================================

Statistics for All Table Factors

Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 = 17.35951      d.f. = 8      p = 0.0266 

        Minimum expected frequency: 10.40559

#加载sjPlot套件
library(sjPlot)
#使用tab_xtab函数做交叉表
tab_xtab(mydata$age.cat5.label, mydata$V1.5.cat3.label, 
        weight.by = mydata$W_Raking,#加权
        show.row.prc = TRUE, # 呈现 row % 
        show.exp = TRUE,#呈现各cell期望值
        show.legend = TRUE,#呈现cell数据说明
        show.summary = TRUE, #呈现卡方检验值,P值,相关系数
        encoding = "utf8") #(for win)防止报表出现中文乱码

age.cat5年龄五分类	V1.5.cat3快不快乐三分类	Total
不快樂	快乐	無法选择
18-29	41 50 10.4 %	341 333 86.5 %	12 11 3 %	394 394 100 %
30-39	61 51 15.1 %	333 342 82.4 %	10 11 2.5 %	404 404 100 %
40-49	61 47 16.4 %	298 315 79.9 %	14 10 3.8 %	373 373 100 %
50-59	47 47 12.6 %	313 315 84.1 %	12 10 3.2 %	372 372 100 %
60及以上	43 58 9.4 %	408 388 88.9 %	8 13 1.7 %	459 459 100 %
Total	253 253 12.6 %	1693 1693 84.6 %	56 56 2.8 %	2002 2002 100 %
χ²=17.360 · df=8 · Cramer's V=0.066 · p=0.027

rm(list = ls()) #清除记忆体中所有物件与资料

本章小结

本章所使用到的R套件包与函数摘录如下表。

套件包	函数	说明
内建	setwd ( )	设置当前工作目录
	getwd ( )	显示当前工作目录
	install.packages ( )	安装套件包
	library ( )	加载套件包
	save ( )	将档案存成R数据资料档格式(*.rda)
	load ( )	载入R数据资料档(*.rda)
	dim ( )	查看R数据中包含几笔资料，几个变量
	names ( )	查看变量名称
	table ( )	呈现频数
	levels ( )	呈现类别变量的水平(有几类)
	data.frame ( )	将数据(例如：向量)转换为数据框
	rm (list = ls( ))	去除R当前工作环境中所有的数据与物件
dplyr	select ( )	选择分析变量
plyr	mapvalues ( )	将变量值进行配对编码
sjlabelled	read_spss ( )	将SPSS格式数据资料档(*.sav)汇入到R
sjmisc	frq ( )	呈现频数分布表
	rec ( )	将变量值进行转换编码
sjPlot	tab_xtab (, show.summary=F)	呈现交叉分析表(不含卡方检验值)
	tab_xtab (, show.summary=T)	呈现交叉分析表(含卡方检验值)
	plot_xtab (, type=“bar”)	呈现交叉分析频数百分比条形图
	plot_xtab (, type=“line”)	呈现交叉分析频数百分比折线图
ggradar	ggradar ( )	呈现交叉分析百分比雷达图
descr	crosstab ( )	呈现交叉分析表