[R-sig-ME] Error with nlme and fixed effects defined with a minus 1

[Phillip Alday]

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
>