分享

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

 树悲风 2019-10-14
本文转载自统计之都,文章翻译自  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(Nmu = 0, sd = 1)
  scale(rnorm(N)) * sd + mu

# 图形风格
theme_axis <- function(Pjitter = 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)
  }
}

一开始,我们简单点,使用三组正态分布数据,且整理为宽(abc)和长(valuegroup)格式:

# Wide format (sort of)
# y = rnorm_fixed(50, mu=0.3, sd=2)  # Almost zero mean.
y <- c(rnorm(15), exp(rnorm(15)), runif(20min = -3max = 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 秩相关系数就是  和  的秩变换后的 Pearson 相关系数
很快我就会介绍  这个概念。现在,线性模型的相关系数等价于 “真正的” Pearson 相关系数,但 P 值是近似值,这个近似值适用于  的情况,并且在  时 几乎完美

很多学生都没有意识到这么漂亮和神奇的等价关系!给线性回归带上数据标签来,我们可立刻看到秩转换过程:


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 的秩的变化,所获得的相应  的秩的变化。我认为这个数字非常有趣。然而,截距更难解释,因为它定义在rank(x=0)的时候,然而这其实是不可能的,因为 x 是从 1 开始的。

查看相同的 r (即这里的 rho) 和 p
# Spearman correlationa <- cor.test(y, x, method = 'spearman') # Built-inb <- 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-testa <- t.test(y)
# Equivalent linear model: intercept-onlyb <- lm(y ~ 1)
结果:


4.1.3 R代码:Wilcoxon符号秩检验

除了一致的p值,lm也提供了符号秩均值,我发现这个数字非常有信息量。
# Built-ina <- wilcox.test(y)
# Equivalent linear modelb <- lm(signed_rank(y) ~ 1) # See? Same model as above, just on signed ranks
# Bonus: of course also works for one-sample t-testc <- t.test(signed_rank(y))
结果


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

4.2.1 理论:作为线性模型

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

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


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

类似地,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-testb <- lm(y - y2 ~ 1) # Equivalent linear model
结果:

4.2.3 R代码:Wilcoxon配对组检验

我们再一次运用符号秩转换技巧。这依然是近似值,但是非常接近:
# Built-in Wilcoxon matched pairsa <- 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 ranksc <- t.test(signed_rank(y - y2))
结果:
对于大样本量(N>>100),这计算方式某种程度上比较接近符号检验。但是本例子中这种近似效果较差。

5 双均值

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

5.1.1 理论:作为线性模型


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

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


Mann-Whitney U 检验(也被称为两个独立组的 Wilcoxon 秩和检验;这次没有符号秩了)是有着非常接近的近似的相同模型,除了它不是在原有值而是在  和  的秩上计算的:
rank(yi)=β0+β1xi    H0:β1=0

对我来说,这种等价性使“非参”统计量更容易理解了。这种近似在每个组样本量大于 11 的时候比较合适,在每个组样本量大于 30 的时候看起来 相当完美


5.1.2 理论:示性变量

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

slope=Δ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'10# 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 Ua <- 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-ina <- t.test(y, y2, var.equal = FALSE)
# As linear model with per-group variancesb <- 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 ‍‍且其余 xi=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 modela <- car::Anova(aov(value ~ group, D)) # Dedicated ANOVA functionb <- 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-inb <- 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+β1X1+β2X2+β3X1X2      H0:β3=0

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

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

# Crossing factorD$mood <- c('happy', 'sad')# Dummy codingD$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 functionsa <- car::Anova(aov(value ~ mood * group, D), type = 'II') # Normal notation. '*' both multiplies and adds main effectsb <- 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 modelnull <- lm(value ~ 1 + group_b + group_c + mood_happy, D) # Without interactionc <- 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

其中,xi 是熟悉的示性变量。β0β0是第一个组在 age=0时候的均值。可以用这个方法把所有的方差分析转化成为协方差分析,比如说,在上一节的双因素方差分析加上 βN⋅age。我们先继续探究单因素协方差分析,从添加 age到模型里开始:
# Update data with a continuous covariateD$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 testnull_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 testnull_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

    本站是提供个人知识管理的网络存储空间,所有内容均由用户发布,不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息,谨防诈骗。如发现有害或侵权内容,请点击一键举报。
    转藏 全屏 打印 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多