R help - Aug 2011 - How set lm() to don't return NA in summary()?

Hi,

I've data from an incomplete fatorial design. One level of a factor
doesn't
has the levels of the other. When I use lm(), the summary() return NA for
that non estimable parameters. Ok, I understant it. But I use
contrast::contrast(), gmodels::estimable(), multcomp::glht() and all these
fail when model has NA estimates. This is becouse vcov() and coef() has
different dimensions. Is possible set lm() to omit NA? Below same toy data
and code.
> # toy data
> adi <- expand.grid(cult=gl(1,3,la=LETTERS[1]), fert=101)
> fat <- expand.grid(cult=gl(2,3,la=LETTERS[2:3]), fert=seq(50,150,50))
> da <- rbind(adi, fat)
> da$y <- da$fert+rnorm(nrow(da),0,10)
>
> # plot
> require(lattice)
> xyplot(y~fert|cult, da)
>
> # adjust
> m0 <- lm(y~cult*fert, da)
> summary(m0)
...
Coefficients: (1 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  7.55401   10.18956   0.741    0.469
cultB       -8.04486   13.54672  -0.594    0.561
cultC       -3.83644    6.74848  -0.568    0.578
fert         0.87486    0.08265  10.586 1.24e-08 ***
cultB:fert   0.13509    0.11688   1.156    0.265
cultC:fert        NA         NA      NA       NA
...
> require(gmodels)
> estimable(m0, cm=c(1,0,0,100,0,0))
Erro em estimable.default(m0, cm = c(1, 0, 0, 100, 0, 0)) :
 Dimension of structure(c(1, 0, 0, 100, 0, 0), .Dim = c(1L, 6L)): 1x6, not
compatible with no of parameters in m0: 5

Thanks.
Walmes.

=========================================================================Walmes
Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: walmes@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
=========================================================================
	[[alternative HTML version deleted]]

On Aug 6, 2011, at 03:23 , Walmes Zeviani wrote:
> Hi,
> 
> I've data from an incomplete fatorial design. One level of a factor
doesn't
> has the levels of the other. When I use lm(), the summary() return NA for
> that non estimable parameters. Ok, I understant it. But I use
> contrast::contrast(), gmodels::estimable(), multcomp::glht() and all these
> fail when model has NA estimates. This is becouse vcov() and coef() has
> different dimensions. Is possible set lm() to omit NA? Below same toy data
> and code.
> 
>> # toy data
>> adi <- expand.grid(cult=gl(1,3,la=LETTERS[1]), fert=101)
>> fat <- expand.grid(cult=gl(2,3,la=LETTERS[2:3]),
fert=seq(50,150,50))
>> da <- rbind(adi, fat)
>> da$y <- da$fert+rnorm(nrow(da),0,10)
>> 
>> # plot
>> require(lattice)
>> xyplot(y~fert|cult, da)
>> 
>> # adjust
>> m0 <- lm(y~cult*fert, da)
>> summary(m0)
> 
> ...
> Coefficients: (1 not defined because of singularities)
>            Estimate Std. Error t value Pr(>|t|)
> (Intercept)  7.55401   10.18956   0.741    0.469
> cultB       -8.04486   13.54672  -0.594    0.561
> cultC       -3.83644    6.74848  -0.568    0.578
> fert         0.87486    0.08265  10.586 1.24e-08 ***
> cultB:fert   0.13509    0.11688   1.156    0.265
> cultC:fert        NA         NA      NA       NA
> ...
>> require(gmodels)
>> estimable(m0, cm=c(1,0,0,100,0,0))
> Erro em estimable.default(m0, cm = c(1, 0, 0, 100, 0, 0)) :
> Dimension of structure(c(1, 0, 0, 100, 0, 0), .Dim = c(1L, 6L)): 1x6, not
> compatible with no of parameters in m0: 5
...and of course (you forgot to say) just removing the 6th element of cm
won't do either.

A generic trick is to modify the model matrix directly:
> M <- model.matrix(m0)[,-6]
> m1 <- update(m0, ~ M - 1)
> summary(m1)
Call:
lm(formula = y ~ M - 1, data = da)

Residuals:
    Min      1Q  Median      3Q     Max 
-27.165  -4.704  -2.446   6.666  25.328 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
M(Intercept)  -1.8618    14.6718  -0.127    0.901    
McultB         2.2102    19.5058   0.113    0.911    
McultC        11.5308     9.7171   1.187    0.253    
Mfert          0.9506     0.1190   7.989 5.65e-07 ***
McultB:fert    0.0758     0.1683   0.450    0.658    
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 

Residual standard error: 14.57 on 16 degrees of freedom
Multiple R-squared: 0.9866,	Adjusted R-squared: 0.9824 
F-statistic: 235.5 on 5 and 16 DF,  p-value: 2.165e-14 
> estimable(m1, cm=c(1,0,0,100,0))
              Estimate Std. Error  t value DF     Pr(>|t|)
(1 0 0 100 0) 93.20357   8.415451 11.07529 16 6.515116e-09


Two things to notice: The coefficient names are not quite the same as the
original, and you need to fit without the intercept (or remove the 1st column as
well).

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com
"D?den skal tape!" --- Nordahl Grieg