categorical-data-analysis.qmd

# 分类数据的分析 {#sec-categorical-data-analysis}

```{r}
#| echo: false

source("_common.R")
```

::: hidden
$$
 \def\bm#1{{\boldsymbol #1}}
$$
:::

1.  最常见的两个广义线性模型：泊松和逻辑回归
2.  理论公式、R 输出及其解释，应用案例
3.  与计数/离散数据的假设检验的关系
4.  辛普森悖论，分类数据处理，高维列联表的压缩和分层，边际和条件
5.  泰坦尼克号 4x2x2x2 高维复杂列联表分析

```{r}
library(MASS)
```

计数数据，通俗来说，对象是一个一个或一份一份的，可数的、离散的数据，比如人数。列联表来组织数据，分二维和多维的情况。

```{r}
#| eval: false
#| echo: false

# 逻辑回归模型 二分 dichotomous logit regression
glm(family = binomial(link = "logit"))
# Probit 回归模型 二分 dichotomous probit regression
glm(family = binomial(link = "probit"))
# 泊松回归
glm(family = poisson(link = "log"))
# 多项逻辑回归模型
nnet::multinom()
# 比例比/比值比/优势比 有序的多分类模型
MASS::polr()
```

## 比例检验 {#sec-prop-test}

### 单样本检验

比例检验函数 `prop.test()` 检验比例是否等于给定的值。单样本的比例检验结果中比例的区间估计与 Wilson 区间估计 [@Wilson1927] 是相关的。区间估计与假设检验是有紧密关系的，对于二项分布比例的 11 种区间估计方法的比较 [@Newcombe1998]。

#### 近似检验

#### 精确检验

函数 `binom.test()` 来做二项检验，函数 `binom.test()` 用来检验伯努利试验中成功概率 $p$ 和给定概率 $p_0$ 的关系，属于精确检验 [@Clopper1934]。

比例 $p$ 的检验，做 $n$ 次独立试验，样本 $X_1,\ldots,X_n \sim b(1, p)$，事件发生的总次数 $\sum_{i=1}^{n}X_i$。

```{r}
# 模拟一组样本
set.seed(20232023)
x <- sample(x = c(0, 1), size = 100, replace = TRUE, prob = c(0.8, 0.2))
```

二项分布中成功概率的检验

```{r}
binom.test(sum(x), n = 100, p = 0.5)
```

检验成功概率 p 是否等于 0.5， P 值 $5.514 \times 10^{-8}$ 结论是拒绝原假设

```{r}
binom.test(sum(x), n = 100, p = 0.2)
```

检验成功概率 p 是否等于 0.2， P 值 0.4534 结论是不能拒绝原假设

切比雪夫不等式（Chebyshev, 1821-1894）。设随机变量 $X$ 的数学期望和方差都存在，则对任意常数 $\epsilon > 0$，有

$$
\begin{aligned}
P(|X - EX| \geq \epsilon) & \leq \frac{Var(X)}{\epsilon^2} \\
P(|X - EX| \leq \epsilon) & \geq 1 - \frac{Var(X)}{\epsilon^2}
\end{aligned}
$$

### 两样本检验

关于两样本的比例检验问题

$$
\begin{aligned}
H_0: P_A = P_B \quad vs. \quad H_1: P_A > P_B \\
H_0: P_A = P_B \quad vs. \quad H_1: P_A < P_B
\end{aligned}
$$

$H_0$ 成立的情况下，暗示着两个样本来自同一总体。

比例检验函数 `prop.test()` 用来检验两组或多组二项分布的成功概率（比例）是否相等。

设随机变量 X 服从参数为 $p$ 的二项分布 $b(n, p)$， $Y$ 服从参数为 $\theta$ 的二项分布 $b(m,\theta)$， $m,n$ 都假定为较大的正整数，检验如下问题

$$
H_0: P_A \geq P_B \quad vs. \quad H_1: P_A < P_B
$$

根据中心极限定理

$$
\frac{\bar{X} - \bar{Y}}{\sqrt{\frac{p(1-p)}{n} + \frac{\theta(1-\theta)}{m}}}
$$

近似服从标准正态分布 $N(0,1)$。如果用矩估计 $\bar{X}$ 和 $\bar{Y}$ 分别替代总体参数 $p$ 和 $\theta$，构造检验统计量

$$
T = \frac{\bar{X} - \bar{Y}}{\sqrt{\frac{\bar{X}(1-\bar{X})}{n} + \frac{\bar{Y}(1-\bar{Y})}{m}}}
$$

根据 Slutsky 定理，检验统计量 $T$ 近似服从标准正态分布，当 $T$ 偏大时，拒绝 $H_0$。该方法的优势在于当 $n,m$ 比较大时，二项分布比较复杂，无法建立统计表，利用标准正态分布表来给出检验所需要的临界值，简便易行！

当 $p$ 和 $\theta$ 都比较小，上述方法检验效果不好，原因在于由中心极限定理对 $\bar{X}$ 和 $\bar{Y}$ 的正态分布近似效果不好，或者间接地导致 $\bar{X}-\bar{Y}$ 的方差偏小，进而 $T$ 的分辨都不好，而且当 $p,\theta$ 很接近 1 时，上述现象也会产生！

下面介绍新的解决办法，办法来自两个二项总体成功概率的比较 [@Song2011]。

上面的检验问题等价于

$$
H_0: \frac{P_A}{P_B} \geq 1 \quad vs. \quad H_1: \frac{P_A}{P_B} < 1
$$

引入检验统计量

$$
T^{\star} = \frac{\bar{X}}{\bar{Y}}
$$

同样由 Slutsky 定理和中心极限定理可知， $\bar{X}/\bar{Y}$ 近似服从 正态分布 $\mathcal{N}(1,\frac{1-\theta}{m\theta})$

当 $(T^\star - 1)/\hat\sigma$ 偏大时接受 $H_0$，临界值可通过 $\mathcal{N}(0, \hat\sigma^2)$ 分布表计算得到， $\hat\sigma^2$ 是对 $\frac{1-\theta}{m\theta}$ 的估计，比如取 $\hat\sigma^2 = \frac{1-\bar{Y}}{m}\cdot \frac{1}{\bar{Y}}$ 或取 $\hat\sigma^2 = \frac{1-\bar{Y}}{m}\cdot \frac{1}{\bar{X}}$

由于渐近方差形如 $\frac{1-\theta}{m\theta}$，因而在 $\theta$ 较小，渐近方差较大，克服了之前 $\bar{X} - \bar{Y}$的方差较小的问题

$p,\theta$ 很接近 1 时，我们取检验统计量

$$
T^{\star\star} = \frac{1-\bar{Y}}{1-\bar{X}}
$$

结论和 $T^\star$ 类似，当 $T^{\star\star}$ 偏大时，拒绝 $H_0$。

### 多样本检验

#### 比例齐性检验

对多组数据的比例检验，可以理解为比例齐性检验。

#### 比例趋势检验

比例趋势检验函数 `prop.trend.test()` 的原假设：四个组里面病人中吸烟的比例是相同的。备择假设：四个组的吸烟比例是有趋势的。

$$
\begin{aligned}
& H_0: P_1 = P_2 = P_3 = P_4 \\
& H_1: P_1 < P_2 < P_3 < P_4 ~\text{或者}~ P_1 > P_2 > P_3 > P_4
\end{aligned}
$$

```{r}
smokers <- c(83, 90, 129, 70)
patients <- c(86, 93, 136, 82)
prop.test(smokers, patients)
prop.trend.test(smokers, patients)
```

## 泊松检验 {#sec-poisson-test}

泊松分布是 1837年由法国数学家泊松 （Poisson, 1781-1840） 首次提出。

$$
p(x) = \frac{\lambda^x\exp(-\lambda)}{x!}, x = 0, 1, \cdots .
$$

泊松分布的期望和方差都是 $\lambda$ ，一般要求 $\lambda > 0$。

### 单样本

`poisson.test()` 泊松分布的参数 $\lambda$ 的精确检验，适用于单样本和两样本。

```{r}
#| eval: false
#| echo: true

poisson.test(x,
  T = 1, r = 1,
  alternative = c("two.sided", "less", "greater"),
  conf.level = 0.95
)
```

参数 `T` 数据的时间单位

### 两样本

## 列联表描述

泰坦尼克号乘客生存死亡统计数据，Titanic 数据集

```{r}
Titanic
```

### 行列分组表格

```{r}
#| label: tbl-titanic-data
#| tbl-cap: "泰坦尼克号乘客生存死亡统计数据"
#| code-fold: true
#| echo: true

# 长格式转宽格式
titanic_data <- reshape(
  data = as.data.frame(Titanic), direction = "wide",
  idvar = c("Class", "Sex", "Age"),
  timevar = "Survived", v.names = "Freq", sep = "_"
)

# 制作表格
gt::gt(titanic_data) |> 
  gt::cols_label(
    Freq_Yes = "存活",
    Freq_No = "死亡",
    Class = "船舱",
    Sex = "性别",
    Age = "年龄"
  )
```

### 百分比堆积图

泰坦尼克号处女航乘客数量按船舱、性别、年龄和存活情况分层， [**ggstats**](https://github.com/larmarange/ggstats/) 包绘制百分比堆积柱形图展示多维分类数据。

```{r}
#| label: fig-titanic-ggstats
#| fig-cap: "百分比堆积柱形图展示多维分类数据"
#| fig-showtext: true
#| fig-width: 7
#| fig-height: 5
#| code-fold: true
#| echo: true

library(ggplot2)
library(ggstats)
ggplot(as.data.frame(Titanic)) +
  aes(x = Class, fill = Survived, weight = Freq, by = Class) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::label_percent()) +
  geom_text(stat = "prop", position = position_fill(.5)) +
  facet_grid(~Sex) +
  labs(x = "船舱", y = "比例", fill = "存活")
```

**ggstats** 包提供的图层 `stat_prop()` 是 `stat_count()` 的变种， `as.data.frame(Titanic)` 中 Age 一列会自动聚合吗？ by = Class 按 Class 分组聚合，统计 Survived 的比例，提供 prop 计算的变量，传递给 `geom_text()` 以添加注释，position 设置将注释放在柱子的中间

### 桑基图

用 [**ggalluvial**](https://github.com/corybrunson/ggalluvial/) 包[@Brunson2020]绘制桑基图展示多维分类数据。

```{r}
#| label: fig-titanic-alluvial
#| fig-cap: "桑基图展示多维分类数据"
#| fig-showtext: true
#| fig-width: 7
#| fig-height: 5
#| code-fold: true
#| echo: true

library(ggplot2)
library(ggalluvial)
ggplot(
  data = as.data.frame(Titanic),
  aes(axis1 = Class, axis2 = Sex, axis3 = Age, y = Freq)
) +
  scale_x_discrete(limits = c("Class", "Sex", "Age")) +
  geom_alluvium(aes(fill = Survived)) +
  geom_stratum() +
  geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
  theme_classic() +
  labs(
    x = "分层维度", y = "人数", fill = "存活",
    title = "泰坦尼克号处女航乘客分层情况"
  )
```

### 马赛克图

```{r}
#| label: fig-titanic-mosaicplot
#| fig-cap: "马赛克图展示多维分类数据"
#| fig-showtext: true
#| fig-width: 7
#| fig-height: 5
#| code-fold: true
#| collapse: true
#| echo: true

op <- par(mar = c(2.5, 2.5, 1.5, 0.5))
mosaicplot(~ Class + Sex + Age + Survived,
  data = Titanic, # shade = TRUE, 
  color = TRUE, border = "white",
  xlab = "船舱", ylab = "性别", main = "泰坦尼克号")
par(op)
```

[**vcd**](https://cran.r-project.org/package=vcd) 包针对分类数据做了很多专门的可视化工作，内置了很多数据集和绘图函数，在 Base R 绘图基础上，整合了许多统计分析功能，提供了一个统一的可视化框架[@Meyer2006; @Zeileis2007]，更多细节见著作《Discrete Data Analysis with R: Visualization and Modeling Techniques for Categorical and Count Data》及其附带的 R 包 [**vcdExtra**](https://github.com/friendly/vcdExtra)[@Friendly2016]。

```{r}
#| label: fig-titanic-mosaic
#| fig-cap: "马赛克图展示多维分类数据"
#| fig-showtext: true
#| fig-width: 7
#| fig-height: 6
#| code-fold: true
#| echo: true

library(grid)
library(vcd)
mosaic(~ Class + Sex + Age + Survived,
  data = Titanic, shade = TRUE, legend = TRUE
)
```

## 列联表分析 {#sec-chisq-test}

是否应该按照列联表的维度分类？还是应该从分析的目的和作用出发？比如我的目的是检验独立性。二者似乎也并不冲突。

列联表中的数据服从多项分布，关于独立性检验，有如下几种常见类型：

1.  相互独立 Mutual independence 所有变量之间相互独立，$X \perp Y \perp Z$ 。
2.  联合独立 Joint independence 两个变量的联合与第三个变量独立，$XY \perp Z$ 。
3.  边际独立 Marginal independence 当忽略第三个变量时，两个变量是独立的。列联表压缩
4.  条件独立 Conditional independence 当固定第三个变量时，两个变量是独立的，$X \perp Y | Z$。

本节数据来自著作《An Introduction to Categorical Data Analysis》[@Agresti2007] 的第2章习题 2.33，探索 1976-1977 年美国佛罗里达州的凶杀案件中被告肤色和死刑判决的关系。

```{r}
#| label: tbl-florida-ethnicity
#| tbl-cap: "佛罗里达州的凶杀案件统计数据"
#| code-fold: true
#| echo: !expr knitr::is_html_output()

tbl <- expand.grid(
  Death = c("Yes", "No"), # 判决结果 是否死刑
  Defend = c("白人", "黑人"),  # 被告 肤色
  Victim = c("白人", "黑人")   # 原告 （被害人）肤色
)
ethnicity <- data.frame(tbl, Freq = c(19, 132, 11, 52, 0, 9,  6, 97))

# 长格式转宽格式
dat1 <- reshape(
  data = ethnicity, direction = "wide",
  idvar = c("Defend", "Victim"),
  timevar = "Death", v.names = "Freq", sep = "_"
)
# 制作表格
gt::gt(dat1) |> 
  gt::cols_label(
    Freq_Yes = "是",
    Freq_No = "否",
    Victim = "被害人",
    Defend = "被告"
  ) |> 
  gt::tab_spanner(
    label = "死刑",
    columns = c(Freq_Yes, Freq_No)
  ) |> 
  gt::opt_row_striping()
```

### 相互独立性

皮尔逊卡方检验（ Pearson's $\chi^2$ 检验） `chisq.test()` 常用于列联表独立性检验和方差分析模型的拟合优度检验。下面是一个 $2 \times 2$ 的列联表。

|        | 第一列 | 第二列 | 合计      |
|--------|--------|--------|-----------|
| 第一行 | $a$    | $b$    | $a+b$     |
| 第二行 | $c$    | $d$    | $c+d$     |
| 合计   | $a+c$  | $b+d$  | $a+b+c+d$ |

: 卡方独立性检验

```{r}
# Death 死刑与 Defend （被告）独立性检验
m <- xtabs(Freq ~ Death + Defend, data = ethnicity)
m
chisq.test(m, correct = TRUE)
chisq.test(m, correct = FALSE)
```

当被告是白人时，死刑判决 19 个，占总的死刑判决数量的 19/36 = 52.78%，当被告是黑人时，死刑判决 17 个，占总的死刑判决数量的 17/36 = 47.22%。判决结果与被告种族没有显著关系，但与原告（受害人）种族是有关系的，请继续往下看。

```{r}
# Death 死刑与 Victim （原告）独立性检验
m <- xtabs(Freq ~ Death + Victim, data = ethnicity)
chisq.test(m, correct = TRUE)
chisq.test(m, correct = FALSE)
```

当受害人是白人时，死刑判决 30 个，占总的死刑判决数量的 30/36 = 83.33%，当受害人是黑人时，死刑判决 6 个，占总的死刑判决数量的 6/36 = 16.67%。受害人是白人时，死刑判决明显多于黑人。

多维列联表

```{r}
m <- xtabs(Freq ~ Death + Defend + Victim, data = ethnicity)
m
```

判决结果、被告种族、原告种族三者是否存在联合独立性，即考虑 (Victim, Death) 是否与 Defend 独立，(Victim, Defend) 是否与 Death 独立，(Death, Defend) 与 Victim 是否相互独立。

```{r}
fm <- loglin(table = m, margin = list(c(1, 2), c(1, 3), c(2, 3)), print = FALSE)
fm 
# 拟合对数线性模型
# fm <- loglin(m, list(c(1), c(2), c(3)))
# fm
```

似然比检验统计量（Likelihood Ratio Test statistic），皮尔逊 $\chi^2$ 统计量（Pearson X-square Test statistic）

```{r}
1 - pchisq(fm$lrt, fm$df)
```

拟合对数线性模型

```{r}
fit_dvp <- glm(Freq ~ ., data = ethnicity, family = poisson(link = "log"))
```

模型输出

```{r}
summary(fit_dvp)
```

Pearson $\chi^2$ 统计量

```{r}
sum(residuals(fit_dvp, type = "pearson")^2)
```

**MASS** 包计算模型参数的置信区间

```{r}
#| message: false

confint(fit_dvp, trace = FALSE)
```

对于单元格总样本量小于 40 或 T 小于 1 时，需采用费希尔精确检验（ Fisher 's Exact 检验）。

### 边际独立性

费希尔精确检验：固定边际的情况下，检验列联表行和列之间的独立性 `fisher.test()` 。

`fisher.test()` 函数用法，统计原理和公式，适用范围和条件，概念背景和历史。

费舍尔 (Sir Ronald Fisher, 1890.2 -- 1962.7)[^categorical-data-analysis-1] 和一位女士打赌，女士说能品出奶茶中奶和茶的添加顺序。

[^categorical-data-analysis-1]: <https://en.wikipedia.org/wiki/Ronald_Fisher>

`fisher.test()` 针对计数数据，检验列联表中行和列的独立性。

```{r}
TeaTasting <- matrix(c(3, 1, 1, 3),
  nrow = 2,
  dimnames = list(
    Guess = c("Milk", "Tea"),
    Truth = c("Milk", "Tea")
  )
)
TeaTasting
```

```{r}
# 单边 P 值
fisher.test(TeaTasting, alternative = "greater")
# 双边 P 值
fisher.test(TeaTasting, alternative = "two.sided")
# 单边 P 值
sum(dhyper(x = c(3, 4), m = 4, n = 4, k = 4))
```

### 对称性

用于计数数据的 McNemar 卡方检验（ McNemar $\chi^2$ 检验）：检验二维列联表行和列的对称性 `mcnemar.test()`。怎么理解对称性？其实是配对检验。看帮助实例。

```{r}
Performance <- matrix(c(794, 86, 150, 570),
  nrow = 2,
  dimnames = list(
    "1st Survey" = c("Approve", "Disapprove"),
    "2nd Survey" = c("Approve", "Disapprove")
  )
)
Performance
mcnemar.test(Performance)
```

### 条件独立性

用于分层分类数据的 Cochran-Mantel-Haenszel 卡方检验：两个枚举（分类）变量的条件独立性，假定不存在三个因素的交互作用。Cochran-Mantel-Haenszel 检验 `mantelhaen.test()`

```{r}
str(UCBAdmissions)
```

`UCBAdmissions` 数据集是一个 $2\times 2 \times 6$ 的三维列联表，R 语言中常用 table 类型表示。实际上，table 类型衍生自 array 数组类型，当把 `UCBAdmissions` 当作一个数组操作时，1、2、3 分别表示 Admit、Gender、Dept 三个维度。

```{r}
mantelhaen.test(UCBAdmissions)
```

没有证据表明院系与性别之间存在关联。在给定院系的情况下，是否录取和性别没有显著关系。

```{r}
# 按系统计
apply(UCBAdmissions, 3, function(x) (x[1, 1] * x[2, 2]) / (x[1, 2] * x[2, 1]))

woolf <- function(x) {
  x <- x + 1 / 2
  k <- dim(x)[3]
  or <- apply(x, 3, function(x) (x[1, 1] * x[2, 2]) / (x[1, 2] * x[2, 1]))
  w <- apply(x, 3, function(x) 1 / sum(1 / x))
  1 - pchisq(sum(w * (log(or) - weighted.mean(log(or), w))^2), k - 1)
}
woolf(UCBAdmissions)
```

## 加州伯克利分校的录取情况 {#sec-ucb-admissions}

1973 年加州伯克利分校 6 个最大的院系的录取情况见下 @tbl-ucb-admissions ，研究目标是加州伯克利分校在招生录取工作中是否有性别歧视？

```{r}
#| label: tbl-ucb-admissions
#| tbl-cap: "加州伯克利分校的录取情况"
#| echo: false

# knitr::kable(as.data.frame(UCBAdmissions),
#   col.names = c("录取与否", "学生性别", "院系", "人数")
# )
# 长格式转宽格式
dat1 <- reshape(
  data = as.data.frame(UCBAdmissions), direction = "wide",
  idvar = c("Dept", "Admit"),
  timevar = "Gender", v.names = "Freq", sep = "_"
)
# 再一次长格式转宽格式
dat2 <- reshape(
  data = dat1, direction = "wide",
  idvar = "Dept", timevar = "Admit",
  v.names = c("Freq_Male", "Freq_Female"), sep = "_"
)
# 面对 HTML 和 PDF 输出 kableExtra 需要两套代码
# kableExtra 不支持 DOCX 输出
# knitr::kable(dat2, booktabs = TRUE, row.names = F, col.names = NULL, align = "ccccc") |>
#   kableExtra::kable_styling(
#     bootstrap_options = "basic", full_width = TRUE, position = "center"
#   ) |> 
#   kableExtra::add_header_above(c("院系" = 1, "男性" = 1, "女性" = 1, "男性" = 1, "女性" = 1)) |>
#   kableExtra::add_header_above(c(" " = 1, "录取" = 2, "拒绝" = 2))

gt::gt(dat2) |> 
  gt::cols_label(
    Dept = "院系",
    Freq_Male_Admitted = "男性",
    Freq_Female_Admitted = "女性",
    Freq_Male_Rejected = "男性",
    Freq_Female_Rejected = "女性"
  ) |> 
  gt::tab_spanner(
    label = "录取",
    columns = c(Freq_Male_Admitted, Freq_Female_Admitted)
  ) |> 
  gt::tab_spanner(
    label = "拒绝",
    columns = c(Freq_Male_Rejected, Freq_Female_Rejected)
  ) |> 
  gt::opt_row_striping()
```

借助马赛克图 @fig-ucb-admissions-mosaicplot 可以更加直观的看出数据中的比例关系。

```{r}
#| label: fig-ucb-admissions-mosaicplot
#| fig-width: 7
#| fig-height: 5
#| fig-cap: "加州伯克利分校院系录取情况"
#| fig-showtext: true
#| echo: false

# plot(UCBAdmissions, col = "lightblue", border = "white",
#      main = "", xlab = "性别", ylab = "院系")
op <- par(mar = c(0.5, 2, 0.5, 0.5))
mosaicplot(~ Admit + Dept + Gender,
  data = UCBAdmissions, color = TRUE, border = "white",
  main = "", xlab = "", ylab = "院系"
)
on.exit(par(op), add = TRUE) 
```

接下来进行定量的分析，首先，按性别和录取情况统计人数，如下：

```{r}
m <- xtabs(Freq ~ Gender + Admit, data = as.data.frame(UCBAdmissions))
m
```

可以看到，申请加州伯克利分校的女生当中，只有 $557 / (557 + 1278) = 30.35\%$ 录取了，而男生则有 $1198 / (1198 + 1493) = 44.52\%$ 的录取率。根据皮尔逊 $\chi^2$ 检验：

```{r}
# 不带耶茨矫正
chisq.test(m, correct = FALSE)
```

可知 $\chi^2$ 统计量的值为 $92.205$ 且 P 值远远小于 0.05， 差异达到统计显著性，不是随机因素导致的。因此，加州伯克利分校被指控在招生录取工作中存在性别歧视。然而，当我们细分到各个院系去看录取率（录取人数 / 申请人数），结果显示院系 A 的录取率为 64.41%，院系 B 的录取率为 63.24%，依次类推，各院系情况如下：

```{r}
proportions(xtabs(Freq ~ Dept + Admit,
  data = as.data.frame(UCBAdmissions)
), margin = 1)
```

```{r}
#| label: fig-ucb-admissions-fourfoldplot
#| fig-width: 7
#| fig-height: 5
#| fig-cap: "加州伯克利分校各院系录取情况"
#| fig-showtext: true
#| echo: false
#| par: true

fourfoldplot(aperm(UCBAdmissions, c(2, 1, 3)), mfcol = c(2, 3))
```

对每个院系，单独使用皮尔逊 $\chi^2$ 检验，发现只有 A 系的男、女生录取率的差异达到统计显著性，其它系的差异都不显著。辛普森悖论在这里出现了，在分类数据的分析中，常常遇到。

```{r}
# 以 A 系为例
ma <- xtabs(Freq ~ Gender + Admit,
  subset = Dept == "A",
  data = as.data.frame(UCBAdmissions)
)
chisq.test(ma, correct = FALSE)
```

为了经一步说明此现象的原因，建立对数线性模型来拟合数据，值得一提的是皮尔逊卡方检验可以从对数线性模型的角度来看，而对数线性模型是一种特殊的广义线性模型，针对计数数据建模。

```{r}
fit_ucb0 <- glm(Freq ~ Dept + Admit + Gender,
  family = poisson(link = "log"),
  data = as.data.frame(UCBAdmissions)
)
summary(fit_ucb0)
```

添加性别和院系的交互效应后，对数线性模型的 AIC 下降一半多，说明模型的交互效应是显著的，也就是说性别和院系之间存在非常强的关联。

```{r}
fit_ucb1 <- glm(Freq ~ Dept + Admit + Gender + Dept * Gender,
  family = poisson(link = "log"),
  data = as.data.frame(UCBAdmissions)
)
summary(fit_ucb1)
```

此辛普森悖论现象的解释是女生倾向于申请录取率低的院系，而男生倾向于申请录取率高的院系，最终导致整体上，男生的录取率显著高于女生。至于为什么女生会倾向于申请录取率低的院系？这可能要看具体的院系是哪些，招生政策如何？这已经不是仅仅依靠招生办的统计数字就可以完全解释得了的，更多详情见文献 @Bickel1975 。

::: callout-tip
对数线性模型的皮尔逊 $\chi^2$ 检验的统计量

```{r}
sum(residuals(fit_ucb1, type = "pearson")^2)
```

比较多个广义线性模型的拟合效果，除了看 AIC，还可以看对数似然，它越大越好。可以看到添加性别和院系的交互效应后，对数似然增加了一倍多。

```{r}
# 基础模型
logLik(fit_ucb0)
# 添加交互效应
logLik(fit_ucb1)
```

```{r}
#| eval: false
#| echo: false

# 似然比、得分和卡方检验，三大检验在此处等价
anova(fit_ucb1, test = "LRT")
anova(fit_ucb1, test = "Rao")
anova(fit_ucb1, test = "Chisq")
anova(fit_ucb0, fit_ucb1, test = "Chisq")
# 似然比是渐近卡方分布的
anova(fit_ucb0, fit_ucb1, test = "LRT")
```
:::

```{r}
#| eval: false
#| echo: false

# https://github.com/nhs-r-community/FunnelPlotR
library(FunnelPlotR)
library(ggplot2)

ucb_dat <- as.data.frame(UCBAdmissions)

fit_ucb <- glm(Freq ~ Dept + Admit + Gender,
  family = poisson(link = "log"),
  data = ucb_dat
)
summary(fit_ucb)

ucb_dat$preds <- predict(fit_ucb, type = "response")

funnel_plot(
  numerator = ucb_dat$Freq,
  denominator = ucb_dat$preds, # 预测值
  group = ucb_dat$Gender, # 分组
  title = "Length of Stay Funnel plot for `medpar` data",
  data_type = "SR", # 标准 standard ratio
  limit = 99, # 置信限
  draw_unadjusted = TRUE,
  draw_adjusted = FALSE, label = "outlier"
)
# 泊松分布的均值等于方差，如果不是，则意味着 overdispersion
# 方差是均值的 125 倍
sum(fit_ucb$weights * fit_ucb$residuals^2) / fit_ucb$df.residual
```

## 分析泰坦尼克号乘客生存率 {#sec-titanic}

分析存活率的影响因素。

除了从条件独立性检验的角度，下面从逻辑回归模型的角度分析这个高维列联表数据，由此，我们可以知道假设检验和广义线性模型之间的联系，针对复杂高维列联表数据进行关联分析和解释。

响应变量是乘客的状态，存活还是死亡，titanic_data 是按船舱 Class、性别 Sex 和年龄 Age 分类汇总统计的数据，因此，下面的逻辑回归模型是对乘客群体的建模。

```{r}
# 建立模型
fit_titanic <- glm(cbind(Freq_Yes, Freq_No) ~ Class + Sex + Age,
  data = titanic_data, family = binomial(link = "logit")
)
```

接着，我们查看模型输出的情况

```{r}
# 模型输出
summary(fit_titanic)
```