[R-sig-ME] population-level predict glmmtmb with poly()

[John Wilson]

Hello,

I'm working on a glmmtmb() model with multiple continuous and categorical
predictors. Two of the predictors are orthogonal polynomials (I just saw
that the package was updated yesterday (!) to correctly handle those). One
of the polynomials has an interaction with another covariate.

Since predict(re.form = 0) doesn't work just yet and one has to use
the model.matrix(lme4::nobars(formula(mod1)[-2]), newdata) approach - how
do I get the correct polynomial predictions out? It looks like my results
depend on how I structure the newdata data frame - when I use
expand.grid(), the predictions are wrong, but when I subset the original
data, the predictions are correct.

Thanks so much!
John

Here's an example:

library(ggplot2)
library(glmmTMB)

set.seed(0)
x <- 1:20
z <- sample(c("a", "b"), length(x), replace = TRUE)
y <- round(5 * 2*x + 3 * x^2 + 0.1 * x^3 + rnbinom(length(x), 10, 0.03))
group <- sample(c("i", "ii"), length(x), replace = TRUE)
df <- data.frame(x = x, y = y, z = z, group = group)
m <- glmmTMB(y ~ poly(x, 3) * z +
(1 | group),
family = nbinom2,
data = df)
# prediction on a new grid
newdata <- expand.grid(x = 1:20, z = unique(df$z))
X.cond = model.matrix(lme4::nobars(formula(m)[-2]), newdata)
beta.cond = fixef(m)$cond
newdata$Pred1 = as.numeric(X.cond %*% beta.cond)
# prediction in original data frame
X.cond = model.matrix(lme4::nobars(formula(m)[-2]))
beta.cond = fixef(m)$cond
df$Pred1 = as.numeric(X.cond %*% beta.cond)

# the newdata preds are obviously off
ggplot(df) +
geom_point(aes(x = x, y = y, colour = z)) +
geom_line(data = newdata, aes(x = x, y = exp(Pred1), colour = z), size = 2)
+
geom_line(aes(x = x, y = exp(Pred1), colour = z, linetype = group))

# if the new grid is defined like this, then the predictions are ok
newdata <- unique(select(df, x, z))
X.cond = model.matrix(lme4::nobars(formula(m)[-2]), newdata)
beta.cond = fixef(m)$cond
newdata$Pred1 = as.numeric(X.cond %*% beta.cond)
ggplot(df) +
geom_point(aes(x = x, y = y, colour = z)) +
geom_line(data = newdata, aes(x = x, y = exp(Pred1), colour = z), size = 2)
+
geom_line(aes(x = x, y = exp(Pred1), colour = z, linetype = group))

	[[alternative HTML version deleted]]