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

[John Wilson]

So far it's looking just grand! +1 for "it  wasn't documented until about
120 seconds ago".
Just to push my luck here - in my full model, I have autocorrelation
structures. When I only set group = NA, the predict() function gave me an
error due to lack of the autocorrelation variable in the new data. I then
crossed my fingers and did the same thing as for group - a column with NAs.
The predict() function had no problem with that. However, I wanted to see
if that's a (semi) legit way to doing it.

On Tue, Oct 23, 2018 at 8:22 PM Ben Bolker <bbolker using gmail.com> wrote:

>
>   Naive prediction with data-dependent bases **will not work** with new
> input variables. (This is a general R/model.matrix() thing ... not
> glmmTMB's fault.)  The current development version of glmmTMB does some
> magic (see ?makepredictcall,
> https://developer.r-project.org/model-fitting-functions.html for more
> detail) that makes this work.
>
>   It wasn't documented until about 120 seconds ago, but in order to do
> population-level predictions with predict() all you need to do is set
> the group variable to NA.  This has less flexibility than the re.form=
> argument in lme4 (which allows you to drop a *subset* of the random
> effects terms for a given grouping variables), but it does handle the
> most common use case ...
>
>   Does this work for you (with devel version of glmmTMB) ?
>
> # prediction on a new grid
> newdata <- expand.grid(x = 1:20, z = unique(df$z), group=NA)
> newdata$y <- predict(m,newdata=newdata, type="response")
>
> (ggplot(df, aes(x,y, colour=z))
>     + geom_point()
>     + geom_line(data = newdata, size =2)
> )
>
>
> On 2018-10-23 6:48 p.m., John Wilson wrote:
> > Thank you, Dimitris. I poked around a bit but didn't see mention of
> > autocorrelation structures - are those supported as well?
> >
> > I was wondering if this issue was related to the problem logged here:
> > https://github.com/glmmTMB/glmmTMB/issues/402. Can anyone comment on
> > whether what I'm seeing is a bug or just a user (=me) mistake? When I
> > change the toy example to be a simple linear form (rather than poly()),
> the
> > two predictions are identical, so I'm pretty sure it's a poly() issue...
> >
> > On Tue, Oct 23, 2018 at 4:31 PM D. Rizopoulos <d.rizopoulos using erasmusmc.nl
> >
> > wrote:
> >
> >> 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/
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-mixed-models using r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]