如何解释 t 检验以外的统计检验中的成对观察（例如回归）？
How to account for paired observations in statistical tests other than t-test (e.g. regression)?
如何解释 t-test 以外的统计检验中的成对观察？下面我将讨论两个示例，其中我尝试使用混合效果方法来实现但失败了。
示例1：如何在lm()中重现t.test(..., paired=T)?
# simulate data
set.seed(66)
x1 <- rnorm(n=100, mean=30, sd=6)
x2 <- rnorm(n=100, mean=60, sd=6)
# arrange the data in a dataset
dd <- data.frame(ID=rep(paste("ID", seq(1,100, by=1), sep="_"),2),
response=c(x1,x2),
group=c(rep("A", 100), rep("B", 100))
)
t.test(x1,x2, paired=F)
summary(lm(response~group, data=dd)) # same outcome
如果观察结果是成对的，可以用 t.test() 来解释，但如何在 lm() 中解释（如果可能的话）？我尝试使用混合效应模型方法，但是：
summary(lmerTest::lmer(response~group + (1+group|ID), data=dd))
报错：
Error: number of observations (=200) <= number of random effects (=200) for term (1 + group | ID);
the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
同时：
summary(lmerTest::lmer(response~group + (1|ID), data=dd))
运行但固定效应参数估计和相关标准。误差和 t 值与 lm().
产生的相同
示例 2：具有配对观察值的线性回归
让我们假设我创建的数据集中的观察结果来自间隔 30 天测量的受试者——也就是说，100 名受试者中的每一个都在第 0 天测量，然后在第 30 天再次测量——我们想估计随时间变化：
dd$time=c(rep(0,100), rep(30, 100)) # add "time" variable to dd
数据如下所示（黑色线性回归，红线链接成对数据）：
lm1 <- lm(response~time, data=dd)
lm1 没有说明观察的成对性质。
我想到了 运行 一个混合效应模型，允许每对数据在截距和斜率上有所不同，但 R 再次抗议我试图估计太多参数：
lmerTest::lmer(response ~ time + (time | ID), data=dd)
# Error: number of observations (=200) <= number of random effects (=200) for term (time | ID);
# the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
一个更简单的模型，允许数据对的截距不同但斜率不同，即：
lmer(response ~ time + (1 | ID), data=dd)
抱怨：
boundary (singular) fit: see ?isSingular
但运行并产生与 lm().
产生的相同的固定效应估计
[更新]
@Limey 提醒我，配对 t 检验只不过是评估两组之间的成对差异是否不为零的 t 检验。这种成对差异可用于执行 t 检验之外的任何配对统计检验。为了验证这一点，我创建了三个不同的“响应”变量，它们是 x1 和 x2 以不同方式排序的组合（分别是：原始随机顺序；x1 递增和 x2 降序；两者均升序）。
dd$response2 <- c(sort(x1, decreasing = FALSE), sort(x2, decreasing = T))
dd$response3 <- c(sort(x1, decreasing = FALSE), sort(x2, decreasing = F))
我计算了对应的差值：
dd$diff1 <- c((dd$response[1:100]-dd$response[1:100]),
(dd$response[101:200]-dd$response[1:100]))
dd$diff2 <- c((dd$response2[1:100]-dd$response2[1:100]),
(dd$response2[101:200]-dd$response2[1:100]))
dd$diff3 <- c((dd$response3[1:100]-dd$response3[1:100]),
(dd$response3[101:200]-dd$response3[1:100]))
并用它们来执行线性模型：
lm2.diff1 <- lm(diff1~time, data=dd)
lm2.diff2 <- lm(diff2 ~time, data=dd)
lm2.diff3 <- lm(diff3 ~time, data=dd)
我预计他们的斜率估计会有所不同，但他们都是一样的：
summary(lm2.diff1)$coeff[2] # 0.9993754
summary(lm2.diff2)$coeff[2] # 0.9993754
summary(lm2.diff3)$coeff[2] # 0.9993754
它们的斜率估计值与相应的“未配对”线性模型（lm(response~time)、lm(response2~time) 和 lm(response3~time)）的估计值相同。我错过了什么？
配对 t 检验只是测试（两组之间差异的）平均值是否为零。因此，要通过其他方式“模拟”配对 t 检验的结果，只需预先计算差异并将其传递给您选择的检验。使用您的示例：
x1 <- rnorm(n=100, mean=30, sd=6)
x2 <- rnorm(n=100, mean=60, sd=6)
diff <- x1-x2
dd <- data.frame(response=diff)
# Standard paired t-test
t.test(x1,x2, paired=T)
Paired t-test
data: x1 and x2
t = -36.167, df = 99, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-31.93760 -28.61546
sample estimates:
mean of the differences
-30.27653
现在是“模拟”t 检验...
t.test(diff)
One Sample t-test
data: diff
t = -36.167, df = 99, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-31.93760 -28.61546
sample estimates:
mean of x
-30.27653
现在作为线性模型
summary(lm(response~1, data=dd))
Call:
lm(formula = response ~ 1, data = dd)
Residuals:
Min 1Q Median 3Q Max
-18.473 -7.328 0.614 6.101 20.764
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -30.2765 0.8371 -36.17 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 8.371 on 99 degrees of freedom
Their slope estimate is the same estimated from the corresponding
"unpaired" linear models (lm(response~time), lm(response2~time),
and lm(response3~time)). What am I missing?
这三个模型的整体斜率相同是有道理的，是所有成对斜率的平均值（这很容易确认）。
我缺少的是，在 lm2.diff3 的情况下，斜率估计周围的 StdErr 大约比 lm2.diff1 和 lm2.diff2 小一个数量级。
lm2.diff3 是在 Response3 上执行的，与 Response1 和 Response2 相比，每对观测值的行为更加统一，因此其斜率估计值周围的 StdErr 更小。
好问题！这里有一些棘手的问题。
- 我们可以通过
t.test() 和手动计算配对测试开始：
pairedtest1 <- t.test(x1,x2, paired=TRUE)
d <- x1-x2
n <- length(x1)
tstat <- mean(d)/(sd(d)/sqrt(n)) ## -37.58846
pval <- 2*pt(abs(tstat), lower.tail=FALSE, df=n-1) ## 2.065802e-60
all.equal(pairedtest1$p.value,pval) ## TRUE
all.equal(unname(pairedtest1$statistic),tstat) ## TRUE
- 正如您所发现的，在线性混合模型中使用随机截距 和 组间治疗差异的模型会过度拟合，因为每个模型只有一个观察值每组治疗；给出了一个类似的例子 here。如果你愿意，你可以强制
lmer 来适应它：
m0 <- lme4::lmer(response~group+(group|ID), data=dd, REML=TRUE,
control=lmerControl(check.nobs.vs.nRE="ignore", calc.deriv=FALSE))
（请注意，我们还可以通过比较对数似然或 REML 标准来查看两个模型是否给出等效拟合 - 当我们有像这样的过度拟合模型时，不同的模型可以通过模型的不同线性组合获得等效拟合参数...)
如果我们运行
library(lmerTest)
coef(summary(as(m1,"lmerModLmerTest"),ddf="Kenward-Roger"))["groupB",]
（这里默认的Satterthwaite ddf计算失败）我们得到
Estimate Std. Error df t value Pr(>|t|)
2.998126e+01 7.976192e-01 9.900000e+01 3.758844e+01 2.065922e-60
t 统计量和 p 值与上面的结果非常接近（我本可以说 summary() 但想提取系数的这一特定行 table）。
- 现在让我们试试简化模型
m1 <- lme4::lmer(response~group+(1|ID), data=dd, REML=TRUE)
如您所述，拟合是奇异的（如果您检查，RE 方差将列为 0）。这里的 t 统计量和 p 值有点偏差（目前，我不太确定为什么以前的模型有效！）。原因是，对于此数据集，组内方差略 大于 组间方差。一般来说，我们期望 var(between) = sigma^2_between + sigma_2_within/n，这可以达到 asymptotically/in 期望，但在小数据集中，排序可以按照我们在这里观察到的方向进行，在这种情况下，我们需要 negative方差以完美拟合数据。
resids <- with(dd,response-ave(response,group, FUN=mean))
wv <- var(resids-ave(resids, dd$ID, FUN=mean)) ## 15.82
bv <- var(tapply(resids, list(dd$ID), FUN=mean)) ## 14.92
（如果我们用 lme 拟合相同的模型，它 出现 OK 并给我们一个小的 [但非零] 组间截距方差估计，但是如果我们尝试 intervals(m2) 它会抱怨近似的 var-cov 矩阵不是正定的...）
- 正如 Gang Chen 在 2011 post to r-sig-mixed-models 中指出的那样，我们可以通过允许组内两个观察值之间的 负 相关性来拟合我们想要的模型（加法上面显示的模型只允许非负相关）：
library(nlme)
m2 <- lme(response~group,random=list(ID=pdCompSymm(form=~group-1)),
data=dd,method="REML")
all.equal(tstat^2, anova(m2)[["F-value"]][2]) ## TRUE
all.equal(pval, anova(m2)[["p-value"]][2]) ## TRUE
anova() 的 p 值与我们上面的结果匹配，F 统计量与我们的 t 统计量的 平方 匹配。
- 我们也可以在
glmmTMB 中这样做：我们必须要小心一点，因为 glmmTMB 中默认的 cs() 协方差结构允许每个项有不同的方差，而 corCompSymm 施加齐次方差：
m3 <- glmmTMB::glmmTMB(response~group+cs(group-1|ID),
data=dd, REML=TRUE)
m4 <- update(m3, map=list(theta=factor(c(1,1,2))))
（map参数将随机效应参数向量的前两个元素设置为相同，对应于第一组和第二组中变异的log-sd）
系数 table 获得了正确的 t 统计量（它称为 z 值），但没有“分母自由度”的概念，Z 检验也是如此，而不是t 检验 ...
coef(summary(m4))$cond["groupB",]
Estimate Std. Error z value Pr(>|z|)
groupB 29.98126 0.7976217 37.58832 0
如何解释 t-test 以外的统计检验中的成对观察？下面我将讨论两个示例，其中我尝试使用混合效果方法来实现但失败了。
示例1：如何在lm()中重现t.test(..., paired=T)?
# simulate data
set.seed(66)
x1 <- rnorm(n=100, mean=30, sd=6)
x2 <- rnorm(n=100, mean=60, sd=6)
# arrange the data in a dataset
dd <- data.frame(ID=rep(paste("ID", seq(1,100, by=1), sep="_"),2),
response=c(x1,x2),
group=c(rep("A", 100), rep("B", 100))
)
t.test(x1,x2, paired=F)
summary(lm(response~group, data=dd)) # same outcome
如果观察结果是成对的，可以用 t.test() 来解释，但如何在 lm() 中解释（如果可能的话）？我尝试使用混合效应模型方法，但是：
summary(lmerTest::lmer(response~group + (1+group|ID), data=dd))
报错：
Error: number of observations (=200) <= number of random effects (=200) for term (1 + group | ID);
the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
同时：
summary(lmerTest::lmer(response~group + (1|ID), data=dd))
运行但固定效应参数估计和相关标准。误差和 t 值与 lm().
示例 2：具有配对观察值的线性回归
让我们假设我创建的数据集中的观察结果来自间隔 30 天测量的受试者——也就是说，100 名受试者中的每一个都在第 0 天测量，然后在第 30 天再次测量——我们想估计随时间变化：
dd$time=c(rep(0,100), rep(30, 100)) # add "time" variable to dd
数据如下所示（黑色线性回归，红线链接成对数据）：
lm1 <- lm(response~time, data=dd)
lm1 没有说明观察的成对性质。
我想到了 运行 一个混合效应模型，允许每对数据在截距和斜率上有所不同，但 R 再次抗议我试图估计太多参数：
lmerTest::lmer(response ~ time + (time | ID), data=dd)
# Error: number of observations (=200) <= number of random effects (=200) for term (time | ID);
# the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
一个更简单的模型，允许数据对的截距不同但斜率不同，即：
lmer(response ~ time + (1 | ID), data=dd)
抱怨：
boundary (singular) fit: see ?isSingular
但运行并产生与 lm().
[更新]
@Limey 提醒我，配对 t 检验只不过是评估两组之间的成对差异是否不为零的 t 检验。这种成对差异可用于执行 t 检验之外的任何配对统计检验。为了验证这一点，我创建了三个不同的“响应”变量，它们是 x1 和 x2 以不同方式排序的组合（分别是：原始随机顺序；x1 递增和 x2 降序；两者均升序）。
dd$response2 <- c(sort(x1, decreasing = FALSE), sort(x2, decreasing = T))
dd$response3 <- c(sort(x1, decreasing = FALSE), sort(x2, decreasing = F))
我计算了对应的差值：
dd$diff1 <- c((dd$response[1:100]-dd$response[1:100]),
(dd$response[101:200]-dd$response[1:100]))
dd$diff2 <- c((dd$response2[1:100]-dd$response2[1:100]),
(dd$response2[101:200]-dd$response2[1:100]))
dd$diff3 <- c((dd$response3[1:100]-dd$response3[1:100]),
(dd$response3[101:200]-dd$response3[1:100]))
lm2.diff1 <- lm(diff1~time, data=dd)
lm2.diff2 <- lm(diff2 ~time, data=dd)
lm2.diff3 <- lm(diff3 ~time, data=dd)
我预计他们的斜率估计会有所不同，但他们都是一样的：
summary(lm2.diff1)$coeff[2] # 0.9993754
summary(lm2.diff2)$coeff[2] # 0.9993754
summary(lm2.diff3)$coeff[2] # 0.9993754
它们的斜率估计值与相应的“未配对”线性模型（lm(response~time)、lm(response2~time) 和 lm(response3~time)）的估计值相同。我错过了什么？
配对 t 检验只是测试（两组之间差异的）平均值是否为零。因此，要通过其他方式“模拟”配对 t 检验的结果，只需预先计算差异并将其传递给您选择的检验。使用您的示例：
x1 <- rnorm(n=100, mean=30, sd=6)
x2 <- rnorm(n=100, mean=60, sd=6)
diff <- x1-x2
dd <- data.frame(response=diff)
# Standard paired t-test
t.test(x1,x2, paired=T)
Paired t-test
data: x1 and x2
t = -36.167, df = 99, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-31.93760 -28.61546
sample estimates:
mean of the differences
-30.27653
现在是“模拟”t 检验...
t.test(diff)
One Sample t-test
data: diff
t = -36.167, df = 99, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-31.93760 -28.61546
sample estimates:
mean of x
-30.27653
现在作为线性模型
summary(lm(response~1, data=dd))
Call:
lm(formula = response ~ 1, data = dd)
Residuals:
Min 1Q Median 3Q Max
-18.473 -7.328 0.614 6.101 20.764
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -30.2765 0.8371 -36.17 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 8.371 on 99 degrees of freedom
Their slope estimate is the same estimated from the corresponding "unpaired" linear models (
lm(response~time),lm(response2~time), andlm(response3~time)). What am I missing?
这三个模型的整体斜率相同是有道理的，是所有成对斜率的平均值（这很容易确认）。
我缺少的是，在 lm2.diff3 的情况下，斜率估计周围的 StdErr 大约比 lm2.diff1 和 lm2.diff2 小一个数量级。
lm2.diff3 是在 Response3 上执行的，与 Response1 和 Response2 相比，每对观测值的行为更加统一，因此其斜率估计值周围的 StdErr 更小。
好问题！这里有一些棘手的问题。
- 我们可以通过
t.test()和手动计算配对测试开始：
pairedtest1 <- t.test(x1,x2, paired=TRUE)
d <- x1-x2
n <- length(x1)
tstat <- mean(d)/(sd(d)/sqrt(n)) ## -37.58846
pval <- 2*pt(abs(tstat), lower.tail=FALSE, df=n-1) ## 2.065802e-60
all.equal(pairedtest1$p.value,pval) ## TRUE
all.equal(unname(pairedtest1$statistic),tstat) ## TRUE
- 正如您所发现的，在线性混合模型中使用随机截距 和 组间治疗差异的模型会过度拟合，因为每个模型只有一个观察值每组治疗；给出了一个类似的例子 here。如果你愿意，你可以强制
lmer来适应它：
m0 <- lme4::lmer(response~group+(group|ID), data=dd, REML=TRUE,
control=lmerControl(check.nobs.vs.nRE="ignore", calc.deriv=FALSE))
（请注意，我们还可以通过比较对数似然或 REML 标准来查看两个模型是否给出等效拟合 - 当我们有像这样的过度拟合模型时，不同的模型可以通过模型的不同线性组合获得等效拟合参数...)
如果我们运行
library(lmerTest)
coef(summary(as(m1,"lmerModLmerTest"),ddf="Kenward-Roger"))["groupB",]
（这里默认的Satterthwaite ddf计算失败）我们得到
Estimate Std. Error df t value Pr(>|t|)
2.998126e+01 7.976192e-01 9.900000e+01 3.758844e+01 2.065922e-60
t 统计量和 p 值与上面的结果非常接近（我本可以说 summary() 但想提取系数的这一特定行 table）。
- 现在让我们试试简化模型
m1 <- lme4::lmer(response~group+(1|ID), data=dd, REML=TRUE)
如您所述，拟合是奇异的（如果您检查，RE 方差将列为 0）。这里的 t 统计量和 p 值有点偏差（目前，我不太确定为什么以前的模型有效！）。原因是，对于此数据集，组内方差略 大于 组间方差。一般来说，我们期望 var(between) = sigma^2_between + sigma_2_within/n，这可以达到 asymptotically/in 期望，但在小数据集中，排序可以按照我们在这里观察到的方向进行，在这种情况下，我们需要 negative方差以完美拟合数据。
resids <- with(dd,response-ave(response,group, FUN=mean))
wv <- var(resids-ave(resids, dd$ID, FUN=mean)) ## 15.82
bv <- var(tapply(resids, list(dd$ID), FUN=mean)) ## 14.92
（如果我们用 lme 拟合相同的模型，它 出现 OK 并给我们一个小的 [但非零] 组间截距方差估计，但是如果我们尝试 intervals(m2) 它会抱怨近似的 var-cov 矩阵不是正定的...）
- 正如 Gang Chen 在 2011 post to r-sig-mixed-models 中指出的那样，我们可以通过允许组内两个观察值之间的 负 相关性来拟合我们想要的模型（加法上面显示的模型只允许非负相关）：
library(nlme)
m2 <- lme(response~group,random=list(ID=pdCompSymm(form=~group-1)),
data=dd,method="REML")
all.equal(tstat^2, anova(m2)[["F-value"]][2]) ## TRUE
all.equal(pval, anova(m2)[["p-value"]][2]) ## TRUE
anova() 的 p 值与我们上面的结果匹配，F 统计量与我们的 t 统计量的 平方 匹配。
- 我们也可以在
glmmTMB中这样做：我们必须要小心一点，因为glmmTMB中默认的cs()协方差结构允许每个项有不同的方差，而corCompSymm施加齐次方差：
m3 <- glmmTMB::glmmTMB(response~group+cs(group-1|ID),
data=dd, REML=TRUE)
m4 <- update(m3, map=list(theta=factor(c(1,1,2))))
（map参数将随机效应参数向量的前两个元素设置为相同，对应于第一组和第二组中变异的log-sd）
系数 table 获得了正确的 t 统计量（它称为 z 值），但没有“分母自由度”的概念，Z 检验也是如此，而不是t 检验 ...
coef(summary(m4))$cond["groupB",]
Estimate Std. Error z value Pr(>|z|)
groupB 29.98126 0.7976217 37.58832 0