## Example from Sec. 13.13

# One-way ANOVA model:

y <- c(43.7,54.2,51.5,47.7,46.7,58.2,50.5,45.6,58.6,51.1,49.1,54.9,68.4,57.7,56.4)

# Making "sedentary" level 1 since that is the baseline level in R:
group.vec <- c('2A','2A','2A','2A','3W','3W','3W','3W','4J','4J','4J','4J','1S','1S','1S')

#creating dummy variables
x1 <- ifelse(group.vec=='2A',1,0)
x2 <- ifelse(group.vec=='3W',1,0)
x3 <- ifelse(group.vec=='4J',1,0)

full.mod <- lm(y ~ x1 + x2 +x3)

summary(full.mod)

reduced.mod <- lm(y ~ 1)

anova(reduced.mod, full.mod)

##########################################
# The explicit matrix calculations here:
##########################################

# The X matrix:
X <- cbind(rep(1,times=length(y)), x1, x2, x3)

# (X'X) matrix:
t(X) %*% X

# the least squares estimates:
(solve(t(X) %*% X)) %*% (t(X) %*% y)

################################################

# A short way to analyze this one-way ANOVA model:

anova(lm(y ~ factor(group.vec)))

# We see the same test-statistic value and P-value for the F-test.

coef(lm(y ~ factor(group.vec)))

# We see the same estimated coefficients.

################################################

# Comparing each pair of factor level means, using a Bonferroni correction:

pairwise.t.test(y, group.vec, p.adj = "bonf")

# This table shows the adjusted p-values for testing H0: mu_i = mu_j for all i not equal to j.
# Any p-value below the chosen alpha indicates that pair of levels 
# has significantly different mean responses.

################################################

# Alternatively, the Holm approach is an alteration of the Bonferroni method that
# is similar but has power uniformly as high or higher than the Bonferroni method:

pairwise.t.test(y, group.vec, p.adj = "holm")