Long ago, when R's t.test had var.equal=TRUE by default, I wrote some
class notes showing that the result was equivalent to a one predictor
regression model. Because t.test does not default to var.equal=TRUE
these days, I'm curious to know if there is a way to specify weights
in an lm to obtain the same result as the Welch-adjusted values
reported by t.test at the current time. Is there a WLS equivalent
adjustment with lm?
Here's example code to show that lm is same as t.test with var.equal.
The signs come out differently, but otherwise the effect estimate,
standard error, t value are same:
set.seed(234234)
dat <- data.frame(x = gl(2, 50, labels = c("F", "M")))
dat$err <- rnorm(100, 0, 1)
dat$y <- ifelse(dat$x == "F", 40 + dat$err, 44 + dat$err)
m1 <- lm(y ~ x, dat)
summary(m1)
m1.t <- t.test(y ~ x, dat, var.equal = TRUE)
m1.t
## diff matches regression coefficient
(m1.t.effect <- diff(m1.t$estimate))
## standard error matches regression se
m1.t.stderr <- m1.t.effect/m1.t$statistic
If you run that, you see lm output:
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 39.9456 0.1180 338.65 <2e-16 ***
xM 3.9080 0.1668 23.43 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.8341 on 98 degrees of freedom
Multiple R-squared: 0.8485, Adjusted R-squared: 0.8469
F-statistic: 548.8 on 1 and 98 DF, p-value: < 2.2e-16
and t.test:
> m1.t <- t.test(y ~ x, dat, var.equal = TRUE)
> m1.t
Two Sample t-test
data: y by x
t = -23.427, df = 98, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-4.239038 -3.576968
sample estimates:
mean in group F mean in group M
39.94558 43.85358
> (m1.t.effect <- diff(m1.t$estimate))
mean in group M
3.908003
> m1.t.effect/m1.t$statistic
mean in group M
-0.1668129
--
Paul E. Johnson http://pj.freefaculty.org
Director, Center for Research Methods and Data Analysis http://crmda.ku.edu
To write to me directly, please address me at pauljohn at ku.edu.
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
