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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
