# Predicting a new observation: marginal mean, estimated marginal mean, or fixed effects estimator?

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. • Somewhat related, perhaps of interest: What are LS means useful for?. – Richard Hardy Oct 26 '19 at 19:19 • 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. – JTH Oct 26 '19 at 22:31 • 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? – JTH Oct 27 '19 at 14:36 • 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. – Richard Hardy Oct 27 '19 at 14:57 • Thanks, I suppose this is a fictitious parameter being estimated: at best, it is the expected response of a yet-to-be-realized group. – JTH Oct 28 '19 at 2:04 ## 2 Answers 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$$. • 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? – JTH Oct 28 '19 at 20:49 • @JTH have a look at my edit. – Dimitris Rizopoulos Oct 28 '19 at 22:40 • 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. – JTH Oct 28 '19 at 23:21 • Other than being a bit more accurate than the ordinary marginal mean,$\bar{y}\$, is there any reason to prefer this estimator? – JTH Oct 29 '19 at 2:01
• @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/… – Dimitris Rizopoulos Oct 29 '19 at 4:55