机器算法验证 - 相当于 GLS 框架中的 Welch t 检验 - 吾爱随笔录

相当于 GLS 框架中的 Welch t 检验

机器算法验证 r 混合模式 t检验 lme4-nlme 广义最小二乘法

2022-03-24 09:32:19

Welch 的 t 检验如何表示为广义最小二乘模型？

标准独立样本 t 检验（假设被比较的样本来自方差相等的总体）可以表示如下：

Y_{i} = β_{0} + β_{1} X_{i} + ε_{i}

$Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i$

其中是结果，是对应于组成员资格的二元变量。的显着性检验将产生与标准独立样本 t 检验相同的 t 统计量。因此，下面的两个命令产生相同的统计数据（具有相同的自由度）： $Y$ $X$ $\beta_1$

t.test(extra~group, data = sleep, var.equal = TRUE)
lm(extra~group, data = sleep)

因为 Welch 的 t 检验允许被比较的样本之间的方差不相等，所以我的猜测是它等同于广义最小二乘法。那么问题是，什么调用gls（如果这确实是概念化问题的正确方法）会产生与以下相同的结果（包括自由度）：

t.test(extra~group, data = sleep, var.equal = FALSE)

2个回答

这是一个有趣的问题。需要注意的一件事是，如果组的大小不相等，则允许不等方差只会改变如果两组大小相等（即），则 Welch 的 -test （表示为）和 Student 的 -test （表示为）给出相同的检验统计量，因为 $t$ $n_1=n_2=n$ $t$ $t_w$ $t$ $t_s$

t_{w} = \frac{{\bar{y}}_{1} - {\bar{y}}_{2}}{\sqrt{\frac{s_{1}^{2}}{n_{1}} + \frac{s_{2}^{2}}{n_{2}}}} = \frac{{\bar{y}}_{1} - {\bar{y}}_{2}}{\sqrt{\frac{s_{1}^{2} + s_{2}^{2}}{n}}} = \frac{{\bar{y}}_{1} - {\bar{y}}_{2}}{\sqrt{(\frac{s_{1}^{2} + s_{2}^{2}}{2}) (\frac{2}{n})}} = t_{s}

$t_w= \frac{\bar{y}_1-\bar{y}_2}{\sqrt{\frac{s^2_1}{n_1}+\frac{s^2_2}{n_2}}}= \frac{\bar{y}_1-\bar{y}_2}{\sqrt{\frac{s^2_1+s^2_2}{n}}}= \frac{\bar{y}_1-\bar{y}_2}{\sqrt{(\frac{s^2_1+s^2_2}{2})(\frac{2}{n})}}= t_s$ 我指出这一点是因为您在帖子中提供的睡眠研究示例涉及相同的组大小，这就是为什么运行您的示例在所有情况下都会

t

$t$

无论如何，要回答您的问题，可以nlme::gls()通过使用weights与nlme::varIdent(). 下面我生成了一些具有不等组大小和不等方差的数据，然后展示如何使用 t.test 和回归函数（lm 或 gls）来拟合假设或不假设等方差的模型：

# generate data with unequal group sizes and unequal variances
set.seed(497203)
dat <- data.frame(group=rep.int(c("A","B"), c(10,20)),
  y = rnorm(30, mean=rep.int(c(0,1), c(10,20)), sd=rep.int(c(1,2),c(10,20))))

# the t-statistic assuming equal variances
t.test(y ~ group, data = dat, var.equal = TRUE)
summary(lm(y ~ group, data = dat))

# the t-statistic not assuming equal variances
t.test(y ~ group, data = dat, var.equal = FALSE)
library(nlme)
summary(gls(y ~ group, data = dat, weights=varIdent(form = ~ 1 | group)))

# a hack to achieve the same thing in lmer
#    (lmerControl options are needed to prevent lmer from complaining
#    about too many levels of the grouping variable)
dat <- transform(dat,
                 obs=factor(1:nrow(dat)),
                 dummy=as.numeric(group=="B"))
library('lme4')
summary(lmer(y ~ group + (dummy-1|obs), data=dat,      
             control=lmerControl(check.nobs.vs.nlev = "ignore",
                                 check.nobs.vs.nRE  = "ignore")))

您还询问了获得相同的自由度。自由度基于Satterthwaite 近似，并t.test默认应用该近似，因为这是 Welch 描述的解决方案的一部分。但gls不这样做。理论上这是可以做到的，我相信会这样做，所以你应该能够PROC MIXED在. 也许（可能）有一些 R 包可以很容易地为一般回归模型（具有连续预测器）获得 Satterthwaite DF，但我不知道它是什么。SASPROC MIXED

@amoeba 更新： Satterthwaite 近似值作为包中的默认值实现lmerTest，因此要获得 $p$ -value 与 Welch 的 t 检验完全匹配，可以运行：

library('lmerTest')
summary(lmer(y ~ group + (dummy-1|obs), data=dat,      
             control=lmerControl(check.nobs.vs.nlev = "ignore",
                                 check.nobs.vs.nRE  = "ignore")))

如果 m2 是 Jake Westfall 的答案中的 gls 对象，则使用 emmeans 包中的 contrast(emmeans(m2)) 计算 Satterthwaite df 和相关的 p 值。在 Jake 的示例中，未调整和调整后的 df 非常相似，因此 p 值和任何解释实际上是相同的。这是一个重要的示例（我使用较小的 n 并反转具有较大方差的组）。

library(nlme)
library(emmeans)
set.seed(497203)
n1 <- 8
n2 <- 4
dat <- data.frame(group=rep.int(c("A","B"), c(n1,n2)),
  y = rnorm(n1+n2, mean=rep.int(c(0,1), c(n1,n2)), sd=rep.int(c(1,2),c(n1,n2))))

# the t-statistic assuming equal variances
t.student <- t.test(y ~ group, data = dat, var.equal = TRUE)
m1 <- lm(y ~ group, data = dat)

# the t-statistic not assuming equal variances
t.welch <- t.test(y ~ group, data = dat, var.equal = FALSE)
m2 <- gls(y ~ group, data = dat, 
          weights=varIdent(form = ~ 1 | group))

m2.contrast <- contrast(emmeans(m2, specs="group"))

将每个中的 df、t 和 p 收集到一个表中给出：

      method               df                 t                  p
1: Student t               10  -2.3012090076821 0.0441633525165716
2:        lm               10   2.3012090076821 0.0441633525165716
3:   Welch t 3.91598345952776 -1.87308830595972  0.135881711655436
4:       gls               10  1.87308831515667 0.0905471567272453
5:   emmeans  3.9158589862009 -1.87308831515667   0.13588402534431

其它你可能感兴趣的问题

上一篇扩展 Cox 模型和 cox.zph 下一篇ImageNet：前五位错误是什么意思？