>  # balanced.r                 
>  #  jfm  September 07, revised September 08
>  #  show some relationships among projection matrices
>  #   in balanced two-way crossed models  
>  #   Exercise 3.24 of 'A Primer on Linear Models' 
>  #
>   a <- 3            # number of levels of factor A (alpha's)
>   b <- 4            # number of levels of factor B (beta's)
>   N <- a*b          # number of observations
>   One <- matrix(1,N,1)        # intercept column
>   XA  <- matrix(c(rep(1,b),rep(0,N)),N,a)        # factor A cols
>   XB  <- matrix(rep(diag(1,b),a),N,b,byrow=TRUE) # factor B cols
>   X <- cbind(One,XA,XB)
>  #                    write out design matrix
>   X
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
 [1,]    1    1    0    0    1    0    0    0
 [2,]    1    1    0    0    0    1    0    0
 [3,]    1    1    0    0    0    0    1    0
 [4,]    1    1    0    0    0    0    0    1
 [5,]    1    0    1    0    1    0    0    0
 [6,]    1    0    1    0    0    1    0    0
 [7,]    1    0    1    0    0    0    1    0
 [8,]    1    0    1    0    0    0    0    1
 [9,]    1    0    0    1    1    0    0    0
[10,]    1    0    0    1    0    1    0    0
[11,]    1    0    0    1    0    0    1    0
[12,]    1    0    0    1    0    0    0    1
>  #                    do QR factorizations
>   qrOne <- qr(One)
>   qrXA <- qr(XA)
>   qrXB <- qr(XB)
>   qrX <- qr(X)
>   qrX$rank                    # rank of X
[1] 6
>  #                    the questions
>  #  a) show Pxa * XB = Pone * XB
>   ZeroS <- matrix(0,N,b)     # start with all zeros
>   Q1XB <- qr.qty(qrOne,XB)
>   ZeroS[1,] <- Q1XB[1,]      # rows correspond to rank of One
>   qr.qy(qrOne,ZeroS)         # Pone * XB 
      [,1] [,2] [,3] [,4]
 [1,] 0.25 0.25 0.25 0.25
 [2,] 0.25 0.25 0.25 0.25
 [3,] 0.25 0.25 0.25 0.25
 [4,] 0.25 0.25 0.25 0.25
 [5,] 0.25 0.25 0.25 0.25
 [6,] 0.25 0.25 0.25 0.25
 [7,] 0.25 0.25 0.25 0.25
 [8,] 0.25 0.25 0.25 0.25
 [9,] 0.25 0.25 0.25 0.25
[10,] 0.25 0.25 0.25 0.25
[11,] 0.25 0.25 0.25 0.25
[12,] 0.25 0.25 0.25 0.25
>   ZeroS <- matrix(0,N,b)     # start with all zeros
>   Q1XB <- qr.qty(qrXA,XB)
>   ZeroS[1:a,] <- Q1XB[1:a,]  # rows correspond to rank of XA
>   qr.qy(qrXA,ZeroS)          # Pxa * XB 
      [,1] [,2] [,3] [,4]
 [1,] 0.25 0.25 0.25 0.25
 [2,] 0.25 0.25 0.25 0.25
 [3,] 0.25 0.25 0.25 0.25
 [4,] 0.25 0.25 0.25 0.25
 [5,] 0.25 0.25 0.25 0.25
 [6,] 0.25 0.25 0.25 0.25
 [7,] 0.25 0.25 0.25 0.25
 [8,] 0.25 0.25 0.25 0.25
 [9,] 0.25 0.25 0.25 0.25
[10,] 0.25 0.25 0.25 0.25
[11,] 0.25 0.25 0.25 0.25
[12,] 0.25 0.25 0.25 0.25
>  #  b) show Pxa*Pxb = Pone
>   QA <-qr.Q(qrXA)
>   QB <-qr.Q(qrXB)
>   inner <- t(QA) %*% QB      # QB'QA -- product in the middle
>  # Pxa * Pxb, is it all 1/N ?   (multiply & just print first four rows) 
>   N*( QA %*% inner %*% t(QB) )[1:4,]   
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,]    1    1    1    1    1    1    1    1    1     1     1     1
[2,]    1    1    1    1    1    1    1    1    1     1     1     1
[3,]    1    1    1    1    1    1    1    1    1     1     1     1
[4,]    1    1    1    1    1    1    1    1    1     1     1     1
>  #  c) show Px - Pxa = Pxb - Pone  or Px - Pxa - Pxb = - Pone
>   QX <- qr.Q(qrX)[,1:qrX$rank]     # watch rank!
>  #                             again multiply & look for ones
>   (N*( QX %*% t(QX) - QA %*% t(QA) - QB %*% t(QB) ))[1:4,]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,]   -1   -1   -1   -1   -1   -1   -1   -1   -1    -1    -1    -1
[2,]   -1   -1   -1   -1   -1   -1   -1   -1   -1    -1    -1    -1
[3,]   -1   -1   -1   -1   -1   -1   -1   -1   -1    -1    -1    -1
[4,]   -1   -1   -1   -1   -1   -1   -1   -1   -1    -1    -1    -1
>  #  see how rank-deficient QR works
>   qrX
$qr
            [,1]       [,2]          [,3]          [,4]          [,5]
 [1,] -3.4641016 -1.1547005 -1.154701e+00 -8.660254e-01 -8.660254e-01
 [2,]  0.2886751 -1.6329932  8.164966e-01 -2.326698e-16 -1.152507e-16
 [3,]  0.2886751  0.3167969  1.414214e+00  1.665335e-16 -5.551115e-17
 [4,]  0.2886751  0.3167969  9.335840e-17  1.500000e+00 -5.000000e-01
 [5,]  0.2886751 -0.2955755 -3.535534e-01 -2.506214e-01  1.414214e+00
 [6,]  0.2886751 -0.2955755 -3.535534e-01  4.160452e-01 -5.736659e-01
 [7,]  0.2886751 -0.2955755 -3.535534e-01  4.160452e-01  1.334409e-01
 [8,]  0.2886751 -0.2955755 -3.535534e-01  4.160452e-01  1.334409e-01
 [9,]  0.2886751 -0.2955755  3.535534e-01 -4.002651e-01 -1.139136e-01
[10,]  0.2886751 -0.2955755  3.535534e-01  2.664015e-01 -7.774087e-01
[11,]  0.2886751 -0.2955755  3.535534e-01  2.664015e-01 -7.030187e-02
[12,]  0.2886751 -0.2955755  3.535534e-01  2.664015e-01 -7.030187e-02
               [,6]          [,7]          [,8]
 [1,] -8.660254e-01 -1.154701e+00 -8.660254e-01
 [2,] -3.790642e-17  8.164966e-01 -3.790642e-17
 [3,]  0.000000e+00 -1.414214e+00 -8.326673e-17
 [4,] -5.000000e-01  4.205918e-16 -5.000000e-01
 [5,] -7.071068e-01  9.681388e-17 -7.071068e-01
 [6,] -1.224745e+00  8.521513e-17  1.224745e+00
 [7,]  6.854355e-01 -6.221446e-16 -6.087617e-16
 [8,] -1.310610e-01 -1.615489e-01 -4.244480e-16
 [9,] -2.729434e-01 -6.728298e-01 -2.990435e-01
[10,] -1.185242e-02 -4.838408e-01  3.818258e-03
[11,]  3.575361e-01  3.714062e-01 -1.269873e-01
[12,] -4.589605e-01 -1.471997e-01  7.046753e-01

$rank
[1] 6

$qraux
[1] 1.288675 1.316797 1.000000 1.211628 1.089829 1.316047 1.357057 1.630766

$pivot
[1] 1 2 3 5 6 7 4 8

attr(,"class")
[1] "qr"
>  #  done
>   rm(list=ls())
>  q()