# RPD -- Example 8.2 -- Algae Density: CO2 only, Rep 1

algae <- data.frame(day = 1:14,
                    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))


plot(algae$day, algae$density)


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

  # Try a cubic polynomial fit:

X <- cbind(1,x,x^2,x^3)

XtXinv <- solve(t(X)%*%X)

betahat <- XtXinv%*%(t(X)%*%y)

SSRes <- drop(t(y)%*%y - t(betahat)%*%t(X)%*%y)

s2 <- SSRes/(dim(X)[1]-dim(X)[2])

se.betahat <- sqrt(diag(XtXinv*s2))

print(X)
print(betahat)
print(se.betahat)


  # To test whether a quadratic model is adequate ...

t0 <- betahat[4]/se.betahat[4]

pval.t0 <- 2*(1-pt(abs(t0),dim(X)[1]-dim(X)[2]))  # pt(t,d) = c.d.f. at t
                                                  # of t(d)-distribution

print(t0)
print(pval.t0)


  # To test whether a straight-line model is adequate ...

XR <- cbind(1,x)

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

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

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

F0 <- ((SSResR - SSRes)/2)/s2

pval.F0 <- 1-pf(F0,2,dim(X)[1]-dim(X)[2])

print(F0)
print(pval.F0)


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

  # function 'I' ensures that '^' is interpreted as the usual power
  # (exponent) operator, rather than as a formula expansion operator

summary(algae.full.mod)
anova(algae.full.mod)


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

anova(algae.red.mod1, algae.full.mod)  # quadratic adequate?


algae.red.mod2 <- lm(density ~ day, data=algae)

anova(algae.red.mod2, algae.full.mod)  # simple linear adequate?