R help - Jan 2012 - lme model specification problem (Error in MEEM...)

Dear all,

In lme, models in which a factor is fully "contained" in another lead
to
an error. This is not the case when using lm/aov.

I understand that these factors are aliased, but believe that such 
models make sense when the factors are fitted sequentially. For example, 
I sometimes fit a factor first as linear term (continuous variable with 
discrete levels, e.g. 1,2,4,6), and then as factor, with the latter then 
accounting for the deviation from linearity.

If x contains the values 1,2,4,6, the model formula then would be y ~ x 
+ as.factor(x)

A complete example is here:

library(nlme)

d <- data.frame(x=rep(c(1,2,4,6),each=2),subj=factor(rep(1:16,each=2)))
d$y <- rnorm(32) + rnorm(length(levels(d$subj)))[d$subj]

anova(lme(y~x,random=~1|subj,data=d))               ## works
anova(lme(y~as.factor(x),random=~1|subj,data=d))    ## works

anova(lme(y~x+as.factor(x),random=~1|subj,data=d))  ## fails:
# Error in MEEM(object, conLin, control$niterEM) :
#  Singularity in backsolve at level 0, block 1

summary(aov(y~x+as.factor(x)+Error(subj),data=d))   ## works:

# Error: subj
#              Df Sum Sq Mean Sq F value Pr(>F)
# x             1  8.434   8.434   4.780 0.0493 *
# as.factor(x)  2 10.459   5.230   2.964 0.0900 .
# Residuals    12 21.176   1.765
# ... rest of output removed ...

I understand I can to some extent work around this limitation by 
modifying the contrast encoding, but I then still don't get an overall 
test for "as.factor(x)" (in the example above, a test for deviation
from
linearity).

d$xfac <- factor(d$x)
contrasts(d$xfac)<-c(1,2,4,6)
summary(lme(y~xfac,random=~1|subj,data=d))

Is there a way to work around this limitation of lme? Or do I 
mis-specify the model? Thanks for your advice.

Pascal

Pascal A. Niklaus <pascal.niklaus <at> ieu.uzh.ch> writes:
> In lme, models in which a factor is fully "contained" in another
lead to
> an error. This is not the case when using lm/aov.
> 
> I understand that these factors are aliased, but believe that such 
> models make sense when the factors are fitted sequentially. For example, 
> I sometimes fit a factor first as linear term (continuous variable with 
> discrete levels, e.g. 1,2,4,6), and then as factor, with the latter then 
> accounting for the deviation from linearity.
> 
> If x contains the values 1,2,4,6, the model formula then would be y ~ x 
> + as.factor(x)
> 
> A complete example is here:
> 
> library(nlme)
> 
> d <- data.frame(x=rep(c(1,2,4,6),each=2),subj=factor(rep(1:16,each=2)))
> d$y <- rnorm(32) + rnorm(length(levels(d$subj)))[d$subj]
> 
> anova(lme(y~x,random=~1|subj,data=d))               ## works
> anova(lme(y~as.factor(x),random=~1|subj,data=d))    ## works
> 
> anova(lme(y~x+as.factor(x),random=~1|subj,data=d))  ## fails:
> # Error in MEEM(object, conLin, control$niterEM) :
> #  Singularity in backsolve at level 0, block 1
> 
> summary(aov(y~x+as.factor(x)+Error(subj),data=d))   ## works:
> 
> # Error: subj
> #              Df Sum Sq Mean Sq F value Pr(>F)
> # x             1  8.434   8.434   4.780 0.0493 *
> # as.factor(x)  2 10.459   5.230   2.964 0.0900 .
> # Residuals    12 21.176   1.765
> # ... rest of output removed ...
> 
> I understand I can to some extent work around this limitation by 
> modifying the contrast encoding, but I then still don't get an overall 
> test for "as.factor(x)" (in the example above, a test for
deviation from
> linearity).
> 
> d$xfac <- factor(d$x)
> contrasts(d$xfac)<-c(1,2,4,6)
> summary(lme(y~xfac,random=~1|subj,data=d))
> 
> Is there a way to work around this limitation of lme? Or do I 
> mis-specify the model? Thanks for your advice.
  This question would probably be more appropriate for 
the r-sig-mixed-models at r-project.org mailing list.

  A short answer (if you re-post to r-sig-mixed-models I might
answer at greater length) would be that you should be able to quantify
the difference between

  y~x

 and

  y~factor(x)

by an anova comparison of the two *models*, i.e.

anova(m1,m2)

Another thing to consider would be setting using orthogonal
polynomial contrasts for factor(x), where summary() would
give you a result for the constant, linear, quadratic ...
terms

  Ben Bolker