[R-sig-ME] Error with nlme and fixed effects defined with a minus 1
Phillip Alday
me @end|ng |rom ph||||p@|d@y@com
Fri Aug 13 04:08:24 CEST 2021
This has more to do with how R expands formulae and contrast coding than
anything specific to nlme, I think. I suspect that if you did this with
just a normal (non mixed) linear model, the same problem would arise.
When you have -1, the intercept is suppressed, which in the presence of
a categorical variable means that all k levels of that variable are
represented as coefficients instead of k-1 (with the final contrast
thrown in with the intercept). This is of course not ideal behavior and
could probably be addressed or worked around in nlme, but it seems like
a bit of an edge case. Maybe setting the contrasts manually will help?
contrasts(simdata2$treatment) <- contrasts(simdata$treatment)
(Good call on setting the factor levels in simdata2!)
On 12/8/21 6:59 pm, bill using denney.ws wrote:
> Hello,
>
>
>
> I think that I have found a bug in nlme. When I generate a model with fixed
> effects defined with a factor - 1 (see model_minus1 below) and I simulate
> with level=0 in a data.frame that does not have all the levels (see predict
> into simdata2 with model_minus1 below), there is an error. But, if the
> model is defined without the minus 1 or all factor levels are represented
> within the simulation data, there is no error.
>
>
>
> Did I find a bug in nlme?
>
>
>
> Thanks,
>
>
>
> Bill
>
>
>
> ``` r
>
> library(nlme)
>
> simdata <-
>
> merge(
>
> merge(
>
> data.frame(treatment=factor(c("A", "B"))),
>
> data.frame(ID=factor(1:10))
>
> ),
>
> data.frame(time=1:10)
>
> )
>
> set.seed(5)
>
> simdata$obs <- rnorm(nrow(simdata))
>
> model_minus1 <- nlme(obs~e0, data=simdata, fixed=e0~treatment - 1,
> random=e0~1|ID, start=c(e0=c(0, 0)))
>
> model <- nlme(obs~e0, data=simdata, fixed=e0~treatment, random=e0~1|ID,
> start=c(e0=c(0, 0)))
>
>
>
> # Generate new data with the correct treatment factor levels, but only one
> of
>
> # the levels represented within the data.
>
> simdata2 <-
>
> merge(
>
> data.frame(treatment=factor("A", levels=c("A", "B"))),
>
> data.frame(time=1:10)
>
> )
>
> # Generate new data with all factor levels represented
>
> simdata3 <-
>
> merge(
>
> data.frame(treatment=factor(c("A", "B"))),
>
> data.frame(time=1:10)
>
> )
>
> predict(model, newdata=simdata2, level=0)
>
> #> [1] 0.02690558 0.02690558 0.02690558 0.02690558 0.02690558 0.02690558
>
> #> [7] 0.02690558 0.02690558 0.02690558 0.02690558
>
> #> attr(,"label")
>
> #> [1] "Predicted values"
>
> # Without all levels in newdata$treatment, it fails
>
> predict(model_minus1, newdata=simdata2, level=0)
>
> #> Error in pars[, nm] <- f %*% beta[fmap[[nm]]]: number of items to replace
> is not a multiple of replacement length
>
> # With all levels in newdata$treatment, it works
>
> predict(model_minus1, newdata=simdata3, level=0)
>
> #> [1] 0.02690558 0.02123989 0.02690558 0.02123989 0.02690558 0.02123989
>
> #> [7] 0.02690558 0.02123989 0.02690558 0.02123989 0.02690558 0.02123989
>
> #> [13] 0.02690558 0.02123989 0.02690558 0.02123989 0.02690558 0.02123989
>
> #> [19] 0.02690558 0.02123989
>
> #> attr(,"label")
>
> #> [1] "Predicted values"
>
> ```
>
>
>
> <sup>Created on 2021-08-12 by the [reprex
> package](https://reprex.tidyverse.org) (v2.0.0)</sup>
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list