One-Way Between-Subjects ANOVA

Author

Dr. Anisha Singh

The dataset for this week is readtime.csv. We will use this data to determine if reading time (in milliseconds) differs by type of device (1 = phone, 2 = tablet, 3 = laptop).

Packages

# install.packages("psych")       # run once if needed
# install.packages("lawstat")
# install.packages("effectsize")
# install.packages("emmeans")
# install.packages("agricolae")

library(psych)
library(lawstat)
library(effectsize)
library(emmeans)
library(agricolae)

Load and Prepare the Data

reader <- read.csv("readtime.csv")
reader$Device <- as.factor(reader$Device)

Exploring the Data

Before fitting any model, get group-level descriptive statistics using describeBy().

describeBy(reader$ReadingTime, reader$Device)

Descriptive statistics by group
group: 1
   vars  n    mean     sd  median trimmed    mad    min     max   range  skew kurtosis    se
X1    1 20 1252.83 305.62 1233.12 1269.18 284.36 572.08 1797.21 1225.13 -0.29    -0.34 68.34
------------------------------------------------------------
group: 2
   vars  n    mean     sd  median trimmed    mad    min     max   range skew kurtosis    se
X1    1 20 1032.20 307.44 1085.92 1026.23 330.37 551.34 1588.64 1037.30 0.02    -1.18 68.75
------------------------------------------------------------
group: 3
   vars  n    mean     sd  median trimmed    mad    min     max range skew kurtosis    se
X1    1 20 1014.00 288.58  969.52 1009.07 244.22 543.75 1530.75   987 0.21    -1.16 64.53

Phone users had the longest average reading time (M = 1252.83 ms, SD = 305.62), followed by tablet (M = 1032.20 ms, SD = 307.44) and laptop (M = 1014.00 ms, SD = 288.58).

Checking Assumptions

Normality

One-way ANOVA assumes that the residuals are normally distributed. We test this with a Shapiro-Wilk test on the model residuals (after fitting the model below) rather than on raw scores within each group.

fit2 <- aov(ReadingTime ~ Device, data = reader)
shapiro.test(fit2$residuals)

    Shapiro-Wilk normality test

data:  fit2$residuals
W = 0.98523, p-value = 0.6822

The test is non-significant (p = .682), so we have no evidence against normality.

Homogeneity of Variance

ANOVA assumes equal variances across groups (homoscedasticity). We test this with Levene’s test.

levene.test(reader$ReadingTime, reader$Device)

    Modified robust Brown-Forsythe Levene-type test based on the
    absolute deviations from the median

data:  reader$ReadingTime
Test Statistic = 0.061395, p-value = 0.9405

The test is non-significant (p = .941), confirming equal variances. The standard ANOVA is appropriate.

Note

If Levene’s test is significant, use a Welch one-way ANOVA instead: oneway.test(ReadingTime ~ Device, data = reader, var.equal = FALSE) This does not assume equal variances, analogous to Welch’s correction for the independent-samples t-test.

Fitting the Model

fit1 <- lm(ReadingTime ~ Device, data = reader)
anova(fit1)

Analysis of Variance Table

Response: ReadingTime
          Df  Sum Sq Mean Sq F value  Pr(>F)
Device     2  706968  353484  3.9102 0.02562 *
Residuals 57 5152810   90400

Device type has a statistically significant effect on reading time, F(2, 57) = 3.91, p = .026.

Alternative syntax: `aov()`

Some post-hoc functions require an aov object. The results are identical.

fit2 <- aov(ReadingTime ~ Device, data = reader)
summary(fit2)

            Df  Sum Sq Mean Sq F value Pr(>F)
Device       2  706968  353484    3.91 0.0256 *
Residuals   57 5152810   90400

ANOVA and regression are the same model

summary(fit1) returns regression output. The F-statistic and degrees of freedom are identical; the coefficients represent mean differences from the reference group (Device 1).

summary(fit1)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  1252.83      67.23  18.635   <2e-16 ***
Device2      -220.63      95.08  -2.320   0.0239 *
Device3      -238.83      95.08  -2.512   0.0149 *

Residual standard error: 300.7 on 57 degrees of freedom
Multiple R-squared:  0.1206,  Adjusted R-squared:  0.08979
F-statistic:  3.91 on 2 and 57 DF,  p-value: 0.02562

Device 2 (tablet) read 220.63 ms faster than Device 1 (phone), and Device 3 (laptop) read 238.83 ms faster.

Overall Effect Size: Eta Squared

Before drilling into which pairs differ, quantify the overall effect of device type. Eta squared (η²) is the proportion of total variance in reading time explained by device.

eta_squared(fit1, partial = FALSE)

# Effect Size for ANOVA (Type I)

Parameter | Eta2 |       95% CI
-------------------------------
Device    | 0.12 | [0.01, 1.00]

η² = .12 — device type explains about 12% of the variance in reading time, a medium effect by Cohen’s conventions (small ≈ .01, medium ≈ .06, large ≈ .14).

Note

η² vs. partial η²: In a one-way design these are identical. In factorial designs (Week 10+) they diverge — partial η² measures the effect of one factor after removing the other factors, so values from different studies are not directly comparable. Always report which you computed.

Post-hoc Tests

A significant ANOVA tells us that at least one group mean differs — but not which ones. Post-hoc tests make all pairwise comparisons while controlling the familywise Type I error rate.

Tukey’s HSD (Recommended)

Tukey’s Honest Significant Difference uses the studentized range distribution to test all pairwise comparisons simultaneously. It is designed specifically for this purpose and is the standard recommendation for balanced or near-balanced designs.

TukeyHSD(fit2)

  Tukey multiple comparisons of means
    95% family-wise confidence level

$Device
         diff       lwr        upr     p adj
2-1 -220.6260 -449.4258   8.173844 0.0610
3-1 -238.8265 -467.6263 -10.026656 0.0389 *
3-2  -18.2005 -247.0003 210.599344 0.9800

Laptop users read significantly faster than phone users (mean difference = −238.83 ms, p = .039). The tablet vs. phone comparison approached but did not reach significance (p = .061). Tablet and laptop did not differ (p = .980).

Confidence Intervals for Tukey Pairwise Differences

TukeyHSD() already includes confidence intervals in its output. We can also obtain them using the emmeans workflow, which integrates cleanly with pairwise comparisons and effect size calculations.

em_fit <- emmeans(fit2, ~ Device)
tukey_results <- pairs(em_fit, adjust = "tukey")
confint(tukey_results, adjust = "tukey")

 contrast          estimate   SE df lower.CL upper.CL
 Device1 - Device2    220.6 95.1 57    -8.17      449
 Device1 - Device3    238.8 95.1 57    10.03      468
 Device2 - Device3     18.2 95.1 57  -210.60      247

Confidence level used: 0.95
Conf-level adjustment: tukey method for comparing a family of 3 estimates

The 95% CI for Device 1 − Device 3 [10.03, 468] does not include 0, confirming a significant difference. The CI for Device 1 − Device 2 [−8.17, 449] just crosses 0, consistent with the marginal p = .061.

Pairwise Effect Sizes: Cohen’s d

Statistical significance tells us whether groups differ; effect sizes tell us how much. Cohen’s d is a standardized mean difference. For post-hoc comparisons after ANOVA, we use the pooled error term from the model (rather than just the two groups’ SDs), which is more stable and consistent across comparisons.

d_results <- eff_size(emmeans(fit2, pairwise ~ Device),
                      sigma = sigma(fit2),
                      edf = df.residual(fit2),
                      ci = TRUE)
d_results

 contrast            effect.size    SE df lower.CL upper.CL
 (Device1 - Device2)      0.7338 0.324 57   0.0858    1.382
 (Device1 - Device3)      0.7943 0.325 57   0.1438    1.445
 (Device2 - Device3)      0.0605 0.316 57  -0.5728    0.694

sigma used for effect sizes: 300.7

Phone vs. tablet: d = 0.73 (medium-large); phone vs. laptop: d = 0.79 (large); tablet vs. laptop: d = 0.06 (negligible). The significant Tukey comparison (phone vs. laptop) corresponds to a large effect.

Bonferroni Correction

Warning

Planned vs. exploratory comparisons

There are two fundamentally different situations for multiple comparisons:

Planned (a priori) comparisons: You specified exact hypotheses before collecting data. For example, “I predict phone will differ from laptop.” In this case, Bonferroni is appropriate — you divide α by the number of tests you committed to in advance.
Exploratory (post-hoc) comparisons: You run all pairwise tests after finding a significant omnibus ANOVA, scanning for differences. This is the more common situation, and it is why we need stricter corrections. Bonferroni is overly conservative here — it was designed for a small set of independent tests, not the correlated structure of all-pairwise comparisons. Use Tukey’s HSD instead, which is purpose-built for this scenario and offers better power with the same familywise error control.

emmeans(fit2, pairwise ~ Device, adjust = "bonferroni")

$contrasts
 contrast          estimate   SE df t.ratio p.value
 Device1 - Device2    220.6 95.1 57   2.320  0.0718
 Device1 - Device3    238.8 95.1 57   2.512  0.0446
 Device2 - Device3     18.2 95.1 57   0.191  1.0000

P value adjustment: bonferroni method for 3 tests

The p-values shown are already adjusted — compare them to α = .05, not .0167. With Bonferroni, no comparison reaches significance. Compare this to Tukey’s HSD, which identified phone vs. laptop as significant (p = .039): Tukey is the more appropriate tool for exploratory post-hoc comparisons, and it has better power.

To get confidence intervals for the Bonferroni-adjusted comparisons:

bonf_results <- pairs(em_fit, adjust = "bonferroni")
confint(bonf_results, adjust = "bonferroni")

 contrast          estimate   SE df lower.CL upper.CL
 Device1 - Device2    220.6 95.1 57    -13.9      455
 Device1 - Device3    238.8 95.1 57      4.3      473
 Device2 - Device3     18.2 95.1 57   -216.3      253

Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 3 estimates

Notice that the Bonferroni CI for Device 1 − Device 3 [4.3, 473] is slightly wider than the Tukey CI [10.03, 468], and its lower bound is closer to 0. This reflects the added conservatism of the Bonferroni correction — it requires a larger observed difference before excluding 0 from the interval.

Fisher’s LSD

Fisher’s LSD performs standard t-tests between all pairs using the pooled MSE. It appears frequently in older published literature but does not adequately control the familywise error rate when there are more than three groups. You may encounter it in papers, but Tukey’s HSD is preferred for new analyses.

# install.packages("agricolae")   # run once if needed
library(agricolae)
LSD.test(fit2, "Device", console = TRUE)

Treatments with the same letter are not significantly different.

  ReadingTime groups
1    1252.825      a
2    1032.199      b
3    1013.999      b