Dear list members,
After further testing, I found that the following simplified version of
model.matrix.lme(), which omits passing xlev to the default method, is
more robust. The previous version generated spurious warnings in some
circumstances.
model.matrix.lme <- function(object, ...){
data <- object$data
if (is.null(data)){
NextMethod(formula(object), data=eval(object$call$data),
contrasts.arg=object$contrasts)
} else {
NextMethod(formula(object), data=data,
contrasts.arg=object$contrasts)
}
}
Best,
John
On 2024-09-20 12:49 p.m., John Fox wrote:> Caution: External email.
>
>
> 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/
> --
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel