This question just migrated from crossvalidated, as it is indeed more a programming question.
I have tried a lot of things (all kind of way to provide info to the newdata-argument), from which I present some below. The summary: everything works with lme4, but only the first works with nlme.
The question:
I need to use nlme rather than lme4 because I want to be able to account for heteroscedasticity. But for some reason, comparisons in the marginaleffects package keeps producing the same error, no matter what I try: Error: Unable to compute predicted values with this model. You can try to supply a different dataset to the newdata argument.
Below some some toy-data that replicates the problem. Only the first comparison-command works, the others all give the error-message.
base <- expand.grid(time = 0:5,
Treat = LETTERS[1:3], stringsAsFactors = FALSE)
base <- base%>%bind_rows(base, base)%>%group_by(time, Treat)%>%
mutate(repl = row_number(),
replicate = sprintf("%s-%s", Treat, repl))%>%ungroup()
re <- data.frame(replicate = unique(base$replicate),
re = rnorm(9))
datasim <- base%>%left_join(re)%>%
mutate(out = re + rnorm(54, sd = 1 + time))%>%
mutate(ctime = as.character(time))
fit <- nlme::lme(out ~ ctime*Treat,
random = ~ 1| replicate,
weights = varIdent(form = ~ 1 | ctime),
data = datasim)
comparisons(fit)
comparisons(fit,
variables = "Treat")
comparisons(fit,
variables = list(Treat = "pairwise"))
comparisons(fit,
variables = list(Treat = "reference"),
)
comparisons(fit,
variables = list(Treat = "pairwise"),
newdata = datagrid(ctime = "5"))
comparisons(fit,
variables = list(Treat = "pairwise"),
newdata = datagrid(ctime = "5", model = fit))
comparisons(fit,
variables = list(Treat = "pairwise"),
newdata = datagrid(ctime = "5", replicate = "A-1"))
If I replace the lme-call by an lmer-call, all of the comparisons run as they should.
Convert the Treat variable to a factor before fitting your model.
This is (arguably) an upstream problem with the handling of character variables in predict.lme(), which breaks when newdata does not include all levels of the predictor.
Load packages and prep data:
library(tidyverse)
library(nlme)
library(marginaleffects)
base <- expand.grid(
time = 0:5,
Treat = LETTERS[1:3],
stringsAsFactors = FALSE)
base <- base %>%
bind_rows(base, base) %>%
group_by(time, Treat) %>%
mutate(
repl = row_number(),
replicate = sprintf("%s-%s", Treat, repl)) %>%
ungroup()
re <- data.frame(
replicate = unique(base$replicate),
re = rnorm(9))
datasim <- base %>%
left_join(re, by = "replicate") %>%
mutate(out = re + rnorm(54, sd = 1 + time)) %>%
mutate(ctime = as.character(time)) %>%
mutate(Treat = factor(Treat))
Fit model and comparisons:
fit <- nlme::lme(out ~ ctime * Treat,
random = ~ 1 | replicate,
weights = varIdent(form = ~ 1 | ctime),
data = datasim)
comparisons(fit, variables = "Treat")
#
# Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
# Treat B - A -0.0365 1.03 -0.0355 0.972 0.0 -2.05 1.978
# Treat B - A -0.9755 1.52 -0.6408 0.522 0.9 -3.96 2.008
# Treat B - A 0.9695 2.22 0.4363 0.663 0.6 -3.39 5.325
# Treat B - A 0.1212 4.38 0.0277 0.978 0.0 -8.46 8.707
# Treat B - A -6.6345 3.08 -2.1565 0.031 5.0 -12.66 -0.605
# --- 98 rows omitted. See ?avg_comparisons and ?print.marginaleffects ---
# Treat C - A -1.7364 1.52 -1.1406 0.254 2.0 -4.72 1.247
# Treat C - A -1.7145 2.22 -0.7715 0.440 1.2 -6.07 2.641
# Treat C - A 0.2368 4.38 0.0541 0.957 0.1 -8.35 8.823
# Treat C - A -1.4930 3.08 -0.4853 0.627 0.7 -7.52 4.537
# Treat C - A 0.5092 6.14 0.0829 0.934 0.1 -11.53 12.551
# Columns: rowid, term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted, out, ctime, Treat, replicate
# Type: response