2 列联分析
2.1 主要理论
对于二维列联表：
拟合优度检验：样本是从总体的不同类别中分别抽取，研究不同类别的目标量之间是否存在显著性差异。
独立性检验：两个变量间是否存在联系。
两者都是利用Pearson \(\chi^2\)统计量，在R中使用函数chisq.test。
chisq.test(x, y = NULL, correct = TRUE,
p = rep(1/length(x), length(x)), rescale.p = FALSE,
simulate.p.value = FALSE, B = 2000)用法参考：https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/chisq.test 。主要用法为将列联表矩阵赋予参数x。
使用注意：
要求样本量必须足够大，特别是每个单元中的期望频数不能过小。一般需要满足两个准则：
如果只有两个单元，每个单元的期望频数必须大于等于5.在2x2列联表中，有一个期望值小于5且样本量小于40，或几个期望值小于5，或有一个期望值小于1，均不宜作\(\chi^2\)检验。
有两个以上的单元，20%的单元期望频数小于5，不能应用\(\chi^2\)检验。
对于三维列联表：
- 独立性检验：建立对数线性模型。
利用MASS包的loglin()函数，用法详见：https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/loglin 。主要用法为：
# a为一个三维数组
loglin(a,list(1,2,3)) # 检验X Y Z独立性
loglin(a,list(1,c(2,3))) # 检验X与（Y Z）独立性
loglin(a,list(c(1,2),c(1,3))) # 检验(X Y)(X Z)独立性，即给定X检验Y Z独立性并且可以计算Pearson统计量以及LRT统计量，当值较大时拒绝原假设（独立）。
2.2 习题4
4.1
从两台机器生产的产品中随机抽样，来检查两台机器生产的产品次品率是否有差异，列联表如下：
| 次品 | 正品 | 合计 | |
|---|---|---|---|
| 机器1 | 12 | 88 | 100 |
| 机器2 | 16 | 74 | 90 |
| 合计 | 28 | 162 | 190 |
试在显著性水平\(\alpha=0.05\)下，检验两台机器生产的产品次品率是否有差异？
解：
二维列联表的拟合优度检验。
x <- c(12,88,16,74)
x <- matrix(x,2,2,byrow = TRUE)
chisq.test(x)
#>
#> Pearson's Chi-squared test with Yates' continuity correction
#>
#> data: x
#> X-squared = 0.8, df = 1, p-value = 0.4\(p\)值大于0.05，应当接受原假设，即认为两台机器生产的产品次品率没有差异。
4.2
| 赞成该方案 | 反对该方案 | 合计 | |
|---|---|---|---|
| 一分公司 | 68 | 32 | 100 |
| 二分公司 | 75 | 45 | 120 |
| 三分公司 | 57 | 33 | 90 |
| 四分公司 | 79 | 31 | 100 |
| 合计 | 279 | 141 | 420 |
试在显著性水平\(\alpha=0.1\)下，检验四个分公司对改革方案的赞成比例是否一致。
解：
二维列联表的拟合优度检验。
x <- c(68,75,57,79,32,45,33,31)
x <- matrix(x,4,2)
chisq.test(x)
#>
#> Pearson's Chi-squared test
#>
#> data: x
#> X-squared = 3, df = 3, p-value = 0.4\(p\)值大于0.1，应当接受原假设，即认为四个分公司对改革方案的赞成比例一致。
4.3
| 小于等于30（岁） | 31-50（岁） | >50（岁） | 合计 | |
|---|---|---|---|---|
| 体育类 | 83 | 91 | 41 | 215 |
| 电视剧类 | 70 | 86 | 38 | 194 |
| 综艺类 | 45 | 15 | 10 | 70 |
| 合计 | 198 | 192 | 89 | 479 |
试在显著性水平\(\alpha=0.05\)下，检验不同观众对三类节目的关注率是否一样。
解：
二维列联表的拟合优度检验。
x <- c(83,70,45,91,86,15,41,38,10)
x <- matrix(x,3,3)
chisq.test(x)
#>
#> Pearson's Chi-squared test
#>
#> data: x
#> X-squared = 19, df = 4, p-value = 9e-04\(p\)值小于0.05，应当拒绝原假设，即认为不同观众对三类节目的关注率不一样。
4.4
| 得肺癌 | 未得肺癌 | 合计 | |
|---|---|---|---|
| 抽烟 | 8 | 19 | 27 |
| 未抽烟 | 1 | 17 | 18 |
| 合计 | 9 | 36 | 45 |
试在显著性水平\(\alpha=0.05\)下，检验抽烟对得肺癌是否有显著影响。
解：
利用\(\chi^2\)检验，一时间分不清是拟合优度检验还是独立性检验。。
先进行\(\chi^2\)分布的期望值准则验证。
# 输入二维列联表矩阵
# 返回值：
# n 总和
# expectation 期望频数矩阵
# rule 是否满足准则
chisq.expectation <- function(x){
res <- list()
r <- nrow(x)
c <- ncol(x)
rs <- apply(x,1,sum) # 每行求和
cs <- apply(x,2,sum) # 每列求和
as <- sum(x) # 总和
res$n <- as
res$expectation <- rs%*%t(cs)/as # 期望频数矩阵
if(r==2&c==1|r==1&c==2){
res$rule <- !sum(res$expectation<5)
}else if(r==2&c==2){
if(sum(res$expectation<5)==1&as<40){res$rule <-FALSE}
else if(sum(res$expectation<5)>1|sum(res$expectation<1)>0) {res$rule <-FALSE}
else {res$rule <-TRUE}
}else{
res$rule <- (mean(res$expectation<5)<0.2)
}
res
}chisq.expectation(x)
#> $n
#> [1] 479
#>
#> $expectation
#> [,1] [,2] [,3]
#> [1,] 88.9 86.2 39.9
#> [2,] 80.2 77.8 36.0
#> [3,] 28.9 28.1 13.0
#>
#> $rule
#> [1] TRUE满足准则。
x <- c(8,1,19,17)
x <- matrix(x,2,2)
chisq.test(x)
#> Warning in chisq.test(x): Chi-squared approximation may be incorrect
#>
#> Pearson's Chi-squared test with Yates' continuity correction
#>
#> data: x
#> X-squared = 3, df = 1, p-value = 0.1\(p\)值大于0.05，应当接受原假设，即认为抽烟对得肺癌没有显著影响。
4.5
| 肝炎 | 肝硬化 | 对照（非肝病患者） | 合计 | |
|---|---|---|---|---|
| O型 | 98 | 38 | 289 | 425 |
| A型 | 67 | 41 | 262 | 370 |
| B型 | 13 | 8 | 57 | 78 |
| AB型 | 18 | 12 | 30 | 60 |
| 合计 | 196 | 99 | 638 | 933 |
试在显著性水平\(\alpha=0.05\)下，检验血型与肝病之间是否存在关联。
解：
先进行\(\chi^2\)分布的期望值准则检验。
x <- c(98,67,13,18,38,41,8,12,289,262,57,30)
x <- matrix(x,4,3)
x
#> [,1] [,2] [,3]
#> [1,] 98 38 289
#> [2,] 67 41 262
#> [3,] 13 8 57
#> [4,] 18 12 30
chisq.expectation(x)
#> $n
#> [1] 933
#>
#> $expectation
#> [,1] [,2] [,3]
#> [1,] 89.3 45.10 290.6
#> [2,] 77.7 39.26 253.0
#> [3,] 16.4 8.28 53.3
#> [4,] 12.6 6.37 41.0
#>
#> $rule
#> [1] TRUE满足准则，进行独立性检验。
chisq.test(x)
#>
#> Pearson's Chi-squared test
#>
#> data: x
#> X-squared = 15, df = 6, p-value = 0.02\(p\)值小于0.05，应当拒绝原假设，即认为血型与肝病之间存在关联。
4.6
A、B、C是候选学生，列联表展示了各学院投票情况：
| A | B | C | 合计 | |
|---|---|---|---|---|
| 一学院 | 24 | 23 | 12 | 59 |
| 二学院 | 24 | 14 | 10 | 48 |
| 三学院 | 17 | 8 | 13 | 38 |
| 四学院 | 27 | 19 | 9 | 55 |
| 合计 | 92 | 64 | 44 | 200 |
问在显著性水平\(\alpha=0.01\)下，我们能否认为投票结果与投票者在什么学院相互独立？
解：
先进行\(\chi^2\)分布的期望值准则检验。
x <- c(24,24,17,27,23,14,8,19,12,10,13,9)
x <- matrix(x,4,3)
chisq.expectation(x)
#> $n
#> [1] 200
#>
#> $expectation
#> [,1] [,2] [,3]
#> [1,] 27.1 18.9 12.98
#> [2,] 22.1 15.4 10.56
#> [3,] 17.5 12.2 8.36
#> [4,] 25.3 17.6 12.10
#>
#> $rule
#> [1] TRUE满足准则，下面进行独立性检验：
chisq.test(x)
#>
#> Pearson's Chi-squared test
#>
#> data: x
#> X-squared = 7, df = 6, p-value = 0.4\(p\)值大于0.05，应当接受原假设，即认为投票结果与投票者在什么学院相互独立。
4.7
2x2x2列联表，三个变量分别为护士衣着颜色（X）、儿童性别（Y）、护士性别（Z）。
护士性别 女：
| 护士衣着颜色/儿童性别 | 女 | 男 |
|---|---|---|
| 花衣 | 36 | 8 |
| 白衣 | 12 | 4 |
护士性别 男：
| 护士衣着颜色/儿童性别 | 女 | 男 |
|---|---|---|
| 花衣 | 25 | 13 |
| 白衣 | 1 | 0 |
试在显著性水平\(\alpha=0.01\)下，检验这三个变量之间哪些是独立的，哪些是不独立的。
解：
三维列联表的独立性检验：
x <- c(36,12,8,4,25,1,13,0)
x <- array(x,c(2,2,2)) # 转为三维数组
loglin(x,list(1,2,3))
#> 2 iterations: deviation 1.42e-14
#> $lrt
#> [1] 15.3
#>
#> $pearson
#> [1] 12.6
#>
#> $df
#> [1] 4
#>
#> $margin
#> $margin[[1]]
#> [1] 1
#>
#> $margin[[2]]
#> [1] 2
#>
#> $margin[[3]]
#> [1] 3library(tidyverse)
#> -- Attaching packages ------------------
#> √ ggplot2 3.3.2 √ purrr 0.3.4
#> √ tibble 3.0.3 √ dplyr 1.0.2
#> √ tidyr 1.1.1 √ stringr 1.4.0
#> √ readr 1.3.1 √ forcats 0.5.0
#> -- Conflicts -- tidyverse_conflicts() --
#> x dplyr::filter() masks stats::filter()
#> x dplyr::lag() masks stats::lag()
iterfun <- function(l){
x <- c(36,12,8,4,25,1,13,0)
x <- array(x,c(2,2,2))
loglin(x,margin = l,print = FALSE)[1:3]
}
models <- list(list(1,2,3),list(1,c(2,3)),list(2,c(1,3)),list(3,c(1,2)),
list(c(1,2),c(1,3)),list(c(2,1),c(2,3)),list(c(1,3),c(2,3)))
models_name <- c("(X,Y,Z)","(X,YZ)","(Y,XZ)","(Z,XY)","(XY,XZ)","(XY,YZ)","(XZ,YZ)")
df <- map(models,iterfun) %>% # 对所有模型应用loglin函数
unlist() %>%
matrix(ncol = 3,byrow = TRUE) %>%
as_tibble() %>% # 转为tibble
rename(lrt=V1,pearson=V2,df=V3) %>% # 重命名
mutate(models=models_name) %>% # 加入对应模型名称
select(models,everything()) %>% # 调整列顺序
mutate(cv = qchisq(1-0.01,df)) %>% # 计算显著性水平为0.01的临界值
mutate(lrt_test=(lrt>cv),pearson_test=(pearson>cv)) # 检验结果，是否接受原假设
#> Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
#> Using compatibility `.name_repair`.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_warnings()` to see where this warning was generated.
knitr::kable(df,digits = 3)| models | lrt | pearson | df | cv | lrt_test | pearson_test |
|---|---|---|---|---|---|---|
| (X,Y,Z) | 15.26 | 12.638 | 4 | 13.28 | TRUE | FALSE |
| (X,YZ) | 13.06 | 10.213 | 3 | 11.35 | TRUE | FALSE |
| (Y,XZ) | 3.35 | 3.119 | 3 | 11.35 | FALSE | FALSE |
| (Z,XY) | 15.22 | 12.594 | 3 | 11.35 | TRUE | TRUE |
| (XY,XZ) | 3.31 | 3.077 | 2 | 9.21 | FALSE | FALSE |
| (XY,YZ) | 13.03 | 10.370 | 2 | 9.21 | TRUE | TRUE |
| (XZ,YZ) | 1.15 | 0.854 | 2 | 9.21 | FALSE | FALSE |
其中lrt_test,pearson_test表示lrt统计量跟pearson统计量的检验结果是否接受原假设（变量独立），可以看到两者有时并非一致的。
4.8
2x2x3列联表，三个变量分别为近视与否（X）、学生性别（Y）、学校（Z）。
学校 甲：
| 近视与否/学生性别 | 男 | 女 |
|---|---|---|
| 近视 | 55 | 58 |
| 不近视 | 45 | 41 |
学校 乙：
| 近视与否/学生性别 | 男 | 女 |
|---|---|---|
| 近视 | 66 | 85 |
| 不近视 | 87 | 70 |
学校 丙：
| 近视与否/学生性别 | 男 | 女 |
|---|---|---|
| 近视 | 66 | 50 |
| 不近视 | 41 | 39 |
试在显著性水平\(\alpha=0.01\)下检验哪些变量独立，哪些不独立。
解：
仿照上一题做法。
library(tidyverse)
iterfun <- function(l){
x <- c(55,45,58,41,66,87,85,70,66,41,50,39)
x <- array(x,c(2,2,3))
loglin(x,margin = l,print = FALSE)[1:3]
}
models <- list(list(1,2,3),list(1,c(2,3)),list(2,c(1,3)),list(3,c(1,2)),
list(c(1,2),c(1,3)),list(c(2,1),c(2,3)),list(c(1,3),c(2,3)))
models_name <- c("(X,Y,Z)","(X,YZ)","(Y,XZ)","(Z,XY)","(XY,XZ)","(XY,YZ)","(XZ,YZ)")
df <- map(models,iterfun) %>% # 对所有模型应用loglin函数
unlist() %>%
matrix(ncol = 3,byrow = TRUE) %>%
as_tibble() %>% # 转为tibble
rename(lrt=V1,pearson=V2,df=V3) %>% # 重命名
mutate(models=models_name) %>% # 加入对应模型名称
select(models,everything()) %>% # 调整列顺序
mutate(cv = qchisq(1-0.01,df)) %>% # 计算显著性水平为0.01的临界值
mutate(lrt_test=(lrt>cv),pearson_test=(pearson>cv)) # 检验结果，是否接受原假设
knitr::kable(df,digits = 3)| models | lrt | pearson | df | cv | lrt_test | pearson_test |
|---|---|---|---|---|---|---|
| (X,Y,Z) | 12.18 | 12.12 | 7 | 18.5 | FALSE | FALSE |
| (X,YZ) | 10.91 | 10.90 | 5 | 15.1 | FALSE | FALSE |
| (Y,XZ) | 6.36 | 6.35 | 5 | 15.1 | FALSE | FALSE |
| (Z,XY) | 10.85 | 10.93 | 6 | 16.8 | FALSE | FALSE |
| (XY,XZ) | 5.04 | 5.03 | 4 | 13.3 | FALSE | FALSE |
| (XY,YZ) | 9.59 | 9.54 | 4 | 13.3 | FALSE | FALSE |
| (XZ,YZ) | 5.10 | 5.09 | 3 | 11.3 | FALSE | FALSE |
好像…所有变量组合都不独立…