本文转载自统计之都，文章翻译自 Jonas Kristoffer Lindeløv 的 Common statistical tests are linear models (or: how to teach stats)，翻译工作已获得原作授权。本翻译工作首发于统计之都网站和微信公众号上。
y = ax + b 然而很不幸的是，统计入门课程通常把各种检验分开教学，给学生和老师们增加了很多不必要的麻烦。在学习每一个检验的基本假设时，如果不是从线性模型切入，而是每个检验都死记硬背，这种复杂性又会成倍增加。因此，我认为先教线性模型，然后对线性模型的一些特殊形式进行改名是一种优秀的教学策略，这有助于更深刻地理解假设检验。线性模型在频率学派、贝叶斯学派和基于置换的U检验的统计推断之间是相通的，对初学者而言，从模型开始比从 P 值、第 I 类错误、贝叶斯因子或其它地方更为友好。# 加载必要的 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)
signed_rank = function(x) sign(x) * rank(abs(x))在 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
置信区间没有完全一致，但是非常相近。
r （即这里的 rho） 和 p：# Spearman correlationa <- cor.test(y, x, method = 'spearman') # Built-inb <- lm(rank(y) ~ 1 + rank(x)) # Equivalent linear model
这个近似对于样本量大于14的情况已足够好，对于大于50的情况接近完美。
# Built-in t-testa <- t.test(y)
# Equivalent linear model: intercept-onlyb <- lm(y ~ 1)# 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))
y2−y1=β0 H0:β0=0
类似地，Wilcoxon配对组和Wilcoxon符号秩的唯一 差别，就是它对配对的差y-x的符号秩进行检验。
signed_rank（y2-y1）=β0 H0：β0 = 0
a <- t.test(y, y2, paired = TRUE) # Built-in paired t-testb <- lm(y - y2 ~ 1) # Equivalent linear model# 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))
上式中，xi是示性变量（0 或 1），用于示意数据点 ii 是从一个组里采样还是另一个组里采样的。 示性变量 indicator variable 或哑变量 dummy variable 或者 one-hot 编码 存在于很多线性模型当中，我们很快就会看到它有什么用途。
slope=Δy/Δx=Δy/1=Δy=difference
奇迹啊！即使类别之间的差值也可以用线性模型来表达，这真的是一把瑞士军刀！
另一方面，采样自第二个组的数据点有 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。
呼，搞定！对于初学者，理解示性变量需要一些时间，但是只需要懂得加法和乘法就能上手了！
提醒一下，当我们在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'))
# 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
# 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))其中 xi 是示性变量（x=0 或 x=1），且最多只有一个 xi=1 且其余
注意这和我们已做的其它模型“有很大的相同之处”。如果只有两个组，这个模型就是 y=β0+β1∗x，即 独立 t 检验。如果只有一个组，这就是 y=β0，即 单样本 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 codegroup_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 的估计值。
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 modela <- 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))
在一个更大的模型框架里，主效应实际上就是上述 单因素方差分析模型。虽然交互效应只是一些数字乘以另外一些识数字，但是它更难解释。我会把这些解释内容留给课堂上的老师们，而聚焦在等价表达之上。:-)
使用矩阵记号：
这里 βi是 β 的分量，其中之后一个会被示性向量 Xi 所选取。这里出现的
继续探究上文单因素方差分析的数据集，我们加上交互因子 mood（情绪），这样就能够测试 group:mood 的交互效应（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) # 同上，但是我们不需要同时设置现在转向 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)y=β0+β1x1+β2x2+…+β3age
# 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 相关性检验：
#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)
结果：
联系客服
