# RPD -- Example 8.4 -- 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))


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

  # Cubic orthogonal polynomial fit:

Q <- cbind(1/sqrt(length(x)),poly(x,degree=3))

  # Note:  Q is not quite the same as the matrix in the
  # textbook because the columns of Q have length 1.

QtQinv <- solve(t(Q)%*%Q)  # should be the identity matrix

gammahat <- QtQinv%*%(t(Q)%*%y)

SSRes <- drop(t(y)%*%y - t(gammahat)%*%t(Q)%*%y)

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

print(Q)
print(zapsmall(QtQinv))  # zapsmall sets very small values to zero
print(gammahat)
print(s2)


  # To translate back to the usual polynomial coefficient estimates ...

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

R <- t(Q)%*%X

betahat <- solve(R,gammahat)

print(zapsmall(R))
print(betahat)


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

t0 <- gammahat[4]/sqrt(QtQinv[4,4]*s2)

pval.t0 <- 2*(1-pt(abs(t0),dim(Q)[1]-dim(Q)[2]))

print(t0)
print(pval.t0)