如何解决如何解释除 t 检验例如回归以外的统计检验中的配对观察？
如何解释除 t-test 之外的统计检验中的成对观察？下面我将讨论两个示例，其中我尝试使用混合效果方法进行操作但失败了。
示例 1：如何在 t.test(...,paired=T) 中重现 lm()？
# 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,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 递增和 {{1} } 降序；两者都是升序）。
x2
我计算了相应的差异：
dd$response2 <- c(sort(x1,decreasing = FALSE),sort(x2,decreasing = T))
dd$response3 <- c(sort(x1,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)）的估计值相同。我错过了什么？
解决方法
配对 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,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
好问题！这里有一些棘手的地方。
- 我们可以从通过
t.test()和手动计算配对测试开始：
pairedtest1 <- t.test(x1,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()，但想取出系数表的这一特定行）。
- 现在让我们尝试简化模型
m1 <- lme4::lmer(response~group+(1|ID),REML=TRUE)
正如您所指出的，拟合是奇异的（如果您选中，则 RE 方差将列为 0）。这里的 t 统计量和 p 值有点偏离（此时，我不太确定为什么以前的模型有效！）。原因是，对于这个数据集，组内方差比组间方差略大。一般来说，我们期望 var(between) = sigma^2_between + sigma_2_within/n，这在渐近/期望中起作用，但在小数据集中，排序可能与我们在这里观察到的方向一致，在这种情况下，我们需要一个负方差以完美拟合数据。
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 拟合相同的模型，它出现可以，并为我们提供一个小的 [但非零] 组间截距方差估计，但如果我们尝试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)),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 统计量的 square 相匹配。
- 我们也可以在
glmmTMB中这样做：我们必须小心一点，因为cs()中的默认glmmTMB协方差结构允许每个术语有不同的方差，而 {{1 }} 施加同质方差：
corCompSymm（m3 <- glmmTMB::glmmTMB(response~group+cs(group-1|ID),REML=TRUE)
m4 <- update(m3,map=list(theta=factor(c(1,1,2))))
参数将随机效应参数向量的前两个元素设置为相同）
系数表使 t 统计量（它称为 az 值）正确，但没有“分母自由度”的概念，因此 Z 检验而不是 t 检验......
map他们的斜率估计与相应的估计相同
“未配对”线性模型（lm(response~time)、lm(response2~time)、
和lm(response3~time)）。我错过了什么？
三个模型的总斜率相同是有道理的，即所有成对斜率的平均值（这很容易确认）。
我缺少的是，在 lm2.diff3 的情况下，斜率估计周围的 StdErr 大约比 lm2.diff1 和 lm2.diff2 中的小一个数量级。
lm2.diff3 在 Response3 上执行，与 Response1 和 Response2 相比，每对观察的行为更加一致，因此其斜率估计周围的 StdErr 更小。
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。