# 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?. Oct 26, 2019 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, 2019 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, 2019 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. Oct 27, 2019 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, 2019 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, 2019 at 20:49 • @JTH have a look at my edit. Oct 28, 2019 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, 2019 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, 2019 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/… Oct 29, 2019 at 4:55