Confidence Intervals for nlme

[strengejacke]

I also tried ggeffects a while ago and noticed that it does not handle nlme objects (at least as I expected). The function 'nlraa::predict_nlme' produced confidence bands and could be called by 'ggeffects'

Originally posted by @femiguez in easystats/insight#309 (comment)

[femiguez]

Let me know if I can assist with this. I have some thoughts about confidence bands in NLMEs (https://femiguez.github.io/nlraa-docs/Confidence-Bands.html).

[ptompalski]

Not able to help with this, but wanted to say it would be fantastic to have this available!

[femiguez]

@ptompalski Are you willing to give 'nlraa::predict_nlme' a try? You can also construct confidence bands using bootstrapping, which is a bit more computationally intensive.

[strengejacke]

Just revisiting some older issues... If I understand right, getting confidence or prediction intervals can be easily achieved by calling nlraa::predict_nlme()? Does it also work for a data grid only, i.e. intervals for a certain set of values only (and not all observations)?

[femiguez]

Yes, I think nlraa::predict_nlme will work as you expect. It was modeled as close as possible to 'predict.lm' with same argument names. If your second question is about whether it can make predictions for new data, then the answer is 'yes' via the 'newdata' argument

[strengejacke]

confidence intervals seems to be in line with nrlaa, but prediction intervals don't work.

library(ggeffects)
library(nlme)
data(sleepstudy, package = "lme4")
fm2 <- lme(Reaction ~ Days, data = sleepstudy, random = ~ 1 | Subject)
nlraa::predict_lme(fm2, newdata = data_grid(fm2, "Days"), interval = "prediction")
#> Error in simulate_lme_one(object, psim = psim, data = data, ...): psim = 2 should only be used for the deepest level of the hierarchy

^{Created on 2024-05-07 with reprex v2.1.0}

[femiguez]

@strengejacke Thanks for coming back to this! The example above also fails because data_grid does not seem to work for 'lme' objects

library(modelr)
library(nlme)
data(sleepstudy, package = "lme4")
fm2 <- lme(Reaction ~ Days, data = sleepstudy, random = ~ 1 | Subject)
data_grid(fm2, "Days")
Error in UseMethod("expand") : 
  no applicable method for 'expand' applied to an object of class "lme"

The error you are seeing is because by default, predict_lme uses fixed effects only, but 'predictions' in a mixed model make sense only if you include random effects to a certain extent. In 'lm' a prediction is a statement about 'observations' and not 'mean' effects. Of course, the word prediction can be interpreted in many ways, here I'm referring to a 'prediction interval'.

I had a problem in the function for 'newdata' and 'prediciton' which is different from the issue above. I fixed this in github in version of the package 1.9.9. This now works,

library(nlme)
library(nlraa)
library(ggplot2)
data(sleepstudy, package = "lme4")

fm2 <- lme(Reaction ~ Days, data = sleepstudy, random = ~ 1 | Subject)
ndat <- expand.grid(Days = 0:9, Subject = unique(sleepstudy$Subject))
prds <- predict_lme(fm2, newdata = ndat, plevel = 1, interval = "prediction")

ndatA <- cbind(ndat, prds)

ggplot() + 
  facet_wrap(~ Subject) + 
  geom_point(data = sleepstudy, aes(x = Days, y = Reaction)) + 
  geom_line(data = ndatA, aes(x = Days, y = Estimate)) + 
  geom_ribbon(data = ndatA, aes(x = Days, ymin = Q2.5, ymax = Q97.5, fill = "purple"),
              alpha = 0.3)

The default for 'predict_lme' is plevel = 0 (see documentation for predict.lme with confusing argument 'level' which is not a probability level but the 'level' of the hierarchical mixed model at which to return predictions). 'level' in predict.lme is equivalent to 'plevel' in 'predict_lme'. The value of 0 refers to fixed effects only, but again this does not make sense for 'prediction intervals' because we are asking for something at 'deeper' level than the fixed effects. In this particular example, plevel = 1 builds prediction intervals at the level of observation within a subject. I'm open to making improvements/changes to nlraa if this makes it easier to work with ggeffects.