# RPD -- Example 9.2 -- Numerical Example

artificial <- data.frame(group = gl(4,2),
                         y = c( 8.9,  8.76, 11.78, 12.07,
                               14.5, 12.48, 16.79, 16.57))

  # gl(n,k) = vector of factor levels 1 to n, repeated k times each (balanced)


# Cell Means Model:

X1 <- model.matrix(~ 0 + group, data=artificial)

  # matrix for cell means parameterization (one factor, no intercept)

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

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

print(X1)
print(XtXinv1)
print(betahat1)


# Effects Model: Sum-to-Zero

X2 <- model.matrix(~ group, data=artificial,
                   contrasts.arg = list(group="contr.sum"))

  # matrix for sum-to-zero parameterization

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

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

tauhat4 <- - betahat2[2] - betahat2[3] - betahat2[4]

print(X2)
print(XtXinv2)
print(betahat2)
print(tauhat4)


# Effects Model: Reference Cell

X3 <- model.matrix(~ group, data=artificial,
                   contrasts.arg = list(group="contr.SAS"))

  # matrix for reference-cell parameterization (using last level, as in SAS(R))

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

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

print(X3)
print(XtXinv3)
print(betahat3)


artificial.model1 <- lm(y ~ group, data=artificial)
summary(artificial.model1)

  # gives estimates of a type defined by 'options("contrasts")';
  # typically, this is "contr.treatment", which gives "reference cell model"
  # estimates, with the first (NOT the last) level as the reference cell


artificial.model2 <- lm(y ~ C(group,contr.SAS), data=artificial)
summary(artificial.model2)

  # gives "reference cell model" estimates, with the last level
  # as the reference cell, as in SAS(R)


artificial.model3 <- lm(y ~ 0 + group, data=artificial)
summary(artificial.model3)

  # gives "cell means model" estimates;
  # Note: The listed ANOVA F-test is for H0: all means are ZERO,
  # and the R-Squared value is wrong because it doesn't account for the
  # implicit intercept.


artificial.model4 <- lm(y ~ C(group,contr.sum), data=artificial)
summary(artificial.model4)

  # gives "effects model" estimates (except for the last level)
  # that satisfy the sum-to-zero restriction