Three versions of the independent two sample t-test (and R)

Question

This post concerns three versions of the independent two sample t-test:

1. Student's t-test uses a pooled standard deviation in the denominator (all equations are shamelessly copied from Wikipedia - Student's t-test): $$ {\displaystyle t={\frac {{\bar {X}}_{1}-{\bar {X}}_{2}}{s_{p}\cdot {\sqrt {{\frac {1}{n_{1}}}+{\frac {1}{n_{2}}}}}}}}$$ where: $$ {\displaystyle s_{p}={\sqrt {\frac {\left(n_{1}-1\right)s_{X_{1}}^{2}+\left(n_{2}-1\right)s_{X_{2}}^{2}}{n_{1}+n_{2}-2}}}} $$ and: $$d.f. = n_1 + n_2 -2 $$ It is well known that Student's t-test can be seriously biased for unequal sample sizes and variances (Welch 1938 demonstrated this analytically, before computers and the many published simulation studies on this subject).

2. Next is "an alternative criterion that has often been employed," which is less biased for unequal sample sizes and variances than is Student's t (Welch 1938). This version uses the standard error of the difference in the denominator: $$ {\displaystyle t={\frac {{\bar {X}}_{1}-{\bar {X}}_{2}}{s_{\bar {\Delta }}}}} $$ where: $$ {\displaystyle s_{\bar {\Delta }}={\sqrt {{\frac {s_{1}^{2}}{n_{1}}}+{\frac {s_{2}^{2}}{n_{2}}}}}.} $$ Note that if $ n_1 = n_2 $ or $ s_{X_{1}}^{2} = s_{X_{2}}^{2} $ then Student's t-test is identical to this alternative (I'll let the reader do the algebra).

3. Welch (1938) presented an approximate degrees of freedom for the above statistic that further reduces bias (see Wikipedia - Welch's t-test for the formula).

Edit: I replaced the example with one that better demonstrates a Type I error.

The t.test function in R uses versions 1 or 3 above (with var.equal = TRUE or FALSE, respectively). Here's an example (see this post) where Student's and the alternative t-test are calculated from the data and then compared with the t.test function:

set.seed(4)
n1 <- 15
n2 <- 50
dat1 <- rnorm(n1, 100, 20)
dat2 <- rnorm(n2, 100, 5)
m1 <- mean(dat1)
m2 <- mean(dat2)
v1 <- var(dat1)
v2 <- var(dat2)
vp <- ((n1-1)*v1 + (n2-1)*v2) / (n1+n2-2)
t.S <- (m1-m2)/sqrt(vp/n1 + vp/n2) # Student's t (not available in R {stats})
df.S <- n1 + n2 -2 # same for Student's and Alternative t
t.A <- (m1-m2)/sqrt(v1/n1 + v2/n2) # "Alternative" t (for comparison with R)
t.R <- t.test(dat1, dat2, var.equal=TRUE) # R t.test
t.W <- t.test(dat1, dat2, var.equal=FALSE) # Welch t.test
cat(" Student's t = ", round(t.S, 3), "d.f. =" , df.S, "p.value = ",     round(2*pt(abs(t.S), df=df.S, lower.tail=F), 3), "\n",
    "Alternative t = ", round(t.A, 3), "d.f. =" , df.S, "p.value = ", round(2*pt(abs(t.A), df=df.S, lower.tail=F), 3), "\n",
    "R (var.equal=TRUE) t = ", round(t.R$statistic, 3), "d.f. =" , t.R$parameter, "p.value = ", round(t.R$p.value, 3), "\n",
    "R (var.equal=FALSE) t = ", round(t.W$statistic, 3), "d.f. =" , round(t.W$parameter, 1), "p.value = ", round(t.W$p.value, 3), "\n")

Here is the output:

Student's t =  2.956 d.f. = 63 p.value =  0.004 
Alternative t =  1.768 d.f. = 63 p.value =  0.082 
R (var.equal=TRUE) t =  2.956 d.f. = 63 p.value =  0.004 
R (var.equal=FALSE) t =  1.768 d.f. = 14.6 p.value =  0.098  

See that R is calculating the t statistics correctly. The parameters were chosen to demonstrate the serious bias that can occur when using Student's t with unequal sample sizes and variances.

Also of concern is that a linear model/ANOVA is equivalent to Student's t-test:

dat <- data.frame(y=c(dat1, dat2), g=as.factor(c(rep("A", length(dat1)), rep("B", length(dat2)))))
summary(lm(y~g, dat))
summary(aov(y~g, dat))

Here is the (abbreviated) output:

> summary(lm(y~g, dat))
Coefficients:
            Estimate Std. Error t value Pr(>|t|)     
(Intercept)  108.827      2.388  45.569  < 2e-16 ***
gB            -8.049      2.723  -2.956  0.00438 **  

> summary(aov(lm(y~g, dat)))
            Df Sum Sq Mean Sq F value Pr(>F)  
g            1    747   747.5   8.737 0.00438 **
Residuals   63   5390    85.6                  

Note that $ F = t^2 $ for independent samples. In a subsequent paper, Welch (1951) proposed a more robust version of ANOVA (the oneway.test function in R).

OK, some questions:

• What is the source of the "alternative criterion," for which Welch (1938) did not provide a reference? And, what should we call it?
• Why is the standard error of the difference superior to the pooled s.d.?
• It seems that ANOVA (LM) is robust only when sample sizes or variances are equal across the groups? I should read Welch's other paper (1951).

Welch, B.L. (1938). The significance of the difference between two means when the population variances are unequal. Biometrika 29: 350–362.

Welch, B.L. (1951). On the comparison of several mean values: an alternative approach. Biometrika 38: 330–336.

Answer 1

Finally! I found this answer in Lumley et. al. (2002):

• The "unequal variance t-statistic" (i.e. the "alternative criterion") can easily be derived by applying the Central Limit Theorem (see the image at the bottom of this answer and note that n1 is missing in the SE in the second line).
• The CLT guarantees that the unequal variance t-statistic is normally distributed with variance equal to one for large samples. It's presented as a "Z-test" in some textbooks.
• For small samples, the unequal variance t-statistic doesn't exactly follow a t-distribution, even if the outcomes are Normal. Approximate degrees of freedom (Welch 1938) can be used for small samples.
• The CLT guarantees that Student's t (the "equal variance t-statistic") is normally distributed for large samples however the variance may not equal one, i.e. Student's t is sensitive to unequal variances when sample sizes are unequal. Likewise, LMs are sensitive to heteroscedasticity.

Lumley et. al. (2002) concluded that heteroscedasticity is a more serious problem for LMs than is Normality.

I suggest that more flexible models be applied when there is strong heteroscedasticity, including Generalised Least Squares (where the analyst specifies the form of the heteroscedasticity) and Generalised Linear Models (where the error distribution and link function define the heteroscedasticity).

The traditional approach of transforming the response variable and applying LMs can "muddy" the results: it can be difficult to explain estimates on the transformed scale (e.g. what are log dollars?) and even more difficult to back-transform estimates and associated confidence intervals to the original scale. Transformation can also bias the results (which I am currently studying).

Here's GLS applied to the same data in the question:

library(nlme)
summary(gls(y~g, dat, weights=varIdent(form = ~ 1 | g)))

Here is the (abbreviated) output:

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | g 
 Parameter estimates:
        A         B 
1.0000000 0.2763807 

Coefficients:
                Value Std.Error   t-value p-value
(Intercept) 108.82652  4.500116 24.183046  0.0000
gB           -8.04867  4.551386 -1.768399  0.0818

The weight for group B in these results equals the ratio of standard deviations, i.e. sd(dat2)/sd(dat1). The p-value is the same as that for the unequal variance t-test with d.f. = n1 + n2 - 2 (which shouldn't be surprising).

Knowing that the Normality is not an important assumption for large samples and that autocorrelation and/or heteroscedasticity are common, the GLS framework could be used more frequently (e.g. fitting a GLS model to count data with autocorrelated errors could be easier than fitting a GLM with autocorrelated errors).

Lumley T, Diehr P, Emerson S, Chen L (2002) The importance of the normality assumption in large public health data sets. Annu Rev Public Health 23: 151–169.