[R] lme model specification problem (Error in MEEM...)

[Ben Bolker]

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