1

Does anyone know if it is possible to use lmFit or lm in R to calculate a linear model with categorical variables while including all possible comparisons between the categories? For example in the test data created here:

set.seed(25)
f <- gl(n = 3, k = 20, labels = c("control", "low", "high"))
mat <- model.matrix(~f, data = data.frame(f = f))
beta <- c(12, 3, 6)  #these are the simulated regression coefficient
y <- rnorm(n = 60, mean = mat %*% beta, sd = 2)
m <- lm(y ~ f)

I get the summary:

summary(m)
Call:
lm(formula = y ~ f)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.3505 -1.6114  0.1608  1.1615  5.2010 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  11.4976     0.4629  24.840  < 2e-16 ***
flow          3.0370     0.6546   4.639 2.09e-05 ***
fhigh         6.1630     0.6546   9.415 3.27e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.07 on 57 degrees of freedom
Multiple R-squared:  0.6086,    Adjusted R-squared:  0.5949 
F-statistic: 44.32 on 2 and 57 DF,  p-value: 2.446e-12

which is because the contrasts term ("contr.treatment") compares "high" to "control" and "low" to "control".

Is it possible to get also the comparison between "high" and "low"?

duplode
  • 33,731
  • 7
  • 79
  • 150
Makis Ladoykakis
  • 111
  • 2
  • 8

1 Answers1

2

If you use aov instead of lm, you can use the TukeyHSD function from the stats package:

fit <- aov(y ~ f)
TukeyHSD(fit)
#   Tukey multiple comparisons of means
#     95% family-wise confidence level

# Fit: aov(formula = y ~ f)

# $f
#                  diff      lwr      upr    p adj
# low-control  3.036957 1.461707 4.612207 6.15e-05
# high-control 6.163009 4.587759 7.738259 0.00e+00
# high-low     3.126052 1.550802 4.701302 3.81e-05

If you want to use an lm object, you can use the TukeyHSD function from the mosaic package:

library(mosaic)
TukeyHSD(m)

Or, as @ben-bolker suggests,

library(emmeans)
e1 <- emmeans(m, specs = "f")
pairs(e1)
#  contrast        estimate        SE df t.ratio p.value
#  control - low  -3.036957 0.6546036 57  -4.639  0.0001
#  control - high -6.163009 0.6546036 57  -9.415  <.0001
#  low - high     -3.126052 0.6546036 57  -4.775  <.0001

# P value adjustment: tukey method for comparing a family of 3 estimates

With lmFit:

library(limma)
design <- model.matrix(~0 + f)
colnames(design) <- levels(f)
fit <- lmFit(y, design)
contrast.matrix <- makeContrasts(control-low, control-high, low-high,
                                 levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
round(t(rbind(fit2$coefficients, fit2$t, fit2$p.value)), 5)
#                    [,1]     [,2]  [,3]
# control - low  -3.03696 -4.63938 2e-05
# control - high -6.16301 -9.41487 0e+00
# low - high     -3.12605 -4.77549 1e-05

Also see Multiple t-test comparisons for more information.

Weihuang Wong
  • 12,868
  • 2
  • 27
  • 48