# RPD -- Example 4.8 -- Physical Fitness Data

physfit <- read.table("physfit.dat", header=TRUE)
  # create a data frame from the variable names and data in "physfit.dat"

y <- cbind(physfit$oxygen)  # cbind supplies matrix attributes

X <- with(physfit, cbind(1,time,hr.rest,hr.run,hr.max))
  # use variables in physfit to form a matrix, columnwise

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

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

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

s2 <- SSResid/(dim(X)[1] - dim(X)[2])  # dim(A) = vector of dimensions of A

SSMu <- length(y) * mean(y)^2

SSTotalC <- sum(y^2) - SSMu

SSReg <- drop(t(betahat)%*%t(X)%*%y) - SSMu

MSReg <- SSReg/(dim(X)[2]-1)

print(betahat)
print(XtXinv)
print(SSTotalC)
print(SSReg)
print(SSResid)
print(MSReg)
print(s2)


  # To test whether hr.rest and hr.max can be dropped ...

K1 <- cbind(c(0,0,1,0,0), c(0,0,0,0,1))

m1 <- rbind(0,0)

Q1 <- drop(t(t(K1)%*%betahat-m1) %*% solve(t(K1)%*%XtXinv%*%K1) %*%
           (t(K1)%*%betahat-m1))

F1 <- (Q1/dim(K1)[2])/s2

pval1 <- 1-pf(F1,dim(K1)[2],dim(X)[1]-dim(X)[2])  # pf(f,n,d) = c.d.f. at f
                                                  # of F(n,d)-distribution

print(Q1)
print(F1)
print(pval1)


  # To simultaneously test whether the intercept is 90 and 
  # hr.rest and hr.max can be dropped ...

K2 <- cbind(c(1,0,0,0,0), c(0,0,1,0,0), c(0,0,0,0,1))

m2 <- rbind(90,0,0)

K2b.m <- t(K2)%*%betahat-m2

K2XtXiK2 <- t(K2)%*%XtXinv%*%K2

Q2 <- drop(t(K2b.m) %*% solve(K2XtXiK2) %*% K2b.m)

F2 <- (Q2/dim(K2)[2])/s2

pval2 <- 1-pf(F2,dim(K2)[2],dim(X)[1]-dim(X)[2])

print(K2b.m)
print(K2XtXiK2)
print(Q2)
print(F2)
print(pval2)



physfit.model <- lm(oxygen ~ time + hr.rest + hr.run + hr.max, data=physfit)

coef(physfit.model)  # betahat

chol2inv(physfit.model$qr$qr)  # inverse of t(X)%*%X

anova(physfit.model)["Residuals","Sum Sq"]  # SS(Res)

library(gmodels)  # a package containing function glh.test
                  # to test the general linear hypothesis
                  # (on older installations, use "library(gregmisc)" instead)
                  # may need to be installed using "install.packages()"

summary(glh.test(physfit.model, t(K1), m1))

summary(glh.test(physfit.model, t(K2), m2))

  # Note:  The library "multcomp" also provides a function "glht"
  # for testing the general linear hypothesis.