常见统计检验的本质都是线性模型(或:如何教统计学)上篇

本文转载自统计之都,文章翻译自   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)

3 Pearson相关性和Spearman相关性

3.1 理论:作为线性模型

模型:yy 的形式是一个斜率(β1)乘以 xx 加上一个截距(β0),也就是一条直线。
y=β0+β1x     H0:β1=0
以上模型,实际上是我们熟悉的旧公式 y=ax+b(这里书写顺序变为 y=b+ax) 的一个更数学化的表达。在 R 软件里面,我们比较懒,所以写成了 y ~ 1 + x,R 对此理解为 y=1*数 + x*另一个数,而且,t 检验,线性回归等等都只是去寻找能够最好地预测 y 的数字。
无论你怎么书写,截距(β0)和斜率(β1)都可以确定一条直线:

这是所谓的回归模型,当右边有多个 ββ 和自变量相乘的时候,它扩展为多元回归。以下的所有模型,从  单样本 t 检验 到  双因素方差分析 ,都是这个模型的特殊形式,不多不少。
顾名思义, Spearman 秩相关系数 就是  x 和  y 的秩变换后的  Pearson 相关系数
很快我就会介绍   这个概念。现在,线性模型的相关系数等价于 “真正的” Pearson 相关系数,但 P 值是近似值,这个近似值适用于  N > 10  的情况,并且在  N > 20  时  几乎完美
很多学生都没有意识到这么漂亮和神奇的等价关系!给线性回归带上数据标签来,我们可立刻看到秩转换过程:

3.2 理论:秩转换

对于一串数字,秩(rank)意思是使用它们的排序号来替换它们(第 1 最小的,第 2 最小的,第 3 最小的,以此类推)。因此 rank(c(3.6, 3.4, -5.0, 8.2)) 的秩转换结果是 3, 2, 1, 4。从上面的图形里看出来了吗?符号秩是一样的,先根据绝对值排序,再添加上数值前面的符号。所以上面的符号秩是 2, 1, -3, 4。用代码表示如下:
signed_rank = function(x) sign(x) * rank(abs(x))
我希望我说秩很容易理解的时候没有冒犯到其他人。这就是转换大部分 “非参” 数检验到它们的对应参数检验所要做的全部事情!一个重要的推论是很多 “非参” 检验和它们的对应参数检验版本都有一致的参数:均值、标准差、齐次方差等等 —— 区别在于它们是在秩转换后的数据上计算的。这是为什么我把 “非参” 用引号包起来。

3.2.1 R 代码:Pearson 相关系数

在 R 里运行这些模型再容易不过了。它们产生相同的 p 值和 t 值,但是这里有个问题:lm 返回斜率,尽管它通常比相关系数 r 更容易理解和反映了更多信息,但你依然想得到 r 值。幸运地,如果 x 和 y 有相同的标准差,斜率就会变成 r。所以在这里我们使用scale(x)使得SD(x)=1.0和SD(y)=1.0:

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
结果:

置信区间没有完全一致,但是非常相近。

3.2.2  R 代码: Spearman 相关系数

注意,我们可以把斜率解释为:对于每一 x 的秩的变化,所获得的相应  y 的秩的变化。我认为这个数字非常有趣。然而,截距更难解释,因为它定义在rank(x=0)的时候,然而这其实是不可能的,因为 x 是从 1 开始的。
查看相同的  r  (即这里的  rho ) 和  p

# Spearman correlation

a <- cor.test(y, x, method = "spearman") # Built-in

b <- lm(rank(y) ~ 1 + rank(x)) # Equivalent linear model

让我们看一下结果:

4 单均值

4.1 单样本t检验和Wilcoxon符号秩检验

4.1.1  理论:作为线性模型

t检验 模型:单独一个数字来预测y。
y = β0    H0:β0 = 0
换句话说,这是我们所熟悉的y = β0+β1*x,只是最后一项消失了,因为 x不存在了(等价于 x=0,见下方左图) 以上模型,一旦用y的 符号秩 来替换了y本身(见下方右图),就和  Wilcoxon 符号秩检验 非常相近。
signed_rank (y)=β0

这个近似对于样本量大于14的情况已足够好,对于大于50的情况 接近完美

4.1.2  R代码:单样本t检验

尝试运行以下 R 代码,确认线性模型(lm)和内置的 t.test 产生相同的 t、p、r。lm 的输出没有置信区间,但是你可以用 confint(lm(...)) 来确认结果也是相同的:

# Built-in t-test

a <- t.test(y)


# Equivalent linear model: intercept-only

b <- lm(y ~ 1)

4.1.3  R代码:Wilcoxon符号秩检验

除了一致的p值,lm也提供了符号秩均值,我发现这个数字非常有信息量。

# Built-in

a <- wilcox.test(y)


# Equivalent linear model

b <- lm(signed_rank(y) ~ 1) # See? Same model as above, just on signed ranks


# Bonus: of course also works for one-sample t-test

c <- t.test(signed_rank(y))

结果

4.2 配对样本t检验和Wilcoxon配对组检验

4.2.1  理论:作为线性模型

t检验 模型:一个数字(截距)来预测组间之差。

y2−y1=β0   H0:β0=0

这意味着只有一个y=y2-y1需要预测,而且它变成了对于组间之差的 单样本t检 验。 因此可视化效果和单样本 t 检验是相同的。冒着过度复杂化简单作差的风险,你可以认为这些组间之差是斜率(见图的左半部分),我们也可以用  y  的差来表达(见图的右半部分):

类似地,Wilcoxon配对组和Wilcoxon符号秩的唯一 差别,就是它对配对的差y-x的符号秩进行检验。

signed_rank(y2-y1)=β0   H0:β0 = 0

4.2.2  R代码:配对样本t检验

a <- t.test(y, y2, paired = TRUE) # Built-in paired t-test

b <- lm(y - y2 ~ 1) # Equivalent linear model

结果:

4.2.3  R代码:Wilcoxon配对组检验

我们再一次运用符号秩转换技巧。这依然是近似值,但是非常接近:

# Built-in Wilcoxon matched pairs

a <- wilcox.test(y, y2, paired = TRUE)


# Equivalent linear model:

b <- lm(signed_rank(y - y2) ~ 1)


# Bonus: identical to one-sample t-test ong signed ranks

c <- t.test(signed_rank(y - y2))


结果:
对于大样本量(N>>100),这计算方式某种程度上比较接近符号检验。但是本例子中这种近似效果较差。

5 双均值

5.1 独立t检验和Mann-Wnitney U检验

5.1.1  理论:作为线性模型

独立 t 检验模型:两个均值来预测  y。
y i =β0+β1xi    H0:β1=0

上式中,xi是示性变量(0 或 1),用于示意数据点 ii 是从一个组里采样还是另一个组里采样的。 示性变量 indicator variable 或哑变量 dummy variable 或者 one-hot 编码 存在于很多线性模型当中,我们很快就会看到它有什么用途。

Mann-Whitney U 检验 (也被称为两个独立组的  Wilcoxon 秩和检验 ;这次没有 符号 秩了)是有着非常接近的近似的相同模型,除了它不是在原有值而是在  x 和  y 的秩上计算的:
rank(yi)=β0+β1xi    H0:β1=0
对我来说,这种等价性使“非参”统计量更容易理解了。这种近似在每个组样本量大于 11 的时候比较合适,在每个组样本量大于 30 的时候看起来  相当完美

5.1.2  理论:示性变量

示性变量可以用图像来帮助理解。‍‍‍‍‍‍这个变量在 x 轴,所以第一个组的数据点位于 x=0,第二个组的位于x=1。然后 β0是截距(蓝线),β1是两个均值之间的斜率(红线) β 1  。为什么?因为当 Δx=1的时候,斜率等于相差值:

s lope =Δy/Δx=Δy/1=Δy=difference

奇迹啊!即使类别之间的差值也可以用线性模型来表达,这真的是一把瑞士军刀!

5.1.3  理论:示性变量(后续)

如果你觉得你理解了示性变量了,可以直接跳到下一章节。这里是对示性变量更为详细的解释:
如果数据点采样自第一个组,即 xi=0 ,模型就会变成  y= β0+β1⋅0=β0 。换句话 说,模型预测的值是 beta0。这些数据点的均值 β是最好的预测,而 β0是第 1 组的均值。

另一方面,采样自第二个组的数据点有 xi=1 ,所以模型变了 yi=β0+β1⋅1=β0+β1。 换句话说,我们加上了 β1 ,从第一组的均值移动到了第二组的均值。所以 β1 成为了两个组的 均值之差

举个例子,假设第 1 组人是 25 岁 (β0=25) ,第 2 组人 28 岁 (β1=3) ,那么对于第 1 组的人的模 型是 y=25+3⋅0=25,第 2 组的人的模型 是 y=25+3⋅1=28。

呼,搞定!对于初学者,理解示性变量需要一些时间,但是只需要懂得加法和乘法就能上手了!

5.1.4  R代码:独立t检验

提醒一下,当我们在R里写y~1+x,它是y=β0⋅1+β1⋅x的简写,R会为你计算β值。因此y~1+x是R里面表达y=a⋅x+b的形式。

注意相等的t、df、p和估计值。我们可以用confint(lm(...))获得置信区间。

# 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"))
结果:

5.1.5  R代码:Mann-Whitney U  检验

# Wilcoxon / Mann-Whitney U

a <- wilcox.test(y, y2)


# As linear model with our dummy-coded group_y2:

b <- lm(rank(value) ~ 1 + group_y2) # compare to linear model above

5.2  Welth t检验

这等价于以上(学生) 独立t检验 ,除了学生t检验假设同方差,而Welch t检验没有这个假设。所以线性模型是相同的,区别在于,我们对每一个组指定一个方差。我们可用nlme包:

# Built-in

a <- t.test(y, y2, var.equal = FALSE)


# As linear model with per-group variances

b <- nlme::gls(value ~ 1 + group_y2, method = "ML",

weights = nlme::varIdent(form = ~ 1 | group))

结果:

6 三个或多个均值

方差分析ANOVA是只有类别型自变量的线性模型,它们可以简单地扩展上述模型,并重度依赖示性变量。如果你还没准备好,一定要去阅读 示性变量一节

6.1 单因素方差分析和Kruskal-Wallis检验

6.1.1  理论:作为线性模型

模型:每组一个均值来预测 y。

其中‍‍ xi 是示性变量‍‍(x=0 或‍‍ x=1)‍‍,且最多只有一‍‍个 xi=1 ‍‍且其余  x i = 0 。‍

注意这和我们已做的其它模型“有很大的相同之处”。如果只有两个组,这个模型就‍‍是 y=β0+β1∗x,‍‍即 独立 t 检验。如果只有一个组,这就‍‍是 y=β0,‍‍即  单样本 t 检验 。从图中很容易看出来 —— 只要遮盖掉一些组然后看看图像是否对上了其它可视化结果。

单因素方差分析有一个对数线性版本,称为 拟合优度 检验,我们稍后会讲到。顺便一说,因为我们现在对多个 x 进行回归,因此单因素方差分析对应的是 多元回归 模型。
Kruskal-Wallis  检验是 y(value)秩转换的 单因素方差分析
rank(y)=β0+β1x1+β2x2+β3x3+…
这个近似在数据点达到 12 或更多的时候已经 足够好 。同样地,对一个或两个组做这个检验,也已经有对应等式,分别为 Wilcoxon 符号秩检验 Mann-Whitney U 检验

6.1.2  示例数据

首先创建可能取值为 a、b、c 的类别变量,那么 单因素方差分析 基本上成为了 三样本 t 检验 ,然后手动对每个组创建 示性变量

# Three variables in "long" formatN <- 20 # Number of samples per groupD <- 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

伴随着组别 a 的截距全都展示了出来,我们看到每一行有且仅有另一个组 b 或组 c 的参数添加进去,用于预测 value。因此组 b 的数据点永远不会影响到组 c 的估计值。

6.1.3  R代码:单因素方差分析

好的,接下来我们看看一个专用的 方差分析 函数( car::Anova )和手动创建示性变量的  lm  线性模型结果是否一致:

# 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

结果:

实际上,car::Anova 和 aov 就是 lm 包装而来的,所以得到相同的结果一点也不令人感到意外。线性模型的示性变量公式是缩写语法 y ~ factor 背后的模型。实际上,使用 aov 和 car::Anova 而不使用 lm 是为了得到一个漂亮的格式化的方差分析表。
lm 的默认输出包含了参数估计的结果(额外收获!),可以将上述 R 代码展开来看。因为它就是方差分析模型,所以你也可以用 coefficients(aov(...)) 得到参数估计的结果。
注意,我没有使用 aov 函数,是因为它计算了第一类平方和,这种计算方式不提倡。围绕着使用第二类平方和(car::Anova 默认)还是第三类平方和(使用 car::Anova(..., type=3))有着很多争论,我们这里略过不提。

6.1.4  R代码:Kruskal-Wallis检验

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))

结果:

6.2 双因素方差分析

6.2.1  理论:作为线性模型

模型:每组一个均值(主效应)加上这些均值乘以各个因子(交互效应)。

在一个更大的模型框架里,主效应实际上就是上述  单因素方差分析 模型。虽然交互效应只是一些数字乘以另外一些识数字,但是它更难解释。我会把这些解释内容留给课堂上的老师们,而聚焦在等价表达之上。:-)

使用矩阵记号:

y = β 0 + β 1 X 1 + β 2 X 2 + β 3 X 1 X 2       H 0 : β 3 = 0

这里 βi是 β 的分量,其中之后一个会被示性向量 Xi X i 所选取。这里出现的 H 0 是交互效应。注意,截距是 β0是所有因子中第一个水平的均值,而其它所有的  β 都是相对于β0的值。

继续探究上文单因素方差分析的数据集,我们加上交互因子  mood (情绪),这样就能够测试  group:mood  的交互效应( 3 × 2  的 ANOVA)。同样,为了使用线性模型  lm ,将这个因子转为 示性变量

# Crossing factor

D$mood <- c("happy", "sad")

# Dummy coding

D$mood_happy <- ifelse(D$mood == "happy", 1, 0)

# 1 if mood==happy. 0 otherwise.

# D$mood_sad = ifelse(D$mood == 'sad', 1, 0) # 同上,但是我们不需要同时设置

β0 成为了 a 组的开心者!

6.2.2  R代码:双因素方差分析

现在转向 R 里的真实建模。对比专用的 ANOVA 函数(car::Anova;原因参见 单因素方差分析 一节),以及线性模型函数(lm)。注意,在单因素方差分析里,一次性检验全部因子的交互效应,这其中涉及到多个参数(这个例子里有两个参数)。所以我们不能查看总体模型估计结果或者任意一个参数结果。因此,使用 似然比检验(likelihood-ratio test) 来比较带交互效应的双因素方差分析模型和没有交互效应的模型。anova 函数就是用来做这个检验的。这看起来在作弊,实际上,它只是在已经拟合了的模型上计算似然(likelihood)、p 值等等!

# 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

结果:
下面展示了近似的主因素模型。方差分析主效应的精确计算需要 更加繁琐的计算步骤,并且依赖于使用第二型还是第三型平方和来进行推断。
在模型的汇总统计量里可以找到对比以上  Anova  拟合的主因素效应的数值。

# Main effect of group.

e <- lm(value ~ 1 + group_b + group_c, D)


# Main effect of mood.

f <- lm(value ~ 1 + mood_happy, D)

6.2 协方差分析

协方差分析 ANCOVA 是在方差分析的基础上添加连续型回归变量,从而模型里包含了连续型以及(示性变量转化来的)类别型变量。沿用上文 单因素方差分析 的例子,加上 age,就变成 单因素协方差分析 (one-way ANCOVA):

y=β0+β1x1+β2x2+…+β3age

其中, x i  是熟悉的示性变量。 β0 β 0 是第一个组在 age=0时候的均值。 可以用这个方法把所有的方差分析转化成为协方差分析,比如说,在上一节的 双因素方差分析 加上 βN⋅age。 我们先继续探究单因素协方差分析,从添加 age到模型里开始:

# Update data with a continuous covariate

D$age <- D$value + rnorm_fixed(nrow(D), sd = 3) # Correlated to value

比起使用位于 x 轴的位置,最好使用颜色来区分各组数据。β依然是数据的平均yy移动量,唯一的不同是现在用斜率而不是截距来对每一组进行建模。换句话说,单因素方差分析实际上是对于每一组(y=β0)的 单样本 t 检验 ,而 单因素协方差分析 实际上是每一组(yi=β0+βi+β1⋅age)的 Pearson 相关性检验

这里是使用线性模型来做单因素方差分析的 R 代码:

#Dedicated ANCOVA functions. The order of factors matter in pure-aov (type-I variance).

#Use type-II (default for car::Anova) or type-III (set type=3),

a <- car::Anova(aov(value ~ group + age, D))

# a = aov(value ~ group + age, D) # Predictor order matters. Not nice!


#As dummy-coded linear model.

full <- lm(value ~ 1 + group_b + group_c + age, D)


# Testing main effect of age using Likelihood-ratio test

null_age <- lm(value ~ 1 + group_b + group_c, D) # Full without age. One-way ANOVA!

result_age <- anova(null_age, full)


# Testing main effect of groupusing Likelihood-ratio test

null_group <- lm(value ~ 1 + age, D) # Full without group. Pearson correlation!

result_group <- anova(null_group, full)

结果:

(未完,见次条)
译者简介
黄俊文 ,本科中山大学,硕士 University of California, San Diego。现于业界从事数据分析工作。
作者: Jonas Kristoffer Lindeløv
审校: 蔡占锐,黄湘云,谢益辉
编辑: 孙腾飞,任焱

统计之都官网: cosx.org

统计之都公众号: CapStat

我来评几句
登录后评论

已发表评论数()

相关站点

热门文章