# RPD -- Example 4.10 -- Physical Fitness Data

physfit <- read.table("physfit.dat", header=TRUE)

y <- cbind(physfit$oxygen)

X0 <- matrix(1, length(physfit$oxygen), 1)  # column vector of 1s
X1 <- cbind(physfit$time)
X2 <- cbind(physfit$hr.rest)
X3 <- cbind(physfit$hr.run)
X4 <- cbind(physfit$hr.max)

R0 <- drop(t(y)%*%X0%*%solve(t(X0)%*%X0)%*%t(X0)%*%y)  # SS(mu)

X01 <- cbind(X0,X1)

R01 <- drop(t(y)%*%X01%*%solve(t(X01)%*%X01)%*%t(X01)%*%y)

X013 <- cbind(X0,X1,X3)

R013 <- drop(t(y)%*%X013%*%solve(t(X013)%*%X013)%*%t(X013)%*%y)

X0132 <- cbind(X0,X1,X3,X2)

R0132 <- drop(t(y)%*%X0132%*%solve(t(X0132)%*%X0132)%*%t(X0132)%*%y)

X01324 <- cbind(X0,X1,X3,X2,X4)

R01324 <- drop(t(y)%*%X01324%*%solve(t(X01324)%*%X01324)%*%t(X01324)%*%y)

X0124 <- cbind(X0,X1,X2,X4)

R0124 <- drop(t(y)%*%X0124%*%solve(t(X0124)%*%X0124)%*%t(X0124)%*%y)

X0134 <- cbind(X0,X1,X3,X4)

R0134 <- drop(t(y)%*%X0134%*%solve(t(X0134)%*%X0134)%*%t(X0134)%*%y)

X0324 <- cbind(X0,X3,X2,X4)

R0324 <- drop(t(y)%*%X0324%*%solve(t(X0324)%*%X0324)%*%t(X0324)%*%y)

R1324.0 <- R01324 - R0  # SS(Regr)


  # Sequential Sums of Squares in order 1,3,2,4

R1.0 <- R01 - R0

R3.01 <- R013 - R01

R2.013 <- R0132 - R013

R4.0132 <- R01324 - R0132

SSResid <- drop(t(y)%*%y - t(y)%*%X01324%*%solve(t(X01324)%*%X01324)
                %*%t(X01324)%*%y)

print(R1324.0)
print(R1.0)
print(R3.01)
print(R2.013)
print(R4.0132)
print(SSResid)


  # Partial Sums of Squares

R1.0324 <- R01324 - R0324

R3.0124 <- R01324 - R0124

R2.0134 <- R01324 - R0134

print(R1.0324)
print(R3.0124)
print(R2.0134)
print(R4.0132)

  # Note:  Practical computer algorithms for regression use much 
  # more efficient methods for computing these sums of squares.


physfit.model <- lm(oxygen ~ time + hr.run + hr.rest + hr.max, data=physfit)
  # Note the order in which the independent variables are listed!

anova(physfit.model)  # gives sequential sums of squares

drop1(physfit.model)  # gives partial sums of squares