Skip to content

predict() with varying betas (newparams)

2 messages · Hans Ekbrand, Ben Bolker

#
Dear list,

I would like to predict() in lme4 using new betas. In order to try out
the option newparams I tried the following:
[1] 0.1445534
theta parameter vector not named: assuming same order as internal vector
beta parameter vector not named: assuming same order as internal vector
[1] 0.05619105

I don't understand what theta is, but I did find a slot in the fitted
object, my.fit, so I just passed that vector and hoped for the best :-)

I would have expected the results to be identical. What am I missing?

The long term objective is to calculate average marginal effects from
a glmer() object by using predict() with betas sampled from a
multivariate distribution, like

  library(MASS) ## 
  n.samples <- 100
  betas <- matrix(mvrnorm(n = n.samples, mu = my.fit at beta, Sigma = vcov(my.fit)), ncol = n.samples, byrow=TRUE)
  new.df <- my.fit at frame
  new.df$Some.var <- my.fit at frame$Some.var - 0.5
  mean(predict(my.fit, newdata = my.fit at frame, type = "response", newparams = list(theta = my.fit at theta, beta = betas[, 1])))
  new.df$Some.var <- my.fit at frame$Some.var + 0.5
  mean(predict(my.fit, newdata = my.fit at frame, type = "response", newparams = list(theta = my.fit at theta, beta = betas[, 1])))
  ... and again with betas[, 2] and so on.
#
Hans Ekbrand <hans.ekbrand at ...> writes:
[snip]

Good question.  I'll look into it, I don't understand what's going
on here since for the most part we've been using newparams=
and it has looked sensible -- and your application looks sensible.

A reproducible example:

library(lme4)
gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
             data = cbpp, family = binomial)
mean(predict(gm1,newdata=model.frame(gm1),type="response"))
mean(predict(gm1,newdata=model.frame(gm1),
             newparams=getME(gm1,c("theta","beta")),type="response"))

mean(p0 <- predict(gm1,newdata=model.frame(gm1)))
mean(p1 <- predict(gm1,newdata=model.frame(gm1),
                   newparams=list(beta=fixef(gm1))))