Equivalent to Welch's t-test in GLS framework

Question

How can Welch's t-test be expressed as a generalized least squares model?

A standard independent samples t-test (where it is assumed that the samples being compared are drawn from populations with equal variance) can be expressed as follows:

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

where $Y$ is the outcome and $X$ is a binary variable corresponding to group membership. The significance test of $\beta_1$ will produce the same t statistic as the standard independent samples t-test. Thus, the two commands below produce the same statistics (with the same degrees of freedom):

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

Because a Welch's t-test allows for unequal variances between the samples being compared, my guess is that it would be equivalent to a generalized least squares procedure. The question then is, what call to gls (if this indeed the correct way of conceptualizing the problem) would produce the same results (including degrees of freedom) as:

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

Answer 1

It's an interesting question. One thing to note is that allowing unequal variances will only change the $t$-statistic if the groups are of unequal size. If the two groups are of equal size (i.e., $n_1=n_2=n$), then Welch's $t$-test (denoted $t_w$) and Student's $t$-test (denoted $t_s$) give the same test statistic, since $$ 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 $$ I point this out because the sleep study example you give in your post involves equal group sizes, which is why running your example returns the same $t$-statistic in all cases.

Anyway, to answer your question, this can be done in nlme::gls() by using the weights parameter combined with nlme::varIdent(). Below I generate some data with unequal group sizes and unequal variances, then show how to fit the models assuming or not assuming equal variance, using both t.test and a regression function (lm or 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")))

You also asked about getting the same degrees of freedom. The degrees of freedom are based on the Satterthwaite approximation, and t.test applies the approximation by default as that is part of the solution described by Welch. But gls does not do so. Theoretically this could be done, and I believe PROC MIXED in SAS will do so, so you should be able to reproduce the results exactly in PROC MIXED. Maybe (probably) there is some R package that will make it easy to get the Satterthwaite DFs for general regression models (with continuous predictors), but I don't know what it is.

Update by @amoeba: The Satterthwaite approximation is implemented as the default one in the lmerTest package, so to get $p$-value exactly matching Welch's t-test one can run:

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

Answer 2

if m2 is the gls object from Jake Westfall's answer, then the Satterthwaite df and associated p-value are computed using contrast(emmeans(m2)) from the emmeans package. In Jake's example, the unadjusted and adjusted df are very similar, so the p-value and any interpretation is effectively the same. Here is an example where it matters (I use a smaller n and reverse the group with the larger variance).

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"))

Gathering the df, t, and p from each into a single table gives:

      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