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

[D. Rizopoulos]

You could fit the same model also using the GLMMadaptive package in 
which the predict() method can produce predictions for the mean subject 
(i.e., the one with random effects values equal to 0), marginal 
predictions (i.e., averaged over the subjects), and subject-specific 
predictions.

For more info you can have a look at:

https://drizopoulos.github.io/GLMMadaptive/

https://drizopoulos.github.io/GLMMadaptive/articles/Methods_MixMod.html#predictions

and potentially also

https://drizopoulos.github.io/GLMMadaptive/articles/Dynamic_Predictions.html

if you're interested in dynamic predictions.

Best,
Dimitris



On 10/23/2018 7:46 PM, John Wilson wrote:
> 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]]
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 

-- 
Dimitris Rizopoulos
Professor of Biostatistics
Department of Biostatistics
Erasmus University Medical Center

Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands
Tel: +31/(0)10/7043478
Fax: +31/(0)10/7043014
Web (personal): http://www.drizopoulos.com/
Web (work): http://www.erasmusmc.nl/biostatistiek/
Blog: http://iprogn.blogspot.nl/