Stat 505 Assignment 3

1. Find a 4 by 6 matrix A and a particular generalized inverse: A^g so that A A^g = I₄, but no element of A is one.
2. Find a 4 x 6 matrix B and B^g so that B^g B=I₆ or prove none exist.
3. Use the warpbreaks data set, included with R. Use this model:
   y _ijk = μ + α _i + τ _j + γ _ij + ε _ijk
   Use these plotting functions or something better.

   with(warpbreaks, interaction.plot(tension, wool, breaks))
   plot(breaks ~ interaction( tension,wool), data = warpbreaks, col = "gray")
   

   1. Write up a report of the effects of tension and wool type (A or B) on the number of breaks. Include the following points:
      If you could choose one combination of tension and wool type to minimize breaks, which would it be?
      Is your favorite wool type always better?
      Is your favorite tension tension always best?
      Do the residuals give any indication of a problem with the assumptions of ordinary least squares regression? If so, fix it.
      Fit the full interaction model and interpret each line of the table printed by the anova command in R.
   2. What size is the full model matrix? Write out the X matrix for this model showing each unique row.
   3. What size is X^T X? Use R to create it. Start with:

       myX <- with(warpbreaks, cbind(1, wool=="A",wool=="B", tension=="L", ...)
      

      Interaction columns really are elementwise products of main effects columns. Don't show the whole X, but do print out X^T X. Find its rank.

      (XtX <- crossprod(myX))
      qr(XtX)$rank
      

   4. Look at the model matrix R uses (apply the model.matrix function to your full fitted model) and its crossproduct. Verify or refute the claim: by deleting certain rows and the same columns of the full X^T X we can obtain the one R gives. Verify or refute: the R crossproduct model matrix is of full column rank. Use the ginv function in the MASS package to find a generalized inverse of X^TX.
   5. Create a zero matrix the same size as full X^T X. In the rows and columns which are not removed to get R's smaller matrix, set in the values of the generalized inverse above.

        RsXtXinv <- MASS::ginv( crossprod(model.matrix(my.fitted.model)))
        myXtXinv <- matrix(0,12,12)                     ## or ## 0 * myXtX
        omit.rows <- c(....)                    ## fill this in from above
        myXtXinv[-omit.rows, -omit.rows] <-  RsXtXinv
      

      Demonstrate (to machine precision) that the resulting matrix is a generalized inverse of X^T X. Note: in HW1 we used the zapsmall function to convert things very close to zero to zero. You'll want to subtract the answer you should get, and use zapsmall on the differences to demonstrating equality. Alternatively, you could do it by hand to get exact answers (factor out any 9's or 1/9ths).
   6. Compute X^T y and demonstrate that premultiplying by your generalized inverse above gives the same solutions that R spits out, with some estimates (which ones?) set to 0.
   7. Build a table with a row for each estimated coefficient from R. In terms of the greek letters in the model, what is that coefficient estimating?
   8. This is getting long, but we need to see that answers are not unique. Set these options and refit your model:

       options(contrasts=c("contr.SAS","contr.poly"))
      

      Again, find a ginv for R's crossproduct model matrix, and substitute its values into a matrix of 0's to get a different (confirm that) generalized inverse of X^TX (confirm it is a ginverse). Compare the coefficient estimates you get from multiplication of (X^TX)^g X^T y with the R estimates for this new model.
   9. Again, build a table with a row for each estimated coefficient from R. In terms of the greek letters in the model, what is that coefficient estimating? These should be different from the ones above.
   10. For each estimated coefficient vector, demonstrate that we get the same fitted values, Xb (You found two b's using the same full X).