0

I am running a linear mixed effects models using the "nlme" package looking at stress and lifestyle as predictors of change in cognition over 4 years in a longitudinal dataset. All variables in the model are continuous variables. I am able to create the model and get the summary statistics using this code:

mod1 <- lme(MS ~ age + sex + edu + GDST1*Time + HLI*Time + GDST1*HLI*Time, random= ~ 1|ID, data=NuAge_long, na.action=na.omit)
summary(mod1)

I am trying to use the "interactions" package to probe the 3-way interaction:

sim_slopes(model = mod1, pred = Time, modx = GDST1, mod2 = HLI, data = NuAge_long)

but am receiving this error:

Error in if (tcol == "df") tcol <- "t val." : argument is of length zero

I am also trying to plot the interaction using the same "interactions" package:

interact_plot(model = mod1, pred = Time, modx = GDST1, mod2 = HLI, data = NuAge_long)

and am receiving this error:

Error in UseMethod("family") : no applicable method for 'family' applied to an object of class "lme"

I can't seem to find what these errors mean and why I'm getting them. Any help would be appreciated!

danielledamico
  • 1
  • 1

1 Answers1

0

From ?interactions::sim_slopes:

The function is tested with ‘lm’, ‘glm’, ‘svyglm’, ‘merMod’, ‘rq’, ‘brmsfit’, ‘stanreg’ models. Models from other classes may work as well but are not officially supported. The model should include the interaction of interest.

Note this does not include lme models. On the other hand, merMod models are those generated by lme4::[g]lmer(), and as far as I can tell you should be able to fit this model equally well with lmer():

library(lme4)
mod1 <- lmer(MS ~ age + sex + edu + GDST1*Time + HLI*Time + GDST1*HLI*Time 
               + (1|ID), data=NuAge_long)

(things will get harder if you want to specify correlation structures, e.g. correlation = corAR1(), which works for lme() but not lmer() ...)

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • Thank you - this worked and I'm able to get the simple slopes and a nice graph from the 3-way interaction. In addition to the 9 simple slopes of Time on 3MS at each level of the moderators (-1SD, mean, +1SD HLI and -1SD, mean, +1SD GDST1), I want to get the effect of GDST1*Time at each level of HLI. Essentially I wanted to know if GDST1 is associated with change in 3MS (the DV) over time at -1SD, mean, and +1SD HLI --- is there a way to do this with the "interaction" package? Thanks so much again, I hope this makes sense. – danielledamico Aug 10 '22 at 00:21
  • I don't know. I would imagine that either this package, or the `emmeans` package, would do it ... – Ben Bolker Aug 10 '22 at 14:39