Post-Hoc Tests Exercises in R: 18 Tukey & Bonferroni Problems

Explanation: TukeyHSD() consumes an aov fit and returns a list keyed by each factor in the model. Each row is one pairwise difference with a Tukey-adjusted p-value, computed from the studentized range distribution. The adjustment is sharper than Bonferroni for all-pairwise comparisons because it knows the comparisons share a residual variance. Only trt2-trt1 survives at the 5% family-wise level.

Explanation: Only Diet 3 versus Diet 1 clears the family-wise 5% bar at day 21. Diet 4 versus Diet 1 looks promising on the raw difference (60.8 g) but the CI crosses zero once you spread the alpha across six comparisons. This is the classic Tukey lesson: large raw effects can fail post-hoc when the family of tests is broad. A common mistake is to run six separate t-tests at α = 0.05 and pick winners; that inflates Type I error to roughly 26%.

Explanation: The TukeyHSD output is a list of plain matrices, so subscripting works the standard way. The drop = FALSE argument is the subtle bit: without it, a single-row filter collapses to a numeric vector and loses the row name, breaking any downstream table-building code. This idiom matters whenever you slice a matrix by a condition that might return one row. An alternative for tidyverse users is to coerce with as.data.frame() and use dplyr::filter().

Explanation: Calling plot() on a TukeyHSD object dispatches to plot.TukeyHSD, which draws horizontal CIs with a reference line at zero. The las = 1 argument flips axis labels horizontal so they are readable. This visual is more honest than a p-value table because the eye sees both effect size and uncertainty. A common mistake is to use boxplot() instead, which shows distributions but hides the Tukey-adjusted family-wise interval.

Explanation: Bonferroni controls family-wise error by multiplying every raw p-value by $m$, the number of comparisons. With three pairs and a raw trt2-trt1 p of about 0.004, the Bonferroni-adjusted value is roughly 0.013, still significant. Bonferroni is simple and works for any test family, not just ANOVA pairs. Its weakness is power: when comparisons are correlated (as ANOVA pairs are, sharing a pooled SD), Bonferroni overcorrects and Tukey or Holm beats it.

Explanation: Holm sorts raw p-values from smallest to largest and multiplies the $i$-th by $m - i + 1$ instead of by $m$. The smallest raw p still gets multiplied by $m$ (so trt2-trt1 stays at 0.013), but later ones are multiplied by smaller integers, so they shrink. Net effect: the smallest p is unchanged, but middle and large p-values are less harshly penalized than Bonferroni. Holm is the default safe choice when you want FWER control without giving up power.

Explanation: Three different goals, three different answers. Bonferroni and Holm both control FWER (probability of any false positive); BH controls FDR (expected fraction of false positives among rejections). FWER is right when even one false claim is costly (drug approval); FDR is right when you screen many hypotheses and accept a small share of misses (genomics, A/B test dashboards). Holm dominates Bonferroni on power for FWER, and BH is dramatically less conservative when you have many tests.

Explanation: All three pairs use the same residual SD and sample sizes in this balanced design, so the widths are essentially identical (around 1.382). In unbalanced designs the widths differ because each pair uses its own pair-specific standard error scaled by the studentized range critical value. A common mistake is to assume the narrowest CI must contain the smallest p-value, but that is only true when the point estimates differ by similar amounts; precision and effect size are separate axes.

Explanation: Raw trt2-trt1 is 0.004; Tukey adjusts it to 0.012, Bonferroni to 0.013, Holm to 0.013. The trt2-ctrl raw p (0.088) is interesting: it looks borderline raw but is clearly not significant after any family-wise adjustment. Showing raw and adjusted side by side in a report is honest practice; never publish raw post-hoc p-values without their corrected partners, since a reviewer will assume FWER control is missing.

Explanation: At $k = 10$ the chance of at least one Type I error is already 40%, which is why running an uncorrected pairwise scan of even modest size is dangerous. The formula assumes independence, which post-hoc ANOVA comparisons are not, but it is a useful upper bound. The takeaway is that the family-wise rate grows fast even for modest $k$, so methods like Bonferroni, Holm, and Tukey are not optional ceremony; they are the price of doing many comparisons honestly.

Explanation: Treating dose as numeric would let aov() fit a linear trend instead of estimating a mean per dose level, which TukeyHSD() cannot then partition into pairwise contrasts. Converting to a factor is the required step. The supp:dose interaction matrix has 15 rows (6 choose 2 across the six supp x dose cells); some are within-supp dose comparisons, some are between-supp at the same dose, and some are cross-cell. Reading them needs care because not every row answers a clinically interesting question.

Explanation: Subsetting first and then running a one-way Tukey inside each supplement is one of two defensible ways to interrogate a significant interaction. The other is to use the full supp:dose matrix from Exercise 4.1 and filter to the rows you care about. The subset approach uses the within-supplement residual variance (smaller pool) so its CIs are slightly different from the full-model interaction CIs; report which one you used. The VC dose-response here is monotone: more is more, every step is significant.

Explanation: The double-square-bracket indexing ex_4_1[["supp:dose"]] is required because the component name contains a colon and would not survive $-style access. Tightening from 0.05 to 0.01 is a small but useful filter when interaction matrices are large and you want to lead the report with the rock-solid differences. The remaining seven rows above tell a clear story: every comparison involving the highest dose against the lowest dose passes; VC at low dose loses to OJ; the two supplements tie at the high dose.

Explanation: pairwise.wilcox.test() is the rank-based analogue of pairwise.t.test(). Use it when the residuals fail normality (Shapiro-Wilk rejects, QQ plot shows heavy tails) or when the variances are wildly unequal across groups. All three iris species differ on sepal width even after BH correction. A common pitfall is to use pairwise.t.test() then add a footnote saying "data was not normal"; if you suspect non-normality, switch the test, do not just note the violation.

Explanation: With only three pairs the BH and Bonferroni outputs are nearly identical here because the raw p-values are tiny. The differences become large when many tests sit in a borderline zone (raw p around 0.01 to 0.05) or when the family is large. BH ranks p-values and adjusts step-up; Bonferroni multiplies every raw p by $m$. For exploratory work with many groups, BH is almost always the right choice; for confirmatory work with stakes (clinical, regulatory), Bonferroni or Holm preserves the strict FWER guarantee.

Explanation: Wool A at low tension is the worst cell on average (most breaks), and most comparisons against that cell are significant or borderline. The QC actionable insight is to either avoid wool A at low tension or to recheck the operating procedure for that cell. The TukeyHSD output is the right tool here because the analyst wants every pairwise verdict at one go; running 15 separate t-tests uncorrected would inflate FWER to roughly 54%.

Explanation: In VC, every step up in dose produces a significant gain (2 beats 1, 1 beats 0.5, 2 beats 0.5), so no lower dose ties with the top dose; the recommendation defaults to 2.0 itself. The code pattern is portable: filter Tukey rows to those involving the reference level, look for the first non-significant tie, fall back gracefully when none exists. In a real submission you would also want to overlay effect sizes and CIs, not just p-values, before issuing a clinical recommendation.

Explanation: Merging Tukey and Bonferroni outputs requires care because the two functions return different shapes: Tukey gives a named matrix with one row per pair; pairwise.t.test gives a lower-triangular matrix indexed by group names. The bridge is a named lookup keyed by the Tukey row names. A common mistake is to rely on positional ordering, which silently breaks when factor levels are reordered. Showing both corrections in one table is honest reporting: readers can see whether the verdict depends on the choice of adjustment.