5
$\begingroup$

I'm interesting in making predictions using a random-effects model on new data that occur in new groups. Which estimator is most appropriate?

I fit a Poisson random effects model on some fake data:

$$y_i \sim Poisson(\exp(\beta_0 + b_i)), \quad i = 1, \ldots, 100$$ $$b_i \sim N(0, \sigma_b^2)$$

library(lme4)
set.seed(666)

y <- rpois(100, exp( rnorm(100) ))
df <- data.frame(1, y, group = as.character(1:100))
fit <- glmer(y ~ 1 + (1 | group), family = poisson(), data = df)

In the fake data simulation, $\beta_0 = 0$ and $\sigma_b^2 = 1$.

> fit
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: poisson  ( log )
Formula: y ~ 1 + (1 | group)
   Data: df
      AIC       BIC    logLik  deviance  df.resid 
 338.3075  343.5179 -167.1538  334.3075        98 
Random effects:
 Groups Name        Std.Dev.
 group  (Intercept) 0.9156  
Number of obs: 100, groups:  group, 100
Fixed Effects:
(Intercept)  
   -0.02414  

I can think of four ways to estimate the next observation:

  1. (e1) Ordinary marginal mean: Just take the average of $y$. Ignores the model.
  2. (e2) The "unconditional" fixed effect estimate: I set the estimated random effects to $0$ and use $\exp(\hat\beta_0)$ from the model fit. (lme4 permits this if you set re.form = NA in the predict() function)
  3. (e3) Instead of setting the random effects to $0$, I take the average of the observed random effects. The estimator is $\exp(\hat\beta_0 + \frac{1}{100}\sum_{i=1}^{100} \hat{b}_i)$.
  4. (e4) The estimated marginal mean: $\frac{1}{100} \sum_{i=1}^{100} \exp(\hat\beta_0 + \hat{b}_i)$.

Okay, those are easy enough to compare in R:

mean(y) #"ordinary" marginal mean
#> 1.46
exp(fit@beta) #The "Set random effects to 0" estimator
#> 0.97
exp(fit@beta + apply(ranef(fit)$group, 2, mean)) #ave. random effect
#> 1.07
mean(predict(fit, df, type = "response")) # Estimated Marginal Mean
#> 1.47

There's quite a bit of difference between these predictions. I can write a small simulation that compares them. Let's say the best predictor is the one with lowest MSE.

l1 <- l2 <- l3 <- l4 <- c()

for (i in 1:1000){
  y <- rpois(100, exp( rnorm(100) ))
  ynew <- rpois(1, exp(rnorm(1)))
  df <- data.frame(1, y, group = as.character(1:100))
  fit <- glmer(y ~ 1 + (1 | group), family = poisson(), data = df)
  e1 <- mean(y)
  e2 <- exp(fit@beta)
  e3 <- exp(fit@beta + apply(ranef(fit)$group, 2, mean))
  e4 <- mean(predict(fit, df, type = "response"))
  l1[i] <- (ynew - e1) ^ 2
  l2[i] <- (ynew - e2) ^ 2
  l3[i] <- (ynew - e3) ^ 2
  l4[i] <- (ynew - e4) ^ 2
}

mean(l1) # 4.68
mean(l2) # 4.89
mean(l3) # 4.79
mean(l4) # 4.66


Is this correct? Is the estimated marginal mean the 'right' estimator of a new observation in a new group? I'm struggling to understand when, if ever, these unique estimators should be used when dealing with a random effects model.

Edit: My sim was not correct! I took @dimitris' estimator and @mavery's simulation oracle, and it seems like the new estimator is actually the best predictor of a new $y$ from a new group:

ynew <- mean(rpois(10000, exp(rnorm(10000))))

for (i in 1:1000){
  y <- rpois(100, exp( rnorm(100) ))
  df <- data.frame(1, y, group = as.character(1:100))
  fit <- glmer(y ~ 1 + (1 | group), family = poisson(), data = df)
  e1 <- mean(y)
  e2 <- exp(fit@beta)
  e3 <- exp(fit@beta + apply(ranef(fit)$group, 2, mean))
  e4 <- mean(predict(fit, df, type = "response"))
  e5 <- exp(fit@beta + 0.5 * sigma.hat(fit)$sigma$group ^ 2)
  l1[i] <- (ynew - e1) ^ 2
  l2[i] <- (ynew - e2) ^ 2
  l3[i] <- (ynew - e3) ^ 2
  l4[i] <- (ynew - e4) ^ 2
  l5[i] <- (ynew - e5) ^ 2
}

mean(l1) #0.06
mean(l2) #0.44
mean(l3) #0.30
mean(l4) #0.08
mean(l5) #0.05

Which makes me wonder, why does lme4 use unconditional "population level" estimators (via reform = NA, as in $\exp(\hat\beta_0$)), instead of calculating $\exp(\hat\beta_0 + \hat\sigma_b^2/2)$ when a new level is encountered???

Aside, the ordinary marginal mean and the estimated marginal mean seems to be about as good as @dimitis' estimator.

$\endgroup$
7
  • $\begingroup$ Somewhat related, perhaps of interest: What are LS means useful for?. $\endgroup$ – Richard Hardy Oct 26 '19 at 19:19
  • $\begingroup$ Thanks, I've seen that post before. Still thinking about which estimator is useful for this situation. I found Gelman and Hill discuss prediction for non-linear models in their book on pg 274. They do something like estimated marginal means, but with weights. $\endgroup$ – JTH Oct 26 '19 at 22:31
  • $\begingroup$ This post (stackoverflow.com/questions/29259750/…) shows that some have used the fixed effects estimator to make predictions on new levels. It seems that the fixed effects estimator (e2) and the estimated marginal mean (e4) may be the most popular. If that is the case, when should one prefer one over the other? $\endgroup$ – JTH Oct 27 '19 at 14:36
  • $\begingroup$ I do not have the time to dig deeper into your problem, but in general I would be curious to see any problem for which estimated marginal means actually make sense. In my understanding, they correspond to parameters in a fictitious population but not any real one, making their relevance suspect. $\endgroup$ – Richard Hardy Oct 27 '19 at 14:57
  • $\begingroup$ Thanks, I suppose this is a fictitious parameter being estimated: at best, it is the expected response of a yet-to-be-realized group. $\endgroup$ – JTH Oct 28 '19 at 2:04
3
$\begingroup$

If you know that this next observation comes from a specific group $i$, then the best prediction would be $$\exp(\hat \beta_0 + \hat b_i).$$

However, if you do not know from which group it comes from, then it would make more sense to use as prediction the average over the groups, which in your specific random intercepts model would be $$\exp \Bigl ( \hat \beta_0 + \frac{\hat \sigma_b^2}{2} \Bigr ),$$ where $\hat \sigma_b^2$ is the estimated variance of the random intercepts.


EDIT: The marginal mean is:

$$ \begin{eqnarray} E(Y) & = & E\{E(Y \mid b)\} = \int \exp(\beta_0 + b) \; f(b) \, db\\ & = & \exp(\beta_0) \int \exp(b) \frac{1}{\sqrt{2\pi \sigma_b^2}} \exp\Big(-\frac{b^2}{2\sigma_b^2}\Big) \, db\\ & = & \exp(\beta_0) \int \frac{1}{\sqrt{2\pi \sigma_b^2}} \exp\Big(-\frac{b^2 - 2 \sigma_b^2 b}{2\sigma_b^2}\Big) \, db\\ & = & \exp\Big(\beta_0 + \frac{\sigma_b^2}{2}\Big) \int \frac{1}{\sqrt{2\pi \sigma_b^2}} \exp\Big(-\frac{b^2 - 2 \sigma_b^2 b + \sigma_b^4}{2\sigma_b^2}\Big) \, db\\ & = & \exp\Big(\beta_0 + \frac{\sigma_b^2}{2}\Big) \times 1. \end{eqnarray} $$

In the pre-last line we have the integral of the probability density function of a normal distribution with mean and variance equal to $\sigma_b^2$.

$\endgroup$
8
  • $\begingroup$ Thanks, can you help me understand why "averaging over the groups" results in your estimator $\exp(\hat\beta_0 + \hat\sigma^2_b / 2)$, rather than the estimated marginal mean, $\frac{1}{100} \sum_{i=1}^{100} \exp(\hat\beta_0 + \hat{b}_i)$. Or point to a textbook chapter? $\endgroup$ – JTH Oct 28 '19 at 20:49
  • $\begingroup$ @JTH have a look at my edit. $\endgroup$ – Dimitris Rizopoulos Oct 28 '19 at 22:40
  • $\begingroup$ Ah, makes more sense now. I'm confused about why lme4 doesn't do this calculation when it encounters a new level. It will just calculate $\exp(\hat\beta_0)$, the fixed effects estimator. $\endgroup$ – JTH Oct 28 '19 at 23:21
  • $\begingroup$ Other than being a bit more accurate than the ordinary marginal mean, $\bar{y}$, is there any reason to prefer this estimator? $\endgroup$ – JTH Oct 29 '19 at 2:01
  • $\begingroup$ @JTH this estimator works only when you have random intercepts. In more complex random-effects structures you cannot solve the corresponding integral over the random effects. You could still though derive the coefficients that will have a population-average interpretation. This is provided by the marginal_coefs() function of the GLMMadaptive package: drizopoulos.github.io/GLMMadaptive/articles/… $\endgroup$ – Dimitris Rizopoulos Oct 29 '19 at 4:55
1
$\begingroup$

To address the "when is this appropriate" part of your question:

If you fit a mixed model, it may still be appropriate to do inference for the grand mean, provided you interpret it correctly. Specifically, prediction intervals around the grand mean that include the relevant variance components are useful tools for inference. But, given that each observation will come from some group (either one you've already observed or a new one that you haven't), it's unclear what value there actually is from looking at the point estimate; there is no group whose mean will be equal to the grand mean. So the point estimate itself isn't that meaningful, but you can use it to do things that are meaningful. See this tutorial for a discussion of when and why to do estimation with or without accounting for the random effect in a mixed model. This is within the context of an R package for estimating confidence intervals (and other things) for fitted values.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.