在t检验以外的统计检验(例如回归)中，如何解释配对观测？

Question

在t-test以外的统计检验中，如何解释成对观测？下面，我将讨论两个例子，其中我试图用混合效应的方法来做到这一点，但是失败了。

t.test(..., paired=T)示例1:如何在 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))

运行，但固定效果参数估计和相关的Std。错误和t值与lm()产生的错误和t值相同。

示例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))估计的。我遗漏了什么？

Answer 1

问得好！这里有几个棘手的地方。

1. 我们可以从通过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

1. 正如你所发现的，在线性混合模型中，当每组治疗只有一个观察值时，具有随机截距和组间处理变化的模型是过度拟合的；给出了一个类似的这里例子。如果您需要，可以强制lmer安装它：

m0 <- lme4::lmer(response~group+(group|ID), data=dd, REML=TRUE,
                 control=lmerControl(check.nobs.vs.nRE="ignore", calc.deriv=FALSE))

(请注意，我们还可以通过比较log-似然或REML准则-当我们有像这样的过度拟合模型时，不同的模型可以得到模型参数的不同线性组合的等价拟合。)

如果我们跑

library(lmerTest)
coef(summary(as(m1,"lmerModLmerTest"),ddf="Kenward-Roger"))["groupB",]

(默认的Satterthwaite计算在这里失败)

    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()，但是想要提取系数表的这一行)。

1. 现在让我们来试试简化的模型

m1 <- lme4::lmer(response~group+(1|ID), data=dd, REML=TRUE)

正如您注意到的，fit是单数的(如果您检查，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矩阵不是正定的.)

1. 正如陈刚在2011年张贴至r-sig-混合模型中指出的那样，我们可以通过允许群内两个观测值之间的负相关性(上面所示的加性模型只允许非负相关性)来拟合我们想要的模型：

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统计量的平方相匹配.

1. 我们也可以在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 )设置为相同)

系数表得到t-统计量(它称之为z值)正确，但没有“分母自由度”的概念，所以Z检验而不是t检验.

coef(summary(m4))$cond["groupB",]
            Estimate Std. Error  z value Pr(>|z|)
groupB      29.98126  0.7976217 37.58832        0

Answer 2

配对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

Answer 3

它们的斜率估计是由相应的“未配对”线性模型(lm(response~time)、lm(response2~time)和lm(response3~time))估计的。我遗漏了什么？

这是合理的，总斜率是相同的三个模型，是所有两两斜率的平均值(这将是很容易确定的)。我所缺少的是，在lm2.diff3的情况下，斜率估计的StdErr比lm2.diff1和lm2.diff2小一个数量级。lm2.diff3是在Response3上进行的，与Response1和Response2相比，每对观测的行为更为一致，因此其斜率估计的StdErr较小。