[R] generalized inverse using matinv (Design)

[David Winsemius]

On Aug 16, 2011, at 10:15 AM, Frank Paetzold wrote:

> i am trying to use matinv from the Design package

Which has been replaced by the rms package (about 2 years ago).

> to compute the generalized inverse of the normal equations
> of a 3x3 design via the sweep operator.
>
> That is, for the linear model
> y = µ + x1 + x2 + x1*x2
> where x1, x2 are 3-level factors and dummy coding is being used
>
> the matrix to be inverted is
>
> X'X =
>
> 9 3 3 3 3 3 3 1 1 1 1 1 1 1 1 1
> 3 3 0 0 1 1 1 1 0 0 1 0 0 1 0 0
> 3 0 3 0 1 1 1 0 1 0 0 1 0 0 1 0
> 3 0 0 3 1 1 1 0 0 1 0 0 1 0 0 1
> 3 1 1 1 3 0 0 1 1 1 0 0 0 0 0 0
> 3 1 1 1 0 3 0 0 0 0 1 1 1 0 0 0
> 3 1 1 1 0 0 3 0 0 0 0 0 0 1 1 1
> 1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0
> 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0
> 1 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0
> 1 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0
> 1 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0
> 1 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0
> 1 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0
> 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0
> 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1

It would have been courteous to supply paste()-able code. If you use  
this code and then paste in the above matrix, then others would not  
need to do it themselves:

XtX = matrix(scan(), byrow=TRUE, ncol=16)

>
> this matrix has rank=9

The rankMatrix::Matrix function agrees. And in essence, so does the  
eigen function, since the last 7 eigenvalues are effectively zero:

 > eigen(XtX)
$values
  [1]  1.600000e+01  4.000000e+00  4.000000e+00  4.000000e+00   
4.000000e+00
  [6]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00   
1.065814e-14
[11]  4.440892e-15  4.173627e-16  2.220446e-16  1.717376e-16  
-1.193166e-16
[16] -9.593521e-16


>
> however, matinv(X'X) falsely returns rank=4,

Well, the current version of matinv returns a different value (= 6) ,  
although it is not 9.

> no matter what the tolerance threshold eps is set to.
>
> also the defining property of the
> generalized inverse
>             _
> X'X %*% (X'X)  %*% X'X = X'X

See the FAQ regarding issues related to "==" and numerical accuracy.

>
> is not satisfied.
>
> if i use qr (from the base package) the rank is
> correctly determined as 9.

Since matinv is not described as being designed to deal with rank- 
deficient matrices, why have you chosen it over a function designed  
for the sort of problem you are dealing with. You remember the old  
joke about the guy who goes to the doctor and complains: "Doc, it  
hurts a lot when I do ....".

I guess the more modern Zen-oriented question might be: "And how has  
that method been working out for you?"

??"generalized inverse"

I see two "generalized inverse"  functions in the packages I have  
installed: ginv::MASS and pinv::pracma

require(MASS)
 > all.equal(XtX %*% ginv(XtX) %*% XtX ,  XtX)
[1] TRUE


-- 
David Winsemius, MD
West Hartford, CT