翻译：常见统计检验的本质是线性模型（或：如何教统计学）

本文翻译自  Jonas Kristoffer Lindeløv  的  Common statistical tests are linear models (or: how to teach stats)，翻译工作已获得原作授权。本翻译工作首发于统计之都网站和微信公众号上。

本文将常见的参数和 “非参” 数检验统一用线性模型来表示，在同一个框架下， 我们可以看到不同检验之间的许多相似之处，极富思考性和启发性。

1 常见检验的本质

大部分常见的统计模型（t 检验、相关性检验、方差分析（ANOVA）、卡方检验等） 是线性模型的特殊情况或者是非常好的近似。这种优雅的简洁性意味着我们学习起来不需要掌握太多的技巧。具体来说，这都来源于大部分学生从高中就学习的模型：y = ax + b 然而很不幸的是，统计入门课程通常把各种检验分开教学，给学生和老师们增加了很多不必要的麻烦。在学习每一个检验的基本假设时，如果不是从线性模型切入，而是每个检验都死记硬背，这种复杂性又会成倍增加。因此，我认为先教线性模型，然后对线性模型的一些特殊形式进行改名是一种优秀的教学策略，这有助于更深刻地理解假设检验。线性模型在频率学派、贝叶斯学派和基于置换的U检验的统计推断之间是相通的，对初学者而言，从模型开始比从 P 值、第 I 类错误、贝叶斯因子或其它地方更为友好。

在入门课程教授“非参”数检验的时候，可以避开 骗小孩 的手段，直接告诉学生“非参”检验其实就是和秩相关的参数检验。对学生来说，接受秩的概念比相信你可以神奇地放弃各种假设好的多。实际上，在统计软件 JASP 里，“非参”检验的贝叶斯等价模型就是使用 潜秩（Latent Rank）来实现的。频率学派的“非参”检验在样本量 N > 15 的时非常准确。

在来源和教材两章节，有很多类似（尽管更为散乱）的材料。我希望你们可以一起来提供优化建议，或者直接在 Github 提交修改。让我们一起来使本文章变得更棒！

2 设置和示例数据

如果你想查看函数和本笔记的其它设置的话，可以查看这段代码：

# 加载必要的 R 包用于处理数据和绘图
library(tidyverse)
library(patchwork)
library(broom)

# 设置随机数种子复现本文结果
set.seed(40)

# 生成已知参数的服从正态分布的随机数
rnorm_fixed <- function(N, mu = 0, sd = 1)
  scale(rnorm(N)) * sd + mu

# 图形风格
theme_axis <- function(P, jitter = FALSE,
                       xlim = c(-0.5, 2),
                       ylim = c(-0.5, 2),
                       legend.position = NULL) {
  P <- P + theme_bw(15) +
    geom_segment(
      x = -1000, xend = 1000,
      y = 0, yend = 0,
      lty = 2, color = "dark gray", lwd = 0.5
    ) +
    geom_segment(
      x = 0, xend = 0,
      y = -1000, yend = 1000,
      lty = 2, color = "dark gray", lwd = 0.5
    ) +
    coord_cartesian(xlim = xlim, ylim = ylim) +
    theme(
      axis.title = element_blank(),
      axis.text = element_blank(),
      axis.ticks = element_blank(),
      panel.border = element_blank(),
      panel.grid = element_blank(),
      legend.position = legend.position
    )

  # 是否设置抖动
  if (jitter) {
    P + geom_jitter(width = 0.1, size = 2)
  }
  else {
    P + geom_point(size = 2)
  }
}


一开始，我们简单点，使用三组正态分布数据，且整理为宽（a、b、c）和长（value、group）格式：

# Wide format (sort of)
# y = rnorm_fixed(50, mu=0.3, sd=2)  # Almost zero mean.
y <- c(rnorm(15), exp(rnorm(15)), runif(20, min = -3, max = 0)) # Almost zero mean, not normal
x <- rnorm_fixed(50, mu = 0, sd = 1) # Used in correlation where this is on x-axis
y2 <- rnorm_fixed(50, mu = 0.5, sd = 1.5) # Used in two means

# Long format data with indicator
value <- c(y, y2)
group <- rep(c("y1", "y2"), each = 50)

a <- cor.test(y, x, method = "pearson") # Built-in
b <- lm(y ~ 1 + x) # Equivalent linear model: y = Beta0*1 + Beta1*x
c <- lm(scale(y) ~ 1 + scale(x)) # On scaled vars to recover r

结果：

# Built-in independent t-test on wide data
a <- t.test(y, y2, var.equal = TRUE)

# Be explicit about the underlying linear model by hand-dummy-coding:
group_y2 <- ifelse(group == "y2", 1, 0) # 1 if group == y2, 0 otherwise
b <- lm(value ~ 1 + group_y2) # Using our hand-made dummy regressor

# Note: We could also do the dummy-coding in the model
# specification itself. Same result.
c <- lm(value ~ 1 + I(group == "y2"))

结果：

结果：

# Three variables in "long" format
N <- 20 # Number of samples per group
D <- data.frame(
  value = c(rnorm_fixed(N, 0), rnorm_fixed(N, 1), rnorm_fixed(N, 0.5)),
  group = rep(c("a", "b", "c"), each = N),

  # Explicitly add indicator/dummy variables
  # Could also be done using model.matrix(~D$group)
  # group_a = rep(c(1, 0, 0), each=N),  # This is the intercept. No need to code
  group_b = rep(c(0, 1, 0), each = N),
  group_c = rep(c(0, 0, 1), each = N)
) # N of each level

# Compare built-in and linear model
a <- car::Anova(aov(value ~ group, D)) # Dedicated ANOVA function
b <- lm(value ~ 1 + group_b + group_c, data = D) # As in-your-face linear model

结果：

a <- kruskal.test(value ~ group, D) # Built-in
b <- lm(rank(value) ~ 1 + group_b + group_c, D) # As linear model
# The same model, using a dedicated ANOVA function. It just wraps lm.
c <- car::Anova(aov(rank(value) ~ group, D)) 

# Dedicated two-way ANOVA functions
a <- car::Anova(aov(value ~ mood * group, D), type = "II") # Normal notation. "*" both multiplies and adds main effects
b <- car::Anova(aov(value ~ mood + group + mood:group, D)) # Identical but more verbose about main effects and interaction

# Testing the interaction terms as linear model.
full <- lm(value ~ 1 + group_b + group_c + mood_happy + group_b:mood_happy + group_c:mood_happy, D) # Full model
null <- lm(value ~ 1 + group_b + group_c + mood_happy, D) # Without interaction
c <- anova(null, full) # whoop whoop, same F, p, and Dfs

结果：

# Built-in test
a <- chisq.test(D$counts)

# As log-linear model, comparing to an intercept-only model
full <- glm(counts ~ 1 + mood_happy + mood_sad, data = D, family = poisson())
null <- glm(counts ~ 1, data = D, family = poisson())
b <- anova(null, full, test = "Rao")

# Note: glm can also do the dummy coding for you:
c <- glm(counts ~ mood, data = D, family = poisson())

# Contingency data in long format for linear model
D <- data.frame(
  mood = c("happy", "happy", "meh", "meh", "sad", "sad"),
  sex = c("male", "female", "male", "female", "male", "female"),
  Freq = c(100, 70, 30, 32, 110, 120)
)

# ... and as table for chisq.test
D_table <- D %>%
  tidyr::spread(key = mood, value = Freq) %>% # Mood to columns
  select(-sex) %>% # Remove sex column
  as.matrix()

# Dummy coding of D for linear model (skipping mood=="sad" and gender=="female")
# We could also use model.matrix(D$Freq~D$mood*D$sex)
D$mood_happy <- ifelse(D$mood == "happy", 1, 0)
D$mood_meh <- ifelse(D$mood == "meh", 1, 0)
D$sex_male <- ifelse(D$sex == "male", 1, 0)

# Built-in chi-square. It requires matrix format.
a <- chisq.test(D_table)

# Using glm to do a log-linear model, we get identical results when testing the interaction term:
full <- glm(Freq ~ 1 + mood_happy + mood_meh + sex_male + 
              mood_happy * sex_male + mood_meh * sex_male, data = D, family = poisson())
null <- glm(Freq ~ 1 + mood_happy + mood_meh + sex_male, data = D, family = poisson())
b <- anova(null, full, test = "Rao") # Could also use test='LRT' or test='Chisq'

# Note: let glm do the dummy coding for you
full <- glm(Freq ~ mood * sex, family = poisson(), data = D)
c <- anova(full, test = "Rao")

# Note: even simpler syntax using MASS:loglm ("log-linear model")
d <- MASS::loglm(Freq ~ mood + sex, D)

8 资料来源和更多的等价性模型

下面是本文内容的部分资料来源，还包含了很多本文没有提到的等价性模型：

• Cross Validated 网站上，我的原始想法。

• 对于“非参”检验，我之前提出的疑问 和有用的答案。

• StackOverflow 网站上，关于 t 检验和方差分析的问题和回答。

• Christoph Scheepers 的幻灯片，介绍了卡方检验如何被理解为对数线性模型。

• Philip M. Alday 的笔记，里面包括了卡方、二项、多项、泊松分布作为对数线性模型和 logistic 模型的理解。文中介绍的“等价性”没有我在本文展示的那么精确，因此我没有在本文详细介绍。然而，它们对理解这些检验是有帮助的！

• Kristoffer Magnusson  的文章使用 lme4::lmer 混合模型（ mixed model ）介绍了  RM-ANOVA  和增长模型（ growth model ）。

• Thom Baguley 的文章介绍了 Friedman 检验。这篇文章实际上启发了我开始思考“非参”检验的线性模型等价形式，而且最终推动我写下了本文章。

9 教材和课程大纲 

大部分高等统计书籍（和一些入门书籍）也都同意“所有模型都是 GLMM（广义线性混合效应模型） 的观点”。然而，线性模型部分通常都是概念上提了一下，而没有清晰地指出细节。我想通过简练的方式把线性模型当作工具。幸运地，大部分对初学者友好的教材后来都合并了：

• Russ Poldrack 的开源书籍 "Statistical Thinking for the 21st century"（从关于建模的第 5 章开始）

• Jeff Rouder 的课程笔记，介绍了仅使用 R^2 和 BIC 来对比模型。它避开了所有关于 p 值、F 值等等的繁琐问题。完整的材料和幻灯片可在这里找到。

我说一下对我所做的事情的看法。我已使用了本文的一部分进行教学，并获得了巨大的成功，但是这并不是完整的教学过程，因为我并没有分派到教授整个课程。

我会花费 50% 的时间在数据的线性模型上，因为它包含了学生所需知道的 70%（以下的第 1 点）。剩下来的课程则是关于当你有一个组、两个组等等数据的时候会发生什么事情。

注意，主流统计课程的开始部分都是关于采样和假设检验的理解；我这里把这部分移动到后面，这样，学生可以基于之前学习的知识来进行理解，而不是一上来就面对各种陌生的概念。

1. 回归的基础知识

1. 回想高中的知识：然后获得对斜率和截距的非常好的直觉。理解到这条式子能用所有的变量名称来重写：如 money = profit * time + starting_money或  或去除系数之后可写成 y ~ x + 1。如果听众接受程度高的话，可以探索这些模型是如何解微分方程的，并指出 y 是如何随着 x 的变化而变化的。

2. 扩展到多元回归模型。记得这时候要带有非常多的生活例子和练习，从而使这些概念变得直觉上非常容易理解。让听众感叹于这些简洁的模型都可以用来描述非常大的数据集。

3. 介绍对于非数值型数据如何进行秩转换，并进行各种尝试。

4. 教授三种前提假设：数据点的独立性，残差分布的正态性和方差齐性 （homoscedasticity）。

5. 参数的置信（confidence）/可信（credible）区间。指出极大似然估计（Maximum-Likelihood estimate）很难计算，因此区间估计更为重要。

6. 对以上简单的回归模型，简要地介绍 R^2。顺便提及一下，这就是 Pearson 和 Spearman 相关系数。

2. 特殊情况 #1：一个或两个均值（t 检验、Wilcoxon 检验、Mann-Whitney 检验）：

1. 单均值：当只有一个 x 值的时候，回归模型简化成了 y = b。如果 y 不是数值型的，你可以进行秩转换。应用模型假设（只有一个 x，因此方差齐性不适用于这里）。顺便提及一下，这些仅有截距的模型也分别可称为单样本 t 检验和 Wilcoxon 符号秩检验。

2. 双均值：如果我们把两个变量一起放在 x 轴，两者均值之差就是斜率。很好！这就能用我们称为瑞士军刀的线性模型来解决。应用模型的假设条件，检查两个组的方差是否相等，相等即方差齐性。这模型称为独立 t 检验。构造一些例子，做一些练习，也许还能加上 Welch 检验，再加上秩转换 ---- 变成所谓的 Mann-Whitney U 检验的版本。

3. 配对样本：违反了独立性假设。通过计算配对组的差值，这就转化成了 2.1（单截距）的等价形式，尽管这种情况有另外的名称：配对 t 检验和 Wilcoxon 配对组检验。

3. 特殊情况 #2：三个或多个均值（方差分析（ANOVA））

1. 对类别转化后的示性变量：类别的每一个取值范围对应的回归系数，是如何通过乘以一个二元（binary）示性变量，来对每个取值范围对应的截距来进行建模的。 （ How one regression coefficient for each level of a factor models an intercept for each level when multiplied by a binary indicator.） 这只是我们为了使数据能用线性模型建模，而扩展了在 2.1 所做的事情而已。

2. 一个变量的均值：单因素方差分析（one-way ANOVA）.

3. 两个变量的均值：双因素方差分析（two-way ANOVA）.

4. 特殊情况 #3：三个或多个比率（卡方检验）

1. 对数变换：通过对数变换，把“多元乘法”模型转化成线性模型，从而可以对比率进行建模。关于对数线性模型和对比率的卡方检验的等价性，可以查阅这个非常优秀的介绍。此外，还需要介绍 (log-) odds ratio（一般翻译为“比值比”或“优势比”）。当“多元乘法”模型使用对数变换转化为“加法”模型之后，我们仅加上来自 3.1 的示性变量技巧，就会在接下来发现模型等价于 3.2 和 3.3 的方差分析----除了系数的解释发生了变化。

2. 单变量的比率：拟合优度检验.

3. 双变量的比率：列联表.

5. 假设检验：

1. 视为模型比较的假设检验：假设检验用于全模型和某个参数固定了（通常为 0，也即从模型中去除）的模型进行比较，而不是对模型进行估计。比如说，在 t 检验 把两个均值之一固定为零之后，我们探究单独一个均值（单样本 t 检验）对两个组的数据的解释程度。如果解释程度比较好，那么我们更倾向于这个单均值模型，而不是双均值模型，因为前者更为简单。假设检验其实是比较多个线性模型，来获得更多的定量描述。单参数的检验，假设检验包含的信息更少。但是，同时对多个参数（如方差分析的类别变量）进行检验的话，模型比较就会变得没有价值了。

2. 似然比：似然比是一把瑞士军刀，它适用于单样本 t 检验到 GLMM 等情况。BIC 对模型复杂度进行惩罚。还有，加上先验（prior）的话，你就能得到贝叶斯因子（Bayes Factor）。一个工具，就能解决所有问题。我在上文方差分析中使用了似然比检验。

10 不足之处

一些需要澄清的简化前提：

1. 我没在这里覆盖到前提假设的内容。这会在另一篇文章揭晓！但是所有检验都很可能有三个预定假设：a)  数据点的独立性， b)  残差的正态性， c)  同方差性（homoscedasticity）。

2. 我假定所有的零假设是缺失了效应的情况，但是所有原理都和非 0 的零假设所一致的。

3. 我没有讨论推断内容。因为大家都会关心 p 值，因此我在比较中提到了 p 值，从而简短地展示了背后的模型等价性。参数的估计值也会展示出相同的等价性。如何进行推断则是另一个话题了。我个人是贝叶斯学派的，但是展示贝叶斯学派内容的话，会减少这篇文章的受众。此外，构造稳健模型是更好的选择，但是它无法揭示模型的等价性。

4. 本文列表依然缺失了很多其它著名的检验，有可能在以后添加进来。比如说符号检验（sign test）（要求很大的 N 从而可以有效地使用线性模型来近似），Friedman 检验 -- 即在 rank(y) 上的 RM-ANOVA，McNemar 检验，和二项（Binomial）/多项（Multinomial）检验。在链接一节可查阅更多的等价模型。如果你认为它们需要在本文提及到，欢迎在本文档的 Github 仓库 提交对应说明！

11 附录：翻译稿评论

11.1 译者评论

译者：相对于统计检验来说，线性回归实际上是更符合直觉的。想当年某检验实在让笔者百思不得其解，某师姐指点迷津：“你实在搞不懂可以看成是线性回归对系数的检验，我们如此这般建造一个 X ……”让笔者茅塞顿开。故听朋友推荐本文之后，笔者毛遂自荐承接了翻译任务。希望各位读者能从本文感受统计的威力和它带来的喜悦。如各位读者有指正或建议之处，热烈欢迎于主站或微信文下留言评论。

11.2 审稿人讨论

审稿人：我（黄湘云）看完这篇文章的感受是怀疑自己读了个“假”大学，开个玩笑哈！感觉这篇文章是继《心理学的危机》后又一篇需要找个地方安安静静读几天的，文中很多检验的细节都略过了，更加数学严格的检验介绍估计得去看《数理统计引论》陈希孺著的这本书才能明白。这篇文章的覆盖面起码是一个学期的课，如果把本文没有详述的其他检验补充进来，特别是加上检验的数学推导和一些实际应用案例后，估计能成一本书。我是学线性模型的（此线性模型非大多数人了解的彼线性模型），看完之后有点汗颜和如梦初醒，惊奇于作者独具一格的视角。不足之处是有些地方还不够通俗，比如列联表作为对数线性模型来理解，一点也不直接，作者也略去了！这里的列联表其实是指我们通常教科书上的独立性检验。拟合优度检验和独立性检验的检验统计量的极限分布都是卡方分布，故而都归纳在卡方检验下。

文章介绍了那么多的检验问题，实际上都可以归为统计学三大检验 --- 似然比检验、 Wald 检验、（Rao）Score（得分）检验 --- 在线性模型下的特殊情况。数理统计的教材往往是利用似然比这把瑞士军刀展开介绍的。似然比在假设检验中地位相当于极大似然估计在参数估计中的地位，相当于正态分布在抽样分布中的地位。抽样分布、参数估计和假设检验合称统计推断。学数理统计的人往往不愿去记那么多的检验名称，比如 t 检验、F 检验和卡方检验，特别是诸多名人检验，因为本质上那只是似然比统计量在不同的条件下呈现的极限分布不同而已。三大检验的渐近等价性可参考

• Dennis D. Boos and L. A. Stefanski.  Essential Statistical Inference: Theory and Methods. 2013. Springer, New York. <10.1007 style="font-size:inherit;color:inherit;line-height:inherit">

• “”   

在年给出得分检验的渐近性，文中提及得分检验这个名称只是强调在得分检验中的贡献，有些书籍中提及检验基本等同于得分检验，而拉格朗日乘子检验由和于年在经济学领域独立提出来的，所以在统计学文献中见最多的是得分检验，经济学文献中多描述为拉格朗日乘子检验。这二者是殊途同归，都对检验统计量的得分函数做泰勒展开，其极限分布都为卡方分布。

列联表是分类数据的一种方便的组织形式，与之相关的检验和前面带连续变量的线性回归模型的检验是有本质区别的，列联表是与多项分布联系起来的，这里没有残差，线性回归模型往往对残差做了独立同正态分布的假设。

译者简介

黄俊文，本科中山大学，硕士。现于业界从事数据分析工作。

作者：ø

译者：黄俊文

审校：蔡占锐，黄湘云，谢益辉

编辑：孙腾飞，任焱

统计之都：专业、人本、正直的中国统计学社区。

关注方式：扫描下图二维码。或查找公众号，搜索 统计之都 或即可。

往期推送：进入统计之都会话窗口，点击右上角小人图标，查看历史消息即可。

点击“阅读原文”，阅读原文