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