第十章 方差分析与卡方检验
蒋新宇
2025-12-31
10.1 方差分析
10.1.1 方差分析的基本思想
上一章中，我们通过例 9-5 了解了方差比的 F 检验。F 检验还有一种更为广泛的应用，即比较多个总体均值是否存在显著差异，这在分析分类自变量对数值因变量的影响时尤为常用。这种方法更广为人知的名称是方差分析 (ANOVA)。方差分析的核心在于通过分解数据的变异来源，间接推断总体均值是否存在差异。这与 Z 检验、t 检验等直接比较均值的方法有着本质的区别。
为确保方差分析结论的统计有效性，该方法建立在以下三个基本假设之上：
正态性：每个总体应服从或近似服从正态分布。这是因为 F 统计量的抽样分布是基于正态总体推导得出的，只有满足此假设，检验结果才具有理论效力。
同方差性：各总体的方差应相等（即具有方差齐性）。这是比较均值的前提，因为只有当组内的波动水平保持一致时，若检测到组间存在显著差异，我们才有理由将其归因于均值的差异，而非方差的差异。
独立性：各观测值之间必须相互独立。该假设确保了每个数据点提供的信息不重叠，是进行所有统计推断的基础。
方差分析的核心思想在于对数据的总变异进行分解，将其归属于不同的变异来源（主要为组间变异和组内变异），并通过比较这两部分变异的相对大小进行统计推断。若组间变异显著大于组内变异，则有理由推断不同组别的均值存在显著差异。（在方差分析的语境中，“变异”指代数值的波动程度，其量化指标为“平方和”。）以下将详细阐述变异的分解原理及计算公式。
- 总变异
总变异反映了所有观测值围绕总均值的离散程度，采用总平方和 (SST)进行度量，其计算公式为所有观测值与总均值之差的平方和：
\[SST = \sum_{i=1}^{k} \sum_{j=1}^{n_i} (x_{ij} - \bar{\bar{x}})^2\]
其中，\(k\) 为组数，\(n_i\) 为第 \(i\) 组的样本量，\(x_{ij}\) 表示第 \(i\) 组的第 \(j\) 个观测值，\(\bar{\bar{x}}\) 为总均值。
- 组间变异
组间变异反映了不同处理（或水平）对观测值产生的系统性影响，采用组间平方和 (SSA) 度量，即各组均值与总均值之差的加权平方和：
\[SSA = \sum_{i=1}^{k} n_i (\bar{x_i} - \bar{\bar{x}})^2\]
其中，\(\bar{x_i}\) 为第 \(i\) 组的组内均值。
- 组内变异
组内变异反映了组内随机误差或个体间的差异，采用组内平方和 (SSE) 度量，即组内各观测值与本组均值之差的平方和。其数值等于总平方和减去组间平方和：
\[SSE = \sum_{i=1}^{k} \sum_{j=1}^{n_i} (x_{ij} - \bar{x_i})^2 = SST - SSA\]
在此基础上，方差分析通过比较组间变异与组内变异的相对大小来进行推断。值得注意的是，直接比较平方和是不严谨的，因为平方和的大小会受到样本量（更准确地说是自由度）的影响。数据越多，平方和累积往往越大。为了消除这一影响，必须将平方和转化为平均变异，即“均方”。
均方的计算方法如下：
组间均方 (MSA)：\(MSA = SSA / df_A\)，其中组间自由度 \(df_A = k - 1\)。
组内均方 (MSE)：\(MSE = SSE / df_E\)，其中组内自由度 \(df_E = N - k\)，\(N\) 为总样本量。
最终，构建方差分析的F统计量作为检验标准：
\[F = \frac{MSA}{MSE} = \frac{SSA / (k-1)}{SSE / (N-k)}\]
在原假设（各总体均值相等）成立的条件下，该 F 统计量服从F分布，即 \(F \sim F(df_A, df_E)\)。通过将计算得到的 F 值与 F 分布的临界值进行比较，或计算相应的p值，即可判断组间差异是否具有统计学意义。
这种基于方差比值的间接推断方法，使得方差分析在处理多组均值比较时，相较于多次两两t 检验具有显著优势，能够有效控制第一类错误的膨胀。其原因在于，进行多次两两 t 检验会累积犯“假阳性”错误（弃真）的风险；而方差分析通过一次整体性检验，能够将总的犯错误概率严格控制在预设的显著性水平（如 5%）之内。
10.1.2 单因素方差分析
单因素方差分析是一种用于检验一个自变量（称为“因素”）的不同类别（称为“水平”）是否对某个连续型因变量的均值产生显著影响的统计方法。其分析过程遵循上述标准的统计假设检验流程，这里我们采用方法分析表，进行简化：
- 构建原假设和备择假设。
- 原假设 (H₀)：所有水平对应的总体均值均相等，即因素对因变量无显著影响。公式表示为：μ₁ = μ₂ = … = μₖ（k为水平数）。
- 备择假设 (H₁)：至少有两个水平的总体均值不相等，即因素对因变量有显著影响。
- 依照上一节的公式，计算方差分析所需要的各指标以及F统计量。进而汇总到方差分析表中：
| 变异来源 | 平方和 (SS) | 自由度 (df) | 均方 (MS) | F值 |
|---|---|---|---|---|
| 组间 | SSA | k-1 | MSA = SSA/(k-1) | F = MSA / MSE |
| 组内 | SSE | N-k | MSE = SSE/(N-k) | |
| 总计 | SST | N-1 |
- 统计决策与结论 在零假设成立的条件下，F统计量服从自由度为(df₁ = k-1, df₂ = N-k)的F分布。将计算得到的F值与选定显著性水平α（如0.05）下的F分布临界值进行比较：
- 若 F值 > F临界值，则拒绝H₀，认为因素的不同水平对因变量均值有显著影响。
- 若 F值 ≤ F临界值，则没有充分证据拒绝H₀。
- 事后检验 若方差分析的结果显著（拒绝H₀）。但这个结论仅表明至少存在两个水平的均值差异显著，但无法指明具体是哪几组之间存在差异。此时，可以进行进一步的事后检验，如LSD检验、Tukey HSD检验等，以进行两两之间的详细比较。
例10-1：
某研究机构希望评估三种不同的台风预警信息发布策略（因素：预警策略）对沿海城市居民应急物资储备天数（因变量：储备天数）的影响。
水平A（基础文本）：通过短信和广播发布简单的台风路径和风力预报。
水平B（多媒体增强）：在A基础上，增加图文推送和短视频，直观展示影响范围和避险指南。
水平C（社区联动）：在B基础上，嵌入社区干部上门讲解和邻里互助动员。
我们想知道，不同的台风预警策略，是否会对居民的应急物资储备天数产生显著影响？
收集数据如下：
| 水平A (基础文本) | 水平B (多媒体增强) | 水平C (社区联动) |
|---|---|---|
| 3 | 5 | 7 |
| 4 | 6 | 8 |
| 5 | 7 | 6 |
| 4 | 5 | 9 |
| 3 | 7 | 8 |
| 6 | ||
| \(n_A = 5\) | \(n_B = 6\) | \(n_C = 5\) |
| \(\bar{x}_A = 3.8\) | \(\bar{x}_B = 6.0\) | \(\bar{x}_C = 7.6\) |
- 构建研究假设
H₀: μ_A = μ_B = μ_C （三种策略下的平均储备天数无显著差异）
H₁: 至少有两种策略下的平均储备天数存在显著差异。
- 计算方差分析表：
总样本量 \(N = 5 + 6 + 5 = 16\)，组数 \(k = 3\)。
总均值 \(\bar{\bar{x}} = (3+4+5+...+8+6+9+8) / 16 = 91 / 16 = 5.6875\)
计算平方和 (SS)：
总平方和 (SST)：所有观测值与总均值之差的平方和。
\(SST = (3-5.6875)^2 + (4-5.6875)^2 + ... + (8-5.6875)^2 = 54.4375\)
组间平方和 (SSA)：各组均值与总均值之差的加权平方和。
\(SSA = 5*(3.8-5.6875)^2 + 6*(6.0-5.6875)^2 + 5*(7.6-5.6875)^2 = 36.718\)
组内平方和 (SSE)：总平方和减去组间平方和。
\(SSE = SST - SSA = 54.4375 - 36.718 = 17.7195\)
确定自由度 (df)：
组间自由度：\(df_A = k - 1 = 3 - 1 = 2\)
组内自由度：\(df_E = N - k = 16 - 3 = 13\)
总自由度：\(df_T = N - 1 = 15\)
计算均方 (MS)：
组间均方：\(MSA = SSA / df_A = 36.718 / 2 = 18.359\)
组内均方：\(MSE = SSE / df_E = 17.7195 / 13 ≈ 1.363\)
计算F统计量：
\(F = \frac{MSA}{MSE} = \frac{18.359}{1.363} ≈ 13.47\)
方差分析表：
| 变异来源 | 平方和 (SS) | 自由度 (df) | 均方 (MS) | F值 | P值 |
|---|---|---|---|---|---|
| 组间（策略差异） | 36.718 | 2 | 18.359 | 13.47 | 0.0001 |
| 组内（随机误差） | 17.720 | 13 | 1.363 | ||
| 总计 | 54.438 | 15 |
- 统计决策和结论
设定α = 0.05，查F分布表或者通过软件计算得临界值 \(F_{0.05}(2, 13) ≈ 3.81\)。 由于计算F值 13.47 > 3.81，因此拒绝零假设(H₀)。因此，可以得出结论：在0.05的显著性水平上，有充分证据表明，不同的台风预警策略对居民的应急物资储备天数产生了显著影响。
- 事后检验（如果我们想明确差异来源可以继续做LSD或者Tukey HSD。我们这里省略。只展示结果。事后检验表明：
策略B（多媒体增强）和策略C（社区联动）的效果均显著优于策略A（基础文本）。
策略B与策略C之间的效果差异在统计学上不显著。
因此，我们可以利用这一结果进行防灾减灾规划：将预警信息从“基础文本”升级为“多媒体增强”能显著提升民众准备程度；而进一步投入大量人力进行“社区联动”，其附加效益在此研究中并不明显。我们可据此优化投入，优先推广“多媒体增强”策略。
# 加载必要的包
library(ggplot2)
library(dplyr)
# 创建数据框
strategy_data <- data.frame(
strategy = factor(rep(c("基础文本", "多媒体增强", "社区联动"),
times = c(5, 6, 5))),
days = c(3, 4, 5, 4, 3, # 策略A
5, 6, 7, 5, 7, 6, # 策略B
7, 8, 6, 9, 8) # 策略C
)
# 单因素方差分析
anova_result <- aov(days ~ strategy, data = strategy_data)
{
# 单因素方差分析结果
cat("=====方差分析表=====\n")
print(summary(anova_result))
# 事后检验
tukey_result <- TukeyHSD(anova_result)
cat("=====事后检验=====\n")
print(tukey_result$strategy)
}## =====方差分析表=====
## Df Sum Sq Mean Sq F value Pr(>F)
## strategy 2 36.44 18.219 19.74 0.000115 ***
## Residuals 13 12.00 0.923
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## =====事后检验=====
## diff lwr upr p adj
## 基础文本-多媒体增强 -2.2 -3.73613892 -0.6638611 6.022315e-03
## 社区联动-多媒体增强 1.6 0.06386108 3.1361389 4.096204e-02
## 社区联动-基础文本 3.8 2.19555523 5.4044448 8.183061e-05# 箱线图+均值点+显著性标记
ggplot(strategy_data, aes(x = strategy, y = days, fill = strategy)) +
geom_boxplot(alpha = 0.7, width = 0.5) +
geom_jitter(width = 0.2, size = 2, alpha = 0.8) +
stat_summary(fun = mean, geom = "point", shape = 23, size = 4,
fill = "red", color = "black") +
labs(subtitle = paste("不同预警策略对应急物资储备天数影响的单因素方差分析: F =",
round(summary(anova_result)[[1]][1, "F value"], 2),
", p =", round(summary(anova_result)[[1]][1, "Pr(>F)"], 4)),
x = "预警策略",
y = "应急物资储备天数（天）") +
scale_fill_manual(values = c("#FF9999", "#99CCFF", "#99FF99")) +
theme_minimal() +
theme(legend.position = "none",
plot.title = element_text(face = "bold"))
10.1.3 双因素方差分析
1. 双因素方差分析的基本原理
双因素方差分析同时考察两个因素（因素 A 和因素 B）对因变量的影响。在单因素方差分析的基础上，其核心概念如下：
主效应：指因素 A 或因素 B 单独独立产生的影响。例如，不同台风预警策略的效果（因素 A），或者不同类型居住区域的差别（因素 B）。
交互作用：指因素 A 和因素 B 联合起来产生的影响。例如，某策略在老旧小区特好用，但在新建小区效果一般，这种“因人而异”或“因地制宜”的效果即是交互作用。
变异分解：
\[ SST = SSA + SSB + SSAB + SSE \]
（总变异 = A因素变异 + B因素变异 + 交互变异 + 随机误差）
根据分析中是否包含重复观测以及数据的平衡性，双因素方差分析可以进一步细分为无重复双因素方差分析、有重复双因素方差分析以及非平衡数据的方差分析。
2. 无重复双因素方差分析
如果在两个因素之间不存在交互作用，或者交互作用可以忽略不计，我们可以采用无重复双因素方差分析。这里的“无重复”是指在每个具体的“策略 × 小区”组合中，只做一次实验（\(n=1\)），只有一个观测值。
由于 \(n=1\)，组内没有波动，我们无法像单因素方差分析那样直接计算出组内纯粹的随机误差（\(SSE\)）。因此，在数学处理上，我们假设交互作用不存在，将原本代表交互效应的变异（\(SSAB\)）归入误差项（\(SSE\)）。所以，无重复双因素方差分析只能检验 A 和 B 的主效应，无法检验交互作用。
例10-2：
某研究机构探究 3 种台风预警发布策略（因素A）在 2 类居住区域（因素B）下，居民应急物资储备天数的差异。假设策略效果不受区域类型影响（无交互），每个组合仅抽取 1 个代表性样本观测，数据如下：
| B₁（老旧小区） | B₂（新建小区） | 行平均 \(\bar{x}_{i.}\) | |
|---|---|---|---|
| A₁（基础文本） | 4 | 3 | 3.5 |
| A₂（多媒体增强） | 7 | 5 | 6.0 |
| A₃（社区联动） | 9 | 7 | 8.0 |
| 列平均 \(\bar{x}_{.j}\) | 6.67 | 5.0 | 总均值 \(\bar{\bar{x}} = 5.83\) |
- 研究假设
\(H_{0A}\)：\(\mu_{A1}=\mu_{A2}=\mu_{A3}\)，不同预警策略对储备天数无显著影响；\(H_{1A}\)：至少两种策略的储备天数存在显著差异
\(H_{0B}\)：\(\mu_{B1}=\mu_{B2}\)，不同居住区域对储备天数无显著影响；\(H_{1B}\)：不同居住区域的储备天数存在显著差异
- 计算平方和与自由度
\(r=3\)（A的水平数），\(c=2\)（B的水平数），总观测\(N=rc=6\)
总平方和 SST：
\[ \begin{aligned} SST&=\sum_{i=1}^r\sum_{j=1}^c (x_{ij}-\bar{\bar{x}})^2\\ &=(4-5.83)^2+(3-5.83)^2+(7-5.83)^2+(5-5.83)^2+(9-5.83)^2+(7-5.83)^2\\ &=3.35+8.01+1.37+0.69+10.05+1.37= \boldsymbol{24.84} \end{aligned} \]
因素A（策略）平方和 SSA：
\[ \begin{aligned} SSA &= c\sum_{i=1}^r (\bar{x}_{i.}-\bar{\bar{x}})^2\\ &=2\times\left[(3.5-5.83)^2+(6-5.83)^2+(8-5.83)^2\right]\\ &=2\times(5.43+0.03+4.70)=\boldsymbol{20.32} \end{aligned} \]
因素B（区域）平方和 SSB：
\[ \begin{aligned} SSB &= r\sum_{j=1}^c (\bar{x}_{.j}-\bar{\bar{x}})^2\\ &=3\times\left[(6.67-5.83)^2+(5.0-5.83)^2\right]\\ &=3\times(0.71+0.69)=\boldsymbol{4.20} \end{aligned} \]
误差平方和 SSE（无重复设计中，SSE = SST - SSA - SSB，对应无法分离的交互+随机误差）：
\[SSE = 24.84 - 20.32 - 4.20 = \boldsymbol{0.32}\]
自由度：
\(df_A = r-1 = 2\)
\(df_B = c-1 = 1\)
\(df_E = (r-1)(c-1) = 2\)
\(df_T = N-1 = 5\)
计算均方与 F 统计量
均方（MS）等于平方和除以对应的自由度：
\(MS_A = \frac{SSA}{df_A} = \frac{20.32}{2} = 10.16\)
\(MS_B = \frac{SSB}{df_B} = \frac{4.20}{1} = 4.20\)
\(MS_E = \frac{SSE}{df_E} = \frac{0.32}{2} = 0.16\)
F 统计量等于因素均方除以误差均方：
\(F_A = \frac{MS_A}{MS_E} = \frac{10.16}{0.16} = 63.5\)
\(F_B = \frac{MS_B}{MS_E} = \frac{4.20}{0.16} = 26.25\)
方差分析表
| 变异来源 | 平方和SS | 自由度 | 均方MS | F值 |
|---|---|---|---|---|
| 预警策略（A） | 20.32 | 2 | 10.16 | 63.5 |
| 居住区域（B） | 4.20 | 1 | 4.20 | 26.25 |
| 误差 | 0.32 | 2 | 0.16 | - |
| 总计 | 24.84 | 5 | - | - |
- 统计决策与结论
给定显著性水平 \(\alpha=0.05\)，查表得临界值：\(F_{0.05}(2,2)=19.00\)，\(F_{0.05}(1,2)=18.51\)。
策略：\(F_A=63.5 > 19.00\)，拒绝 \(H_0\)。结论：不同预警策略对居民应急物资储备天数有显著影响。
区域B：\(F_B=26.25 > 18.51\)，拒绝 \(H_0\)。结论：老旧小区的储备天数显著高于新建小区。
# 加载必要的包
library(ggplot2)
library(dplyr)
# --- 创建/准备数据 ---
data_10_2 <- data.frame(
strategy = factor(rep(c("A1(基础文本)", "A2(多媒体增强)", "A3(社区联动)"), times = 2)),
region = factor(rep(c("B1(老旧小区)", "B2(新建小区)"), each = 3)),
days = c(4, 7, 9, # B1 老旧小区
3, 5, 7) # B2 新建小区
)
# --- 交互作用图 ---
# 先通过肉眼观察线条是否平行，以判断“无交互作用”的模型假设是否成立。
ggplot(data_10_2, aes(x = strategy, y = days, color = region, group = region)) +
geom_line(size = 1, linetype = "dashed") +
geom_point(size = 4) +
geom_text(aes(label = days), vjust = -0.8, show.legend = FALSE) +
labs(title = "交互作用图",
subtitle = "两条线平行，支持‘无交互作用’的假设",
x = "预警策略 (A)",
y = "应急物资储备天数（天）",
color = "居住区域 (B)") +
theme_minimal()
# ---统计检验 ---
#由于图形显示平行，我们确信可以使用无交互项方差分析。
anova_result_10_2 <- aov(days ~ strategy + region, data = data_10_2)
{
cat("\n=====方差分析表=====\n")
print(summary(anova_result_10_2))
}##
## =====方差分析表=====
## Df Sum Sq Mean Sq F value Pr(>F)
## strategy 2 20.333 10.167 61 0.0161 *
## region 1 4.167 4.167 25 0.0377 *
## Residuals 2 0.333 0.167
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13. 有重复双因素方差分析
如果在每个组合中进行了多次实验（\(n \ge 2\)），即每个格子里有 2 个或更多的数据，这种设计被称为有重复双因素方差分析。由于有重复数据，我们可以计算出组内差异，从而将真正的随机误差（\(SSE\)）从背景噪声中分离出来。这种设计的核心优势在于能够将交互作用（\(SSAB\)）完全剥离出来，并同时检验主效应 A、主效应 B 以及它们的交互作用 A×B。在进行统计决策时，通常需要先观察交互作用是否显著：如果交互作用显著，说明因素之间会相互干扰，此时不能只看主效应，需要单独分析简单效应；如果交互作用不显著，再进一步分析 A 和 B 哪个因素的主效应更显著。
例10-3： 延续前例，现在研究机构在每个组合中抽取 2 个独立样本（重复实验，\(n=2\)），旨在验证“策略效果是否因区域类型不同而变化”（即检验交互作用）。数据如下：
| B₁（老旧小区） | B₂（新建小区） | 行平均 \(\bar{x}_{i..}\) | |
|---|---|---|---|
| A₁（基础文本） | 4, 5 | 2, 3 | 3.5 |
| A₂（多媒体增强） | 6, 7 | 5, 6 | 6.0 |
| A₃（社区联动） | 8, 9 | 6, 7 | 7.5 |
| 列平均 \(\bar{x}_{.j.}\) | 6.5 | 4.83 | 总均值 \(\bar{\bar{x}}=5.67\) |
每个单元格 \(n=2\)，\(r=3\)，\(c=2\)，总观测 \(N=rcn=12\)。
- 研究假设
\(H_{0A}\)：不同预警策略对储备天数无显著影响；\(H_{1A}\)：至少两种策略有差异
\(H_{0B}\)：不同居住区域对储备天数无显著影响；\(H_{1B}\)：两类区域有差异
\(H_{0AB}\)：预警策略与居住区域无交互作用；\(H_{1AB}\)：两者存在交互作用
- 计算平方和与自由度
总平方和 SST：
\[ SST=\sum_{i=1}^r\sum_{j=1}^c\sum_{k=1}^n (x_{ijk}-\bar{\bar{x}})^2 = 43.33 \]
因素A（策略）平方和 SSA：
\[ \begin{aligned} SSA &= nc\sum_{i=1}^r (\bar{x}_{i..}-\bar{\bar{x}})^2\\ &=2\times2\times\left[(3.5-5.67)^2+(6-5.67)^2+(7.5-5.67)^2\right]=32.67 \end{aligned} \]
因素B（区域）平方和 SSB：
\[ \begin{aligned} SSB &= nr\sum_{j=1}^c (\bar{x}_{.j.}-\bar{\bar{x}})^2\\ &=2\times3\times\left[(6.5-5.67)^2+(4.83-5.67)^2\right]=8.33 \end{aligned} \]
误差平方和 SSE：每个单元格内观测对单元格均值的离差平方和。每个单元格内两个数据相差1，离差平方和恒为0.5（即 \((1^2+(-1)^2)/2\)）。共6个单元格。
$$ SSE = 6 \times 0.5 = 3.0 $$交互作用平方和 SSAB：
\[ SSAB = SST - SSA - SSB - SSE = 43.33 - 32.67 - 8.33 - 3.0 \approx 0 \] 自由度：
\(df_A = r-1 = 2\)
\(df_B = c-1 = 1\)
\(df_{AB} = (r-1)(c-1) = 2\)
\(df_E = rc(n-1) = 3\times2\times1 = 6\)
\(df_T = N-1 = 11\)
计算均方与 F 统计量
\(MS_A = SSA / df_A = 32.67 / 2 = 16.335\)
\(MS_B = SSB / df_B = 8.33 / 1 = 8.33\)
\(MS_{AB} = SSAB / df_{AB} = 0 / 2 = 0\)
\(MS_E = SSE / df_E = 3.0 / 6 = 0.5\)
F值计算：
\(F_A = MS_A / MS_E = 16.335 / 0.5 = 32.67\)
\(F_B = MS_B / MS_E = 8.33 / 0.5 = 16.66\)
\(F_{AB} = MS_{AB} / MS_E = 0 / 0.5 = 0\)
方差分析表
| 变异来源 | 平方和SS | 自由度 | 均方MS | F值 |
|---|---|---|---|---|
| 预警策略（A） | 32.67 | 2 | 16.335 | 32.67 |
| 居住区域（B） | 8.33 | 1 | 8.33 | 16.67 |
| 交互作用（A×B） | 0 | 2 | 0 | 0 |
| 误差 | 3.0 | 6 | 0.5 | - |
| 总计 | 43.33 | 11 | - | - |
- 统计决策与结论
给定 \(\alpha=0.05\)，临界值：\(F_{0.05}(2,6)=5.14\)，\(F_{0.05}(1,6)=5.99\)。
交互作用：\(F_{AB}=0 < 5.14\)，不拒绝 \(H_0\)。结论：不存在显著交互作用，说明策略的相对效果在两类区域中是一致的（社区联动在两种区域都是最优，基础文本都是最差）。
预警策略：\(F_A=32.67 > 5.14\)，拒绝 \(H_0\)。结论：不同策略对储备天数的影响显著。
居住区域：\(F_B=16.67 > 5.99\)，拒绝 \(H_0\)。结论：老旧小区的储备天数显著高于新建小区。
library(ggplot2)
library(dplyr)
# --- 准备数据 ---
data_10_3 <- data.frame(
strategy = factor(rep(c("A1(基础文本)", "A2(多媒体增强)", "A3(社区联动)"), times = 4)),
region = factor(rep(c("B1(老旧小区)", "B2(新建小区)"), each = 6)),
days = c(4, 5, 6, 7, 8, 9, # A1, A2, A3 在 B1 的数据
2, 3, 5, 6, 6, 7) # A1, A2, A3 在 B2 的数据
)
# --- 交互作用图 ---
# 目的：直观观察线条走势，预判是否存在交互作用。误差棒可以展示数据的随机波动范围。
interaction_plot <- ggplot(data_10_3, aes(x = strategy, y = days, color = region, group = region)) +
stat_summary(fun = mean, geom = "point", size = 3) +
stat_summary(fun = mean, geom = "line", size = 1, linetype = "dashed") +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
labs(title = "交互作用图",
subtitle = "线条大致平行，猜测可能不存在显著的交互效应",
x = "预警策略 (A)",
y = "平均储备天数",
color = "居住区域 (B)") +
theme_minimal()
print(interaction_plot)
# --- 统计检验 ---
anova_result_10_3 <- aov(days ~ strategy * region, data = data_10_3)
{
cat("\n=====方差分析表=====\n")
print(summary(anova_result_10_3))
}##
## =====方差分析表=====
## Df Sum Sq Mean Sq F value Pr(>F)
## strategy 2 8.167 4.083 0.875 0.464
## region 1 8.333 8.333 1.786 0.230
## strategy:region 2 0.167 0.083 0.018 0.982
## Residuals 6 28.000 4.667# --- 事后检验 ---
# 只有当交互作用不显著时(P > 0.05)，我们重点关注主效应的事后检验。这里我们模拟查看策略 A 的主效应。
tukey_result_10_3 <- TukeyHSD(anova_result_10_3, which = "strategy")
cat("\n=====事后检验 (策略 A)=====\n")##
## =====事后检验 (策略 A)=====print(tukey_result_10_3)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = days ~ strategy * region, data = data_10_3)
##
## $strategy
## diff lwr upr p adj
## A2(多媒体增强)-A1(基础文本) 0.75 -3.936867 5.436867 0.8781242
## A3(社区联动)-A1(基础文本) 2.00 -2.686867 6.686867 0.4408545
## A3(社区联动)-A2(多媒体增强) 1.25 -3.436867 5.936867 0.70635484. 非平衡数据的拓展
在前面的两节中，我们讨论的都是“平衡数据”，即每个实验组合（单元格）中的样本量 \(n\) 是相等的。平衡设计在统计分析中具有极大的优势，如计算公式简洁、各因素效应正交（互不干扰）、检验功效最高。然而，在实际的社会调查或科学实验中，由于种种不可控因素，我们往往无法获得完美的平衡数据。例如，计划在每个小区调查 10 人，但老旧小区配合度高，收回了 12 份问卷，而新建小区配合度低，仅收回了 5 份问卷；或者实验过程中出现样本损坏、动物死亡、问卷无效等情况导致某些组数据缺失。这种各单元格样本量 \(n_{ij}\) 不完全相同（有的格子2个数据，有的5个数据）的数据，我们称为非平衡数据。
对于非平衡数据的分析，其基本原理虽然仍然是进行变异分解，但在具体计算和统计性质上会变得更为复杂。首先，由于样本量在各个单元格内不等，简单的算术平均值不再能准确代表边际效应，必须使用加权平均来计算主效应。更为重要的是，此时各因素的平方和不再独立，一个因素的计算结果会受到另一个因素数据分布的影响（即效应不正交）。为了解决这一问题，在分析非平衡数据时，需要特别注意平方和类型的选择。
通常推荐使用Type III（结果不受变量顺序影响，计算主效应时会剔除其他因素及交互作用的影响，对非平衡数据最稳健），因为它能修正样本量不等带来的偏差，计算出每个因素在剔除其他因素影响后的独立效应，特别是在存在交互作用时，Type III 是最稳健的方法。相比之下，Type I（序贯型，按模型输入顺序依次“瓜分”变异，先进入模型的变量会抢占重叠部分的解释力）受变量输入顺序影响较大，一般不推荐；Type II（分层型，假设无交互作用，计算每个主效应时只剔除其他主效应，不剔除交互效应）仅适用于没有交互作用的情况。此外，如果数据中出现了某个格子完全没有数据（\(n_{ij}=0\)，即“空单元格”）的极端情况，标准的双因素方差分析模型将无法估计完整的交互作用，这时就需要采用特殊的线性模型方法进行处理，或者需要剔除部分水平才能进行有效分析。
R 的 aov() 默认是 Type I。要在 R 中做 Type III 分析，通常使用 car 包的 Anova() 函数，并指定参数 type = “III”，同时需要注意 R 默认的对比编码 contrasts 设置（通常改为 contr.sum）。
例10-4
非平衡数据示意：假设在例10-3的调研中，A3B1（社区联动在老旧小区）组引起了研究者浓厚兴趣，多收集了3个数据，而其他组保持不变。此时数据结构变为非平衡：
| B₁（老旧小区） | B₂（新建小区） | 样本量 | |
|---|---|---|---|
| A₁（基础文本） | 4, 5 | 2, 3 | \(n=2\) |
| A₂（多媒体增强） | 6, 7 | 5, 6 | \(n=2\) |
| A₃（社区联动） | 8, 9, 8, 7, 9 | 6, 7 | \(n_{31}=5, n_{32}=2\) |
10.2 卡方检验
10.2.1 非参数检验
同学们可能已经注意到，我们前面学习的t检验、Z检验以及F检验等假设检验方法，都有一个共同的前提条件：需要假定总体服从某种特定的分布（最常见的是正态分布）。这类依赖于总体分布假设的检验方法，我们称之为参数检验。与参数检验相对应，还存在另一大类重要的检验方法，称为非参数检验。
非参数检验最核心的特点是，它不依赖于总体分布的具体数学形式。也就是说，在进行检验时，我们不需要事先知道或不要求总体必须服从如正态分布之类的特定分布。它的检验统计量通常不是基于总体参数（如均值、方差）的估计推导而出，而是来源于对数据本身一些更基础、更直观特征的转换和利用，例如将数值转换为秩次、只关注其符号、或直接统计类别出现的频数。这种处理方式绕开了对总体分布形态的具体依赖。其统计推断的基石，也随之转变为该统计量在原假设条件下、依据随机性原理所导出的抽样分布。
因此，非参数检验更适用于各种分布类型的数据，特别是当数据严重偏离正态分布（如存在极端异常值、分布高度偏态）时，其结果比参数检验更稳健。它不仅能处理连续数据，还能处理参数检验难以处理的顺序数据（如满意度等级排名）和分类数据（如性别、产品类型）等。
但由于放弃了数据的部分数值信息（例如，将具体数值转换为排序秩次），可能会导致检验效能降低。也就是说，如果数据实际上满足参数检验的条件，那么非参数检验需要更大的样本量才能检测出同样大小的真实差异。因此，如果已知或有充分理由确信总体分布，且数据满足参数检验的条件，我们就优先使用参数检验。反之，当总体分布不明、数据严重偏离参数检验前提或处理的是等级、分类数据时，我们更偏向使用非参数检验。
10.2.2 卡方检验的基本思想
卡方检验是一种应用于分类数据的非参数检验方法，其核心在于评估观察频数与期望频数之间的差异是否具有统计显著性。所谓期望频数，是指在特定的原假设（如变量间无关联、分布无差异等）成立时，理论上我们预期会观察到的频数。卡方检验的基本思想是：如果实际观察到的频数与期望频数之间的差异很小，可能仅由随机抽样误差导致，则支持原假设；反之，如果差异过大，超出了随机误差可以解释的合理范围，则有理由怀疑原假设的正确性。
卡方检验的统计量定义为： \[\chi^2 = \sum \frac{(O_i - E_i)^2}{E_i}\] 其中：
\(O_i\) 是第 i 个类别或单元格的观察频数。
\(E_i\) 是第 i 个类别或单元格的期望频数。
该统计量在原假设成立且样本量足够大时，近似服从卡方分布。其统计原理源于以下逻辑：对于每个单元格，其标准化残差 \((O_i - E_i)/\sqrt{E_i}\) 在原假设下近似服从标准正态分布。而k个独立的标准正态随机变量的平方和服从自由度为k的卡方分布（回想一下前面说的卡方分布和正态分布的关系）。因此，卡方统计量本质上是这些标准化残差的平方和。不过，由于列联表中的频数受到行合计与列合计等条件的约束，变量并非完全独立，因此其自由度df并非简单的单元格数目，而需要根据具体的检验类型和数据维度（如\(df = (r-1)\times(c-1)\)）来确定。
遵循假设检验的通用框架，卡方检验的实施包含以下标准步骤：
建立假设（H₀ 通常为“无差异”、“独立”或“符合某种分布”）。
计算检验统计量 χ² 的值。
确定显著性水平 α 并计算相应的自由度 df。
查卡方分布临界值表或计算 p 值。
比较统计量与临界值（或比较 p 值与 α），做出拒绝或不拒绝原假设的统计决策。
10.2.3 三种卡方检验
a) 拟合优度检验
目的：检验一个分类变量的观察频数分布与某种理论分布（如正态分布、均匀分布、特定比例分布）是否吻合。
假设：
H₀: 观察频数的分布与理论分布无显著差异（拟合良好）
H₁: 观察频数的分布与理论分布有显著差异（拟合不好）
自由度计算：\(df = k - 1 - m\)
\(k\)：分类的类别数
\(m\)：需要从数据中估计的理论分布的参数个数（例如，若检验是否符合正态分布，且需要从数据估计均值和标准差，则 m=2）
例子：一位学者怀疑一枚骰子可能不均匀。为了检验其均匀性，他将该骰子投掷了600次，并记录下各面朝上的次数。
观察频数 (Oᵢ):
| 骰子面 | 1 | 2 | 3 | 4 | 5 | 6 | 合计 |
|---|---|---|---|---|---|---|---|
| 频数 | 95 | 105 | 110 | 90 | 100 | 100 | 600 |
检验步骤：
- 建立假设
零假设 (H₀)：该骰子是均匀的，即每个面出现的概率均为 1/6。（观察分布与理论均匀分布无显著差异）
备择假设 (H₁)：该骰子是不均匀的，即至少有一个面出现的概率不等于 1/6。（观察分布与理论均匀分布存在显著差异）
- 计算检验统计量 χ²
计算期望频数 (Eᵢ)：在H₀成立的条件下，每个面的期望频数 = 总次数 × 理论概率。
\[ E_i = 600 \times \frac{1}{6} = 100 \quad \text{（对于每一面 i=1,2,...,6）} \]
计算卡方统计量： \[ \begin{align*} \chi^2 &= \sum_{i=1}^{6} \frac{(O_i - E_i)^2}{E_i} \\ &= \frac{(95-100)^2}{100} + \frac{(105-100)^2}{100} + \frac{(110-100)^2}{100} + \frac{(90-100)^2}{100} + \frac{(100-100)^2}{100} + \frac{(100-100)^2}{100} \\ &= \frac{25}{100} + \frac{25}{100} + \frac{100}{100} + \frac{100}{100} + \frac{0}{100} + \frac{0}{100} \\ &= 0.25 + 0.25 + 1 + 1 + 0 + 0 \\ &= 2.5 \end{align*} \]
- 确定显著性水平 α 和自由度 df
设定显著性水平 \(\alpha = 0.05\)。
计算自由度 \(df\)：此为拟合优度检验，类别数 $k=6，未从数据中估计参数 \(m=0\)。
\[ df = k - 1 - m = 6 - 1 - 0 = 5 \]
- 查卡方分布表或计算 p 值
临界值：查表得，在 \(df=5, \alpha=0.05\)水平下，卡方分布的临界值 \(\chi^2_{0.05}(5) \approx 11.07\)。
p值法：根据计算出的 \(\chi^2 = 2.5\) 和 \(df=5\)，可求得 p 值远大于 0.05。
- 做出统计决策
决策：由于检验统计量 \(\chi^2 = 2.5 < 11.07\)（或者说，p 值 > 0.05），我们没有足够的证据拒绝零假设 (H₀)。
结论：在 0.05 的显著性水平下，本次检验结果不支持“该骰子不均匀”的说法。也就是说，从这600次投掷的数据来看，该骰子可以被认为是均匀的。
# 加载必要的包
library(ggplot2)
library(dplyr)
# 拟合优度检验
observed <- data.frame(
face = factor(1:6),
count = c(95, 105, 110, 90, 100, 100)
)
# 执行卡方检验
# p = rep(1/6, 6) 是默认设置，
chi_result <- chisq.test(observed$count, p = rep(1/6, 6))
{
cat("\n=====拟合优度检验=====\n")
print(chi_result)
}##
## =====拟合优度检验=====
##
## Chi-squared test for given probabilities
##
## data: observed$count
## X-squared = 2.5, df = 5, p-value = 0.7765# 绘图：观测值 vs 期望值
# 使用 ggplot2 直接绘制，直观对比
ggplot(observed, aes(x = face, y = count, fill = face)) +
geom_col(alpha = 0.7, color = "white") +
geom_hline(yintercept = 100, linetype = "dashed", color = "red", linewidth = 1) +
labs(title = "卡方检验：骰子的公平性",
subtitle = paste("p-value:", format(chi_result$p.value, digits = 3)),
y = "Frequency",
x = "Dice Face") +
theme_minimal() 
b) 独立性检验（列联分析）
目的：检验两个分类变量之间是否相互独立（即是否存在关联）。
数据形式：列联表（Contingency Table），其基本形式如下
一个 \(r \times c\) 的列联表（\(r\) 行，\(c\) 列）：
| 列1 | 列2 | … | 列c | 行合计 | |
|---|---|---|---|---|---|
| 行1 | O₁₁ / E₁₁ | O₁₂ / E₁₂ | … | O₁c / E₁c | R₁ |
| 行2 | O₂₁ / E₂₁ | O₂₂ / E₂₂ | … | O₂c / E₂c | R₂ |
| … | … | … | … | … | … |
| 行r | Or₁ / Er₁ | Or₂ / Er₂ | … | Orc / Erc | Rr |
| 列合计 | C₁ | C₂ | … | Cc | N (总样本数) |
- \(O_{ij}\)：第 \(i\) 行、第 \(j\) 列单元格的观察频数。
- \(E_{ij}\)：第 \(i\) 行、第 \(j\) 列单元格的期望频数（在变量独立的假设下计算）。
- \(R_i\)：第 \(i\) 行的行合计。
- \(C_j\)：第 \(j\) 列的列合计。
- \(N\)：总样本数。
假设： - H₀: 两个变量相互独立 - H₁: 两个变量不独立（即存在显著关联）
自由度计算：\(df = (r - 1) \times (c - 1)\) - \(r\)：列联表的行数 - \(c\)：列联表的列数
期望频数计算：对于列联表中第 i 行第 j 列的单元格，其期望频数为： \[E_{ij} = \frac{(第i行合计) \times (第j列合计)}{总样本数} = \frac{R_i \times C_j}{N}\]
例子：一位心理学家想研究咖啡饮用量（高/低）与工作压力水平（高/中/低）之间是否存在关联。她对150名办公室职员进行了调查，得到如下数据。
观察频数 (Oᵢⱼ) 表:
| 高压 | 中压 | 低压 | 行合计 (Rᵢ) | |
|---|---|---|---|---|
| 高咖啡 | 35 | 25 | 10 | 70 |
| 低咖啡 | 15 | 30 | 35 | 80 |
| 列合计 (Cⱼ) | 50 | 55 | 45 | 150 (N) |
检验步骤：
- 建立假设
零假设 (H₀)：咖啡饮用量与工作压力水平相互独立。
备择假设 (H₁)：咖啡饮用量与工作压力水平不独立。
- 计算检验统计量 χ²
步骤 2.1：计算每个单元格的期望频数 (Eᵢⱼ)， 使用公式 \(E_{ij} = \frac{R_i \times C_j}{N}\)
- \(E_{11}\) (高咖啡 & 高压): \((70 \times 50) / 150 = 23.33\)
- \(E_{12}\) (高咖啡 & 中压): \((70 \times 55) / 150 = 25.67\)
- \(E_{13}\) (高咖啡 & 低压): \((70 \times 45) / 150 = 21.00\)
- \(E_{21}\) (低咖啡 & 高压): \((80 \times 50) / 150 = 26.67\)
- \(E_{22}\) (低咖啡 & 中压): \((80 \times 55) / 150 = 29.33\)
- \(E_{23}\) (低咖啡 & 低压): \((80 \times 45) / 150 = 24.00\)
将期望频数填入表格，便于对比： 期望频数 (Eᵢⱼ) 表: | | 高压 | 中压 | 低压 | | :— | :—: | :—: | :—: | | 高咖啡 | 23.33 | 25.67 | 21.00 | | 低咖啡 | 26.67 | 29.33 | 24.00 |
步骤 2.2：计算每个单元格的 \((O_{ij} - E_{ij})^2 / E_{ij}\)
- 高咖啡 & 高压: \((35-23.33)^2 / 23.33 \approx 5.84\)
- 高咖啡 & 中压: \((25-25.67)^2 / 25.67 \approx 0.02\)
- 高咖啡 & 低压: \((10-21.00)^2 / 21.00 \approx 5.76\)
- 低咖啡 & 高压: \((15-26.67)^2 / 26.67 \approx 5.11\)
- 低咖啡 & 中压: \((30-29.33)^2 / 29.33 \approx 0.02\)
- 低咖啡 & 低压: \((35-24.00)^2 / 24.00 \approx 5.04\)
步骤 2.3：求和得到卡方统计量 \[ \chi^2 = 5.84 + 0.02 + 5.76 + 5.11 + 0.02 + 5.04 \approx 21.79 \]
- 确定显著性水平 α 和自由度 df
设定 \(\alpha = 0.05\)。
行数 \(r=2\)，列数 \(c=3\)，自由度 \(df = (2-1) \times (3-1) = 2\)。
- 查卡方分布表或计算 p 值
在 \(df=2, \alpha=0.05\) 水平下，临界值 \(\chi^2_{0.05}(2) \approx 5.99\)。
计算出的 \(\chi^2 = 21.79 > 5.99\)，且 p 值远小于 0.001。
- 做出统计决策
决策：由于 \(\chi^2 > \chi^2_{critical}\)（且 p 值 < 0.05），我们拒绝零假设 (H₀)。
结论：在 0.05 的显著性水平下，有充分的统计证据表明咖啡饮用量与工作压力水平之间存在显著关联。观察数据可见，高咖啡饮用者中高压力的比例更高，而低咖啡饮用者中低压力的比例更高。
library(ggplot2)
library(dplyr)
# 独立性检验（咖啡与工作压力）
# 1. 创建数据框
coffee_df <- data.frame(
coffee = rep(c("HCoffee", "LCoffee"), each = 3),
stress = rep(c("High Stress", "Medium Stress", "Low Stress"), 2),
frequency = c(35, 25, 10, 15, 30, 35)
)
# 2. 转换为列联表并执行卡方检验
coffee_table <- xtabs(frequency ~ coffee + stress, data = coffee_df)
chi_result <- chisq.test(coffee_table)
{
cat("=====独立性检验=====\n")
print(chi_result)
}## =====独立性检验=====
##
## Pearson's Chi-squared test
##
## data: coffee_table
## X-squared = 21.774, df = 2, p-value = 1.87e-05# 计算百分比用于绘图
coffee_percent <- coffee_df %>%
group_by(coffee) %>%
mutate(percentage = frequency / sum(frequency) * 100)
# 绘图：堆叠百分比条形图
ggplot(coffee_percent, aes(x = coffee, y = percentage, fill = stress)) +
geom_col(alpha = 0.8, width = 0.6) +
geom_text(aes(label = paste0(round(percentage, 1), "%")),
position = position_stack(vjust = 0.5),
color = "white", size = 3.5) +
scale_fill_brewer(palette = "Set1") +
labs(title = "咖啡消费与压力程度的卡方检验",
subtitle = paste("p-value:", format(chi_result$p.value, digits = 3)),
y = "Percentage within Group (%)",
x = "Coffee Consumption",
fill = "Stress Level") +
theme_minimal()
c) 同质性检验
目的：检验两个或多个总体的某个分类变量的分布是否相同（即是否同质）。同质性检验的数据形式和计算方法与独立性检验完全相同，都是使用列联表和相同的卡方统计量公式。两者的区别在于研究设计和抽样的出发点：
独立性检验：从一个总体中抽取一个样本，然后观测每个个体在两个变量上的表现。
同质性检验：先从多个总体中分别独立抽样，然后观测每个个体在一个变量上的表现。
假设：
H₀: 多个总体的分布是相同的（同质的）
H₁: 至少有一个总体的分布与其他总体不同
例子：一家制药公司研发了一种新疫苗。为了检验疫苗副作用（无/轻微/严重）的分布在不同年龄组（青年/中年/老年）中是否相同，公司分别从三个年龄组的总体中独立招募了各200名参与者，共计600人进行临床试验。数据记录如下：
观察频数 (Oᵢⱼ) 表:
| 无副作用 | 轻微副作用 | 严重副作用 | 行合计 (Rᵢ) | |
|---|---|---|---|---|
| 青年组 | 120 | 50 | 30 | 200 |
| 中年组 | 100 | 60 | 40 | 200 |
| 老年组 | 80 | 70 | 50 | 200 |
| 列合计 (Cⱼ) | 300 | 180 | 120 | 600 (N) |
检验步骤：
- 建立假设
- 零假设 (H₀)：三个年龄组的副作用分布是同质的（即分布相同）。
- 备择假设 (H₁)：至少有一个年龄组的副作用分布与其他组不同。
- 计算检验统计量 χ²
- 步骤 2.1：计算每个单元格的期望频数 (Eᵢⱼ) 使用公式 \(E_{ij} = \frac{R_i \times C_j}{N}\)
- \(E_{11}\) (青年组 & 无): \((200 \times 300) / 600 = 100.00\)
- \(E_{12}\) (青年组 & 轻微): \((200 \times 180) / 600 = 60.00\)
- \(E_{13}\) (青年组 & 严重): \((200 \times 120) / 600 = 40.00\)
- \(E_{21}\) (中年组 & 无): \((200 \times 300) / 600 = 100.00\)
- \(E_{22}\) (中年组 & 轻微): \((200 \times 180) / 600 = 60.00\)
- \(E_{23}\) (中年组 & 严重): \((200 \times 120) / 600 = 40.00\)
- \(E_{31}\) (老年组 & 无): \((200 \times 300) / 600 = 100.00\)
- \(E_{32}\) (老年组 & 轻微): \((200 \times 180) / 600 = 60.00\)
- \(E_{33}\) (老年组 & 严重): \((200 \times 120) / 600 = 40.00\)
- 步骤 2.2：计算每个单元格的 \((O_{ij} -
E_{ij})^2 / E_{ij}\)
- 青年组 & 无: \((120-100)^2 / 100 = 4.00\)
- 青年组 & 轻微: \((50-60)^2 / 60 \approx 1.67\)
- 青年组 & 严重: \((30-40)^2 / 40 = 2.50\)
- 中年组 & 无: \((100-100)^2 / 100 = 0.00\)
- 中年组 & 轻微: \((60-60)^2 / 60 = 0.00\)
- 中年组 & 严重: \((40-40)^2 / 40 = 0.00\)
- 老年组 & 无: \((80-100)^2 / 100 = 4.00\)
- 老年组 & 轻微: \((70-60)^2 / 60 \approx 1.67\)
- 老年组 & 严重: \((50-40)^2 / 40 = 2.50\)
- 步骤 2.3：求和得到卡方统计量 \[ \chi^2 = 4.00 + 1.67 + 2.50 + 0.00 + 0.00 + 0.00 + 4.00 + 1.67 + 2.50 \approx 16.34 \]
- 确定显著性水平 α 和自由度 df
- 设定 \(\alpha = 0.05\)。
- 行数 \(r=3\)（三个年龄组），列数 \(c=3\)（三种副作用），自由度 \(df = (3-1) \times (3-1) = 4\)。
- 查卡方分布表或计算 p 值
- 在 \(df=4, \alpha=0.05\) 水平下，临界值 \(\chi^2_{0.05}(4) \approx 9.49\)。
- 计算出的 \(\chi^2 \approx 16.34 > 9.49\)，且 p 值远小于 0.05。
- 做出统计决策
- 决策：由于 \(\chi^2 > \chi^2_{critical}\)，我们拒绝零假设 (H₀)。
- 结论：在 0.05 的显著性水平下，有充分的统计证据表明疫苗副作用的分布在青年、中年、老年三个年龄组中存在显著差异，即分布不同质。观察数据可见，青年组报告“无副作用”的比例最高，而老年组报告“轻微”和“严重”副作用的比例相对更高。
library(ggplot2)
library(dplyr)
# 同质性检验（疫苗副作用年龄分布）
# 创建数据框
vaccine_df <- data.frame(
age_group = rep(c("Young", "Middle", "Elderly"), each = 3),
side_effect = rep(c("None", "Mild", "Severe"), 3),
frequency = c(120, 50, 30, 100, 60, 40, 80, 70, 50)
)
# 转换为列联表并执行卡方检验
vaccine_table <- xtabs(frequency ~ age_group + side_effect, data = vaccine_df)
chi_result <- chisq.test(vaccine_table)
# 输出结果 (自定义风格)
{
cat("=====同质性检验=====\n")
print(chi_result)
}## =====同质性检验=====
##
## Pearson's Chi-squared test
##
## data: vaccine_table
## X-squared = 16.333, df = 4, p-value = 0.002603# 计算百分比用于绘图
vaccine_percent <- vaccine_df %>%
group_by(age_group) %>%
mutate(percentage = frequency / sum(frequency) * 100)
# 绘图：堆叠百分比条形图
ggplot(vaccine_percent, aes(x = age_group, y = percentage, fill = side_effect)) +
geom_col(alpha = 0.8, width = 0.6) +
geom_text(aes(label = paste0(round(percentage, 1), "%")),
position = position_stack(vjust = 0.5),
color = "white", size = 3.5) +
scale_fill_brewer(palette = "Set1") +
labs(title = "不同年龄组的副作用卡方检验",
subtitle = paste("p-value:", format(chi_result$p.value, digits = 3)),
y = "Percentage within Group (%)",
x = "Age Group",
fill = "Side Effect") +
theme_minimal()
10.3 从单变量统计到多变量分析
本章讨论的方差分析与卡方检验，标志着统计学从单变量推断迈向多变量分析的关键转折。尽管这两种方法所处理的对象均明确涉及两个或多个变量（如方差分析中的因素与结果，卡方检验中的列联表变量），但在深入其原理时，同学们或许已察觉到一种“化繁为简”的逻辑特质：表面上探讨的是变量间的二元关系，但在核心操作层面上，实质上仍回归于“单变量”的统计推断。
若同学们产生了这种直觉，实则已触及了经典多变量分析的本质：所谓的“多变量分析”，在多数情境下，遵循的是“多变量的分析设计，单变量的统计推断”这一基本范式。
无论是方差分析还是卡方检验，其方法论的核心均在于将一部分变量 \(X\)（如实验分组、分类特征）设定为“控制条件”，从而将统计推断的焦点完全聚焦于另一个响应变量 \(Y\)（如测量数值、频次计数）之上。此处的统计推断始终严格遵循单变量的概率法则：
- 在方差分析中：我们将样本按特定属性分组（即给定 \(X\)），所考察的核心统计量，乃是不同组别下因变量 \(Y\) 的数值分布差异。
- 在卡方检验中：我们在既定的列联表结构下（即给定 \(X\) 与 \(Y\) 的交叉框架），所关注的焦点，乃是单元格内“频数”这一计数指标的随机波动与理论期望的偏离程度。
由此可见，真正的统计对象始终仅有一个随机变量（\(Y\) 或频数）。因此，前几章所构建的单变量理论基础（如抽样分布、中心极限定理及假设检验逻辑），不仅是充分够用的，更是构建上述多变量方法的坚实基石。
进一步地，这种思维将成为贯穿后续高级统计课程的通用逻辑主线。从基础的线性回归，到进阶的逻辑回归，乃至广义线性模型等复杂体系，尽管其外在形式涉及错综复杂的变量网络，但剖析其内核，无一不践行着“固定条件、推断单变量”的核心思想。理解这一点，我们才能够从单变量向多变量跨越的过程中，建立起统一且连贯的知识体系。