# RPD -- Example 8.3 -- Algae Density: CO2 only, Both Reps

algae <- data.frame(day = rep(1:14,2),
                    density = c(0.530, 1.183, 1.603, 1.994, 2.708,
                                3.006, 3.867, 4.059, 4.349, 4.699,
                                4.983, 5.100, 5.288, 5.374,
                                0.184, 0.664, 1.553, 1.910, 2.585,
                                3.009, 3.403, 3.892, 4.367, 4.551,
                                4.656, 4.754, 4.842, 4.969))


y <- algae$density
x <- algae$day

  # Reduced model is a quadratic polynomial:

XR <- cbind(1,x,x^2)

XtXinvR <- solve(t(XR)%*%XR)

betahatR <- XtXinvR%*%(t(XR)%*%y)

SSResR <- drop(t(y)%*%y - t(betahatR)%*%t(XR)%*%y)

print(SSResR)


  # Full model is a 13-degree polynomial:

XF <- cbind(1,poly(x,degree=13))

  # poly produces columns appropriate for polynomial regression of the
  # specified degree;  the polynomials are actually orthogonalized, but this
  # has no effect on SS(Res)

XtXinvF <- solve(t(XF)%*%XF)

betahatF <- XtXinvF%*%(t(XF)%*%y)

SSResF <- drop(t(y)%*%y - t(betahatF)%*%t(XF)%*%y)

print(SSResF)


  # To test for lack of fit ...

SSLF <- SSResR - SSResF

F0 <- (SSLF/(13-2))/(SSResF/(length(y)-14))

pval.F0 <- 1-pf(F0,13-2,length(y)-14)

print(SSLF)
print(F0)
print(pval.F0)


  # Another way to test lack of fit:

algae.red.mod <- lm(density ~ day + I(day^2), data=algae)
algae.full.mod <- lm(density ~ factor(day), data=algae)

anova(algae.red.mod, algae.full.mod)


  # Plot the data and the quadratic regression:

plot(algae$day, algae$density)
vals <- data.frame(day = seq(min(algae$day), max(algae$day), seq=100))
lines(vals$day, predict(algae.red.mod, vals))
title("Quadratic Regression Plot")

# Uncomment the following line to save the graph as a
# PostScript file in the current directory:

# dev.print(file="rpd08_03.ps")