To compare the mean of every group with every other group, the contrasts will not be orthogonal. The Tukey method controls for the inflation of Type I errors. $q_{j, k} = \frac{\bar Y_j - \bar Y_k}{\sqrt{\hat \sigma^2 / r}} \sim SRD(J, n-J)$ where $\hat \sigma^2 = RSS/(n-J)$, $r$ is the number of replications, and SRD is the Studentized Range Distribution (the [[cumulative density function|cdf]] given by `ptukey` in [[R]]). In R, use the linear regression method to conduct an ANOVA, store the table, and pass to the `TukeyHSD` method. ```R lmod = lm(y ~ xm, data=data) av = aov(lmod) TukeyHSD(av) ``` You can compare the p-values to the standard significance level (in contrast with the Bonferroni method). A p-value less than your alpha indicates a true difference. To compute this from scratch ```R lmod = lm(y ~ x, data=data) n = length(data$x) J = length(unique(data$x)) rss = sum(resid(lmod)^2) sighat = sqrt(rss/(n-J)) ybar <- tapply(data$y, data$x, mean) diff <- as.numeric(c(ybar[2] - ybar[1], ybar[3] - ybar[1], ybar[3] - ybar[1])) qjk = abs(diff)/(sighat/J) p = round(1-ptukey(qjk, nmeans=J, df=n-K), 7) q = qtukey(1-alpha, nmeans = J, df = n-J) ci = data.frame(cbind((diff) - q*(sighat/3), (diff) + q*(sighat/3))) names(ci) = c("lwr", "uppr") cbind(diff, ci, p) ``` If the number of replications for each group is not the same, use the Tukey-Cramér method.