obs <- c(23, 24, 25, 27, 28, 26, 31, 32, 29, 24, 23, 28, 38, 36, 35,
    34, 36, 39, 36, 35, 36, 34, 38, 39, 33, 34, 35, 37, 39, 35, 34,
    38, 36, 34, 36, 31, 28, 24, 27, 35, 35, 34, 26, 27, 25, 26, 29,
    25, 36, 37, 34, 28, 26, 24)
cycletime <- gl(3, 18, labels = c("40", "50", "60"))
temperature <- gl(2, 9 ,54, labels = c("300", "350"))
operator <- gl(3, 3, 54, labels = c("1", "2", "3"))

fitmodel <- lm(obs ~ cycletime * temperature * operator)
anova(fitmodel)

boxplot(obs ~ cycletime)
boxplot(obs ~ operator)
boxplot(obs ~ temperature)

plot(aggregate(obs ~ as.character(cycletime), FUN = mean), type = "b", ylab = "mean of obs")
plot(aggregate(obs ~ as.character(operator), FUN = mean), type = "b" , ylab = "mean of obs")
plot(aggregate(obs ~ as.character(temperature), FUN = mean), type = "b", ylab = "mean of obs")

interaction.plot(cycletime, temperature, obs)
interaction.plot(cycletime, operator, obs)
interaction.plot(operator, temperature, obs)

res <- residuals(fitmodel)
qqnorm(res)
qqline(res)