# RPD -- Example 9.4 -- Listening-Reading Skills Test

listen <- data.frame(trt = gl(3,4),
                     pretest = c(89,82,88,94,89,90,91,92,89,99,84,87),
                     posttest = c(87,86,94,96,84,94,97,93,96,97,100,98))


y <- listen$posttest

  # Full model --- slopes and intercepts differ by trt:

X1 <- model.matrix(~ 0 + trt + trt:pretest, data=listen)

  # in a formula, the operator : produces a pure interaction

print(X1)

XtXinv1 <- solve(t(X1)%*%X1)

betahat1 <- XtXinv1%*%(t(X1)%*%y)

SSRes1 <- drop(t(y)%*%y - t(betahat1)%*%t(X1)%*%y)

s2.1 <- SSRes1/(dim(X1)[1]-dim(X1)[2])

print(SSRes1)


  # Reduced model with same slopes, different intercepts:

X2 <- model.matrix(~ 0 + trt + pretest, data=listen)

print(X2)

XtXinv2 <- solve(t(X2)%*%X2)

betahat2 <- XtXinv2%*%(t(X2)%*%y)

SSRes2 <- drop(t(y)%*%y - t(betahat2)%*%t(X2)%*%y)

s2.2 <- SSRes2/(dim(X2)[1]-dim(X2)[2])

print(SSRes2)


  # Test for equal slopes:

F0 <- (SSRes2 - SSRes1)/(dim(X1)[2] - dim(X2)[2]) / s2.1

pvalF0 <- 1-pf(F0,dim(X1)[2]-dim(X2)[2],dim(X1)[1]-dim(X1)[2])

print(F0)
print(pvalF0)


  # Further reduced model with same slopes and intercepts:

X3 <-  model.matrix(~ pretest, data=listen)

print(X3)

XtXinv3 <- solve(t(X3)%*%X3)

betahat3 <- XtXinv3%*%(t(X3)%*%y)

SSRes3 <- drop(t(y)%*%y - t(betahat3)%*%t(X3)%*%y)

print(SSRes3)


  # Test for equal intercepts, assuming equal slopes:

F0 <- (SSRes3 - SSRes2)/(dim(X2)[2] - dim(X3)[2]) / s2.2

  # using s^2 from the equal-slopes model as the denominator, but could
  # alternatively use s^2 from the full model (and change df below)

pvalF0 <- 1-pf(F0,dim(X2)[2]-dim(X3)[2],dim(X2)[1]-dim(X2)[2])

print(F0)
print(pvalF0)


  # Joint test for equal intercepts and equal slopes:

F0 <- (SSRes3 - SSRes1)/(dim(X1)[2] - dim(X3)[2]) / s2.1

pvalF0 <- 1-pf(F0,dim(X1)[2]-dim(X3)[2],dim(X1)[1]-dim(X1)[2])

print(F0)
print(pvalF0)



listen.model1 <- lm(posttest ~ trt*pretest, data=listen)
anova(listen.model1)

  # in a formula, the operator * produces both main effects and interaction


listen.model2 <- lm(posttest ~ trt + pretest, data=listen)
anova(listen.model2)


anova(listen.model2, listen.model1)  # equal slopes?


listen.model3 <- lm(posttest ~ pretest, data=listen)


anova(listen.model3, listen.model1)  # equal slopes and intercepts?

anova(listen.model3, listen.model2)  # equal intercepts (given equal slopes)?