[Rd] model.matrix() may be misleading for "lme" models

[John Fox]

Dear r-devel list members,

I'm posting this message here because it concerns the nlme package, 
which is maintained by R-core. The problem I'm about to describe is 
somewhere between a bug and a feature request, and so I thought it a 
good idea to ask here rather posting a bug report to the R bugzilla.

I was made aware (by Ben Bolker) that the car::Anova() method for "lme" 
models reports unreasonable results and warnings for mixed models in 
which non-default contrasts are used for factors. I traced the problem 
to a call to model.matrix() on "lme" objects, which was introduced into 
car:::Anova.lme() a couple of years ago to check for problems in the 
model matrix. That invokes model.matrix.default(), which ignores the 
contrasts defined in a call to lme().

Here's a simple direct example:

------- snip -------

 > library(nlme)
 > m <- lme(distance ~ Sex, random = ~ 1 | Subject,
+          data=Orthodont, contrasts=list(Sex = "contr.sum"))
 > m

Linear mixed-effects model fit by REML
   Data: Orthodont
   Log-restricted-likelihood: -253.629
   Fixed: distance ~ Sex
(Intercept)        Sex1
   23.808239    1.160511

Random effects:
  Formula: ~1 | Subject
         (Intercept) Residual
StdDev:    1.595838 2.220312

Number of Observations: 108
Number of Groups: 27

 > model.matrix(m, data=Orthodont)

     (Intercept) SexFemale
1             1         0
2             1         0
3             1         0

. . .

106           1         1
107           1         1
108           1         1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$Sex
[1] "contr.treatment"

--------- snip ---------

So model.matrix() constructs the model matrix using contr.treatment() 
even though contr.sum() was used by lme() to fit the model.

In contrast (pun intended), model.matrix() works as expected with an 
"lm" model (via model.matrix.lm()):

--------- snip ---------

 > m.lm <- lm(distance ~ Sex, data=Orthodont,
+            contrasts=list(Sex = "contr.sum"))
 > m.lm

Call:
lm(formula = distance ~ Sex, data = Orthodont,
    contrasts = list(Sex = "contr.sum"))

Coefficients:
(Intercept)         Sex1
      23.808        1.161

 > model.matrix(m.lm)
     (Intercept) Sex1
1             1    1
2             1    1
3             1    1

. . .

106           1   -1
107           1   -1
108           1   -1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$Sex
[1] "contr.sum"

--------- snip ---------

I was able to get around this problem by defining a model.matrix.lme() 
method, which is used internally in the car package but not registered:

--------- snip ---------

model.matrix.lme <- function(object, ...){

   data <- object$data
   contrasts <- object$contrasts

   if (length(contrasts) == 0) {
     xlev <- NULL
   } else {
     xlev <- vector(length(contrasts), mode="list")
     names(xlev) <- names <- names(contrasts)
     for (name in names){
       xlev[[name]] <- rownames(contrasts[[name]])
     }
   }

   if (is.null(data)){
     NextMethod(formula(object), data=eval(object$call$data),
                  contrasts.arg=contrasts, xlev=xlev, ...)
   } else {
     NextMethod(formula(object), data=data,
                     contrasts.arg=contrasts, xlev=xlev, ...)
   }
}

--------- snip ---------

This function is a bit awkward, particularly the part that constructs 
the xlev argument, but it does appear to work as intended (note, 
however, that the contrast matrix for Sex rather than "contr.sum" is 
reported in the "contrasts" attribute):

--------- snip ---------

 > model.matrix(m)
     (Intercept) Sex1
1             1    1
2             1    1
3             1    1

. . .

106           1   -1
107           1   -1
108           1   -1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$Sex
        [,1]
Male      1
Female   -1

 > sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: aarch64-apple-darwin20
Running under: macOS Sonoma 14.6.1

Matrix products: default
BLAS: 
/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 

LAPACK: 
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; 
  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Toronto
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] nlme_3.1-165

loaded via a namespace (and not attached):
[1] compiler_4.4.1 tools_4.4.1    grid_4.4.1     lattice_0.22-6

--------- snip ---------

Although this apparently solves the problem for car::Anova(), the 
problem is likely more general. For example, insight::get_modelmatrix() 
also reports the wrong model matrix for the "lme" model m above.

My suggestion: Include a correct model.matrix.lme() method in the nlme 
package. That could be my version, but I expect that R-core could come 
up with something better.

Thank you,
  John
-- 
John Fox, Professor Emeritus
McMaster University
Hamilton, Ontario, Canada
web: https://www.john-fox.ca/
--