6
$\begingroup$

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)
$\endgroup$
6
$\begingroup$

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")))
$\endgroup$
  • 1
    $\begingroup$ +1. There was a follow-up question stats.stackexchange.com/questions/142685 asking how to do this in lmer. I searched around and found a hack by Ben Bolker allowing to do just that: stats.stackexchange.com/a/214007/28666. I tried and it works, I get the same value of t-statistic! I took the liberty to edit it into your answer, but if you don't want it to be here, please feel free to roll back with my apologies. $\endgroup$ – amoeba says Reinstate Monica Jun 29 '17 at 11:42
  • $\begingroup$ Actually, if one loads library('lmerTest'), then lmer output will contain p-values based on the Satterthwaite approximation, and the p-value will exactly correspond to the Welch's test! Pretty neat. Do you want to update your answer accordingly? $\endgroup$ – amoeba says Reinstate Monica Jun 29 '17 at 12:00
  • $\begingroup$ The Satterthwaite df are implemented in the package emmeans. contrast(emmeans(m2, specs="group")) replicates the df and p-value, at least to what looks like rounding error, where m2 is the gls model in Jake Westfall's answer $\endgroup$ – JWalker Oct 28 at 11:40
2
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ I'd suggest a slight modification: contrast(emmeans(m2, specs="group"), method = "pairwise"). The default method is "eff", which produces two contrasts for the two-level factor; "pairwise" will produce just one contrast, but the same t ratio and P value since the dimensionality is 1. $\endgroup$ – rvl Oct 29 at 15:34
  • $\begingroup$ yes, agree! This would be about the only time I use method="pairwise" --- I tend to use "revpairwise" because I want the sign of the contrast to be non-reference minus reference. $\endgroup$ – JWalker Nov 6 at 23:41
  • $\begingroup$ Ok. Revpairwise would be consistent with my comment. The point was that eff produces two contrasts that are the same in absolute value. $\endgroup$ – rvl Nov 7 at 0:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.