Why is the Empirical Bayes estimator not dominating like it's supposed to?

External question — Cross Validated Stack Exchange Author: Huy Pham Original post: https://stats.stackexchange.com/questions/655839 License: CC BY-SA 4.0 — https://creativecommons.org/licenses/by-sa/4.0/ Adaptation: HTML converted to plain text; contact email addresses removed. I am modelling the effect of sites on an outcome. It is precisely the effect of the sites that matter to me, and the other covariates are entered into the model to be controlled for. There is a slight twist where the site can interact with one of the covariates, but I can live with or without that term. The model equation is: $$\hat{Y}=X\beta + \gamma(site + \zeta ) + site\times b +\epsilon$$ Then I am left with the question of how to estimate each ith site's effect $b_i$ . I can simply fit the sites as dummy coded categorical variables and have everything as a fixed effects model. But that is a lot of parameters to estimate, and must surely have an effect on the AIC. So, if I assumed $b \sim N(0,\sigma_\mu)$ and fit it as a mixed effects model, I can thereafter extract $\sigma_\mu$ and $\sigma_\epsilon$ estimate each $b_i$ using the empirical Bayes estimator. Again, slight spanner because there is also a random effect, but just assume the software (lme4) takes care of this. The literature is quite adamant that estimating the site's effects using an empirical Bayes (James-Stein) approach should be better in terms of MSE than the ML (fixed effects regression) method. I should be trading bias for variance. So I attempt both methods. I have made a reproducible example below. I simulate a random intercept and a random effect library(tidyverse) library(lme4) set.seed(314159) #this function myLetters is merely to make as many "names" for levels as I need myLetters <- function(length.out) { a <- rep(letters, length.out = length.out) grp <- cumsum(a == "a") id_container<-vapply(seq_along(a), function(x) paste(rep(a[x], grp[x]), collapse = ""), character(1L)) id_container<-sort(rep(id_container,50)) id_container } gen_data<-function(k){ random_elements <- MASS::mvrnorm(k,mu=c(0,0), #in the sigma which defines the random int and random effs, #the random ints are the first column Sigma = matrix(c(c(20,3), c(3,5)), byrow = T, nrow=2)) id <-myLetters(k) fixed_int <- c(rep(2,k*50)) prev_score <- c(rnorm(k*50,20,10)) gender <- rep(c(1,0,1,0,1),k*50/5) random_int <- c() random_eff <- c() for(i in 1:k){ temp_rand_int <- c(rep(random_elements[i,1],50)) random_int <-c(random_int,temp_rand_int) temp_rand_eff <- c(rep(random_elements[i,2],50)) random_eff <-c(random_eff,temp_rand_eff) } y = (fixed_int + random_int) + prev_score*(10+random_eff) + 4*gender + rnorm(50*k,0,2) sim_data<-data.frame(y=y, id=id, prev_score=prev_score, gender = gender, #put these here for later, #but you won't need them when #feeding it back into lmer random_int = random_int, random_eff = random_eff)} simulated_data<-gen_data(1000) m1<-lmer(y~gender + prev_score + (prev_score|id), REML=F, data=simulated_data) #summary(m1) #this line below will take a solid 5 minutes to run due to 2k parameters being estimated m2<-lm(y~ gender + prev_score + id + prev_score*id, data=simulated_data) #summary(m2) #it's ok if you don't have the performance package, #I use base R's AIC and calculate rmse below #performance::compare_performance(m1,m2) Name | Model | AIC (weights) | AICc (weights) | BIC (weights) | RMSE | Sigma | R2 (cond.) | R2 (marg.) | ICC | R2 | R2 (adj.) m1 | lmerMod | 2.3e+05 (<.001) | 2.3e+05 (<.001) | 2.3e+05 (>.999) | 1.960 | 2.000 | 1.000 | 0.785 | 0.999 | | m2 | lm | 2.1e+05 (>.999) | 2.1e+05 (>.999) | 2.3e+05 (<.001) | 1.959 | 2.000 | | | | 1.000 | 1.000 AIC(m1) [1] 225417.7 AIC(m2) [1] 213156.3 sqrt(mean((fitted.values(m1) - simulated_data$y)^2)) [1] 1.959776 sqrt(mean((fitted.values(m2) - simulated_data$y)^2)) [1] 1.959295 When I compare the models; the fixed effects model is better on RMSE and AIC. Why? I don't really mind the AIC. But on RMSE, I thought the dominance of the Empirical Bayes estimator vs ML estimator was theoretically guaranteed. Note I turned REML off to do the comparison, but it doesn't make a difference if it is turned on. The standard errors of the Empirical Bayes estimates are indeed better, however, that's kind of circular really. The formulas say so, so no surprise there. But if the decrease in variance does not lead to an improvement on prediction then what's the point. And the fixed effects model also produces a decent estimate of the variance of the random effects anyway. So there's not even the value of estimating the random effects using the mixed effects model. store_coefs<-coef(m2)%>%as.data.frame()%>%rownames_to_column() names(store_coefs) <- c("variable", "coef") store_coefs_rand_int<-store_coefs %>% filter(!(variable %in% c("(Intercept)", "gender", "prev_score")) & !grepl(":", variable)) store_coefs_rand_eff<-store_coefs %>% filter(!(variable %in% c("(Intercept)", "gender", "prev_score")) & grepl(":", variable)) var(store_coefs_rand_int$coef) var(store_coefs_rand_eff$coef) At first, I thought it was because I did not have enough sites (using the posterior variance as the standard error only works for sites>>predictors), but I've upped my number of sites and it still stays like this. Changing to a random intercepts only model also doesn't make a difference. Have I simply misunderstood the literature regarding the Empirical Bayes estimator?

External answer — Cross Validated Stack Exchange Author: Huy Pham Original post: https://stats.stackexchange.com/a/655883 License: CC BY-SA 4.0 — https://creativecommons.org/licenses/by-sa/4.0/ Adaptation: HTML converted to plain text; contact email addresses removed. After Christian Henning’s comment I went back and rechecked my simulations. The only coefficient that I really care about is the $b_i$ of each site – and to a lesser degree the random effect $\zeta_i$ . So, I make sure to save them out and subtract back the estimated intercepts (and slopes) for both methods. I still sample them from a normal distribution (defined in my code), but then they become fixed as the value for the site. I had to make adjustments to the fixed effects regression because the intercept would take on a value of a particular site and reference level of any categorical variables, and at the value of 0 for any continuous variables. That will mess up my comparisons. So, I deviance coded all the design matrices – for gender and the sites. I will lose the bottom site, but that’s ok. I had to adjust gender so that it was exactly equal number of boys vs girls. And my site sizes are already balanced. I got rid of the continuous covariate because it is just too tedious to try and adjust the intercept with it still in the model equation. This should ensure that my intercept is the grand mean, and that way all my site’s coefficients are deviances from the grand mean – i.e. the same as the mixed effects model, and the same as the way that I set up the underlying simulated data. Once all that’s done, it still wasn’t clear; sometimes the mixed effects model would have better MSE and sometimes the fixed effects. Then I thought about it a bit more and realized that the advantage of the Empirical Bayes estimator is that, while it also finds the observed average difference between a site and the grand mean, but then instead of stopping there (like the linear regression would), it shrinks that estimate based on the ratio of the variance of sites i.e. the random intercept/slope taken from the mixed effects model $\sigma_\mu$ and the sum of it and variance of the residual error $\sigma_\epsilon$ If that ratio is large (close to 1), then the difference between group and grand mean would dominate and there would be little difference between the mixed effects model and the fixed effects model. However, if that ratio was smaller, i.e. by increasing $\sigma_\epsilon$ then the two sets of estimates would diverge. So, after, changing my residual error’s variance term from 2 to 20; I get the code below, and after having run it a few times it’s clear that the mixed effects model is better at MSE against the site’s true intercepts and slopes, than the fixed effects model. If anybody has any input that would be greatly appreciated. I hope that I have actually solved my problem, and I haven’t just accepted results that affirmed my pre-existing conclusions. library(tidyverse) library(lme4) set.seed(314159) #this function myLetters is merely to make as many "names" for levels as I need myLetters <- function(length.out) { a <- rep(letters, length.out = length.out) grp <- cumsum(a == "a") id_container<-vapply(seq_along(a), function(x) paste(rep(a[x], grp[x]), collapse = ""), character(1L)) id_container<-sort(rep(id_container,60)) id_container } gen_data<-function(k){ random_elements <- MASS::mvrnorm(k,mu=c(0,0), #in the sigma which defines the random int and random effs, #the random ints are the first column Sigma = matrix(c(c(20,3), c(3,5)), byrow = T, nrow=2)) id <-myLetters(k) fixed_int <- c(rep(2,k*60)) prev_score <- c(rnorm(k*60,20,10)) #i probably should have effects coded this when calling the lm to keep it neat, oh well. gender <- rep(c(1,1,1,-1,-1,-1),k*60/6) random_int <- c() random_eff <- c() for(i in 1:k){ temp_rand_int <- c(rep(random_elements[i,1],60)) random_int <-c(random_int,temp_rand_int) temp_rand_eff <- c(rep(random_elements[i,2],60)) random_eff <-c(random_eff,temp_rand_eff) } y = (fixed_int + random_int) + #prev_score*(10+random_eff) + (4+random_eff)*gender + rnorm(60*k,0,20) sim_data<-data.frame(y=y, id=id, prev_score=prev_score, gender = gender, #put these here for later, #but you won't need them when #feeding it back into lmer random_int = random_int, random_eff = random_eff)} #only do 120 clusters to save time simulated_data<-gen_data(120) m1<-lmer(y~ gender+(gender|id), REML=T, data=simulated_data) #summary(m1) m2<-lm(y~ gender+id+gender*id , data=simulated_data, contrasts = list(id = contr.sum(120))) #summary(m2) #setup to compare MSE for random INTERCEPT true_value_lookups<-simulated_data %>% group_by(id)%>% summarise(true_rand_int = unique(random_int), true_rand_eff = unique(random_eff))%>% as.data.frame() m1_rand_ints <- as.data.frame(ranef(m1))%>% pivot_wider(names_from = term, values_from = c(condval, condsd))%>% pull(`condval_(Intercept)`) m1_rand_ints <-m1_rand_ints[-120] m2_output <- as.data.frame(summary(m2)$coefficients) m2_rand_ints <- m2_output%>%rownames_to_column()%>%filter(grepl("id", rowname)) mean((m1_rand_ints - true_value_lookups[-120,2])^2) mean((m2_rand_ints$Estimate - true_value_lookups[-120,2])^2) #setup to compare MSE for random SLOPE m1_rand_effs <- as.data.frame(ranef(m1))%>% pivot_wider(names_from = term, values_from = c(condval, condsd))%>% pull(`condval_gender`) m1_rand_effs <-m1_rand_effs[-120] m2_output <- as.data.frame(summary(m2)$coefficients) m2_rand_effs <- m2_output%>%rownames_to_column()%>%filter(grepl("\\:", rowname)) mean((m1_rand_effs - true_value_lookups[-120,3])^2) mean((m2_rand_effs$Estimate - true_value_lookups[-120,3])^2)