You can use still compute the correct design matrix when using poly() on
newdata by working with a terms objects that has an appropriately
defined the 'predvars' attribute. For an example, check the following:
# data and new data for prediction
DF <- data.frame(y = rnorm(10), x = rnorm(10))
newDF <- data.frame(x = rnorm(10))
# *wrong* design matrix X on new data
X_new1 <- model.matrix(~ poly(x, 3), data = newDF)
# correct design matrix on new data
termsX <- terms(model.frame(y ~ poly(x, 3), data = DF))
# check predvars attribute
attr(termsX, "predvars")
X_new2 <- model.matrix(delete.response(termsX), data = newDF)
head(X_new1)
head(X_new2)
Best,
Dimitris
On 10/30/2018 12:58 AM, John Wilson wrote:
Hello,
I've been using the newly documented predict() with group = NA to predict
population-level values, as per the thread here (
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2018q4/027305.html). A
follow-up question: in the case of a zero-inflated model, how would I go
about to get the 95% CIs for response predictions? Obviously, I can get
them for the count and the zero-component, and without poly(), I would
follow the example in Brooks et al (2017).
However, I have a poly() predictor; how do I get the 95% CIs if I can't
model.matrix naively when I have a poly() in the model? The models take a
while to converge, so I don't want to run a full bootstrapping either, if
at all possible.
Thank you!
John
Here's a toy dataset:
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))
y[c(2, 3, 5, 11, 13, 19)] <- 0
group <- sample(c("i", "ii"), length(x), replace = TRUE)
df <- data.frame(x = x, y = y, z = z, group = group)
ggplot(df) +
geom_point(aes(x = x, y = y, colour = z))
m <- glmmTMB(y ~ poly(x, 3) * z +
(1 | group),
zi = ~ z,
family = nbinom1,
data = df)
# prediction on a new grid
newdata <- expand.grid(x = 1:20, z = unique(df$z), group = NA)
newdata$Pred <- predict(m, type = "response")
### now how to add CIs?
ggplot(df) +
geom_point(aes(x = x, y = y, colour = z)) +
geom_line(data = newdata, aes(x = x, y = Pred, colour = z, group = z),
= 2) +
facet_wrap(~ z)
[[alternative HTML version deleted]]