function (alpha = 0.05) 
{
    diets = unique(diet)
    mu.vec = NULL
    nvec = NULL
    mean.vec = NULL
    for (i in 1:length(diets)) {
        mu.vec = c(mu.vec, mean(ctime[diet == diets[i]]))
        nvec = c(nvec, length(ctime[diet == diets[i]]))
        mean.vec = c(mean.vec, rep(mu.vec[i], nvec[i]))
    }
    tr = length(nvec)
    N = sum(nvec)
    MSE = sum((ctime - mean.vec)^2/(N - tr))
    s = sqrt(MSE)
    intervals = NULL
    intervalsLSD = NULL
    intervalsBonf = NULL
    intervalsScheffe = NULL
    for (i in 1:3) {
        for (j in (i + 1):4) {
            nijstar = 1/(0.5 * (1/nvec[i] + 1/nvec[j]))
            qTK = qtukey(1 - alpha, tr, N - tr)
            Diff = mu.vec[i] - mu.vec[j]
            lower = Diff - qTK * s/sqrt(nijstar)
            upper = Diff + qTK * s/sqrt(nijstar)
            lowerLSD = Diff - qt(1 - alpha/2, N - tr) * s * sqrt(1/nvec[i] + 
                1/nvec[j])
            upperLSD = Diff + qt(1 - alpha/2, N - tr) * s * sqrt(1/nvec[i] + 
                1/nvec[j])
            lowerBonf = Diff - qt(1 - alpha/(2 * 6), N - tr) * 
                s * sqrt(1/nvec[i] + 1/nvec[j])
            upperBonf = Diff + qt(1 - alpha/(2 * 6), N - tr) * 
                s * sqrt(1/nvec[i] + 1/nvec[j])
            lowerScheffe = Diff - sqrt((tr - 1) * qf(1 - alpha, 
                tr - 1, N - tr)) * s * sqrt(1/nvec[i] + 1/nvec[j])
            upperScheffe = Diff + sqrt((tr - 1) * qf(1 - alpha, 
                tr - 1, N - tr)) * s * sqrt(1/nvec[i] + 1/nvec[j])
            intervals = rbind(intervals, c(lower, upper))
            intervalsLSD = rbind(intervalsLSD, c(lowerLSD, upperLSD))
            intervalsBonf = rbind(intervalsBonf, c(lowerBonf, 
                upperBonf))
            intervalsScheffe = rbind(intervalsScheffe, c(lowerScheffe, 
                upperScheffe))
        }
    }
    cols = c("orange", "magenta", "blue", "black")
    m = min(intervals, intervalsLSD, intervalsBonf, intervalsScheffe)
    M = max(intervals, intervalsLSD, intervalsBonf, intervalsScheffe)
    plot(0, 0, xlab = "", ylab = "", type = "n", axes = F, xlim = c(0, 
        7), ylim = c(-15, 15))
    Delta = 0.1
    for (i in 1:6) {
        segments(i - 1.5 * Delta, intervals[i, 1], i - 1.5 * 
            Delta, intervals[i, 2], col = cols[1])
        segments(i - 0.5 * Delta, intervalsLSD[i, 1], i - 0.5 * 
            Delta, intervalsLSD[i, 2], col = cols[2])
        segments(i + 0.5 * Delta, intervalsBonf[i, 1], i + 0.5 * 
            Delta, intervalsBonf[i, 2], col = cols[3])
        segments(i + 1.5 * Delta, intervalsScheffe[i, 1], i + 
            1.5 * Delta, intervalsScheffe[i, 2], col = cols[4])
    }
    axis(2)
    mtext(expression(mu[1] - mu[2]), 1, 1, at = 1)
    mtext(expression(mu[1] - mu[3]), 1, 2, at = 2)
    mtext(expression(mu[1] - mu[4]), 1, 1, at = 3)
    mtext(expression(mu[2] - mu[3]), 1, 2, at = 4)
    mtext(expression(mu[2] - mu[4]), 1, 1, at = 5)
    mtext(expression(mu[3] - mu[4]), 1, 2, at = 6)
    abline(h = seq(-15, 15, 1), col = "grey")
    abline(h = 0, lwd = 2, col = "red")
    title(substitute("Pairwise Comparisons of Means (Coagulation Data): " ~ 
        1 - alpha == conf, list(conf = 1 - alpha)))
    legend(c(3, 6.5), c(-7, -15), c("Tukey-Kramer pairwise comparisons", 
        "Fisher's protected LSD", "Bonferroni intervals", "Scheffe's intervals for all contrasts"), 
        col = cols, lty = rep(1, 4), bg = "white")
    out = list(intervals = intervals, intervalsLSD = intervalsLSD, 
        intervalsBonf = intervalsBonf, intervalsScheffe = intervalsScheffe)
    out = cbind(intervals, intervalsLSD, intervalsBonf, intervalsScheffe)
    round(out, 2)
}