[R] Wrong contrast matrix for nested factors in lm(), rlm(), and lmRob()

[Saul Sternberg]

This message also reports wrong estimates produced by lmRob.fit.compute()
for nested factors when using the correct contrast matrix.

And in these respects, I have found that S-Plus behaves the same way as R.

Using the three available contrast types (sum, treatment, helmert)
with lm() or lm.fit(), but just contr.sum with rlm() and lmRob(),
and small examples, I generated contrast matrices for four models
involving nested factors with fixed effects.  For three of the models
the matrices were incorrect.

-----------------------------------------------------------------
Details
---------
For lm() and rlm() I used two data frames:

In "same.df" the nested factor, D, has the same number of levels for
each level of the nesting factor, G:
   G  D  T1
1 g1 d1 -60
2 g1 d2 -50
3 g1 d3 -40
4 g2 d1  30
5 g2 d2  50
6 g2 d3  70

In "diff.df" the nested factor, D, has a different number of levels for
the two levels of the nesting factor, G:
   G  D  T2
1 g1 d1 -60
2 g1 d2 -50
3 g1 d3 -40
4 g2 d1  20
5 g2 d2  40
6 g2 d3  60
7 g2 d4  80

(G, D are factors; T1, T2 are numeric)

For lmRob() I expanded the data frames to 600 or 700 rows by replicating
them 100 times and adding error to the observations.

For lm() the models and commands were
(1) model.matrix(lm(T1 ~ G + D%in%G, same.df))
(2) model.matrix(lm(T1 ~ D%in%G, same.df))
(3) model.matrix(lm(T2 ~ G + D%in%G, diff.df))
(4) model.matrix(lm(T2 ~ D%in%G, diff.df))

Using (1), all three types of contrast matrix were correctly generated
Using (2), the same incorrect contrast matrix was generated for all
        three contrast types.
Using (3), an incorrect contrast matrix was generated for each of the
        three contrast types.  For contr.treatment the error was an extra
        column of zeros; for the others the error was more serious.
Using (4), the same incorrect contrast matrix was generated for all
        three contrast types.

I used only contr.sum with rlm() and lmRob().  For model (1) the programs
worked correctly, but for models (2), (3), and (4) with the formulas
above, rlm() and lmRob() both reported that x was singular.

When x was the correct contrast matrix and T was the observation vector,
rlm(x,T) worked correctly for (2), (3), and (4), as did lm.fit(x,T).
However, whereas lmRob.fit.compute(x2=NULL,y=T,x1=x) worked correctly for
(3), the estimates it produced for (2) and (4) were radically wrong
(and were the same for different random seeds and initial algorithms).

-------------------------------------------------------------------
Questions:
-------------

(1) If there is a way to use lm(), rlm(), and lmRob() in such cases
so that they generate the correct contrast matrices (and the desired
parameter estimates), what is it?

(2) If there is no way to do this, is the best alternative for the user to
create the desired model matrices "by hand" and provide them as arguments
to lim.fit(), rlm(), and lmRob.fit.compute()?  This would also require
that lmRob.fit.compute() generate the correct estimates.

(3) If one uses lm.fit() and lmRob.fit.compute() directly in this way,
then, given that one is warned against doing so, what are the dangers?

(4) According to cran.r-project.org/web/views/Robust.html, lmRob()
"makes use of the M-S algorithm of Maronna and Yohai (2000),
automatically when there are factors among the predictors (where
S-estimators (and hence MM-estimators) based on resampling typically
badly fail)."  Is there an alternative program that uses the M-S
algorithm, if lmRob() or lmRob.fit.compute() cannot be made to work? 

----------------------------------------------------------------
R%%>sessionInfo()
R version 2.12.0 (2010-10-15)
Platform: x86_64-redhat-linux-gnu (64-bit)
----------------------------------------------------------------
I'll be very grateful for any help.
Saul Sternberg, Psychology
University of Pennsylvania