Is a random intercept model exactly the same as a linear model with dummy variable? I have a large household dataset for 28 different countries. Upto now I have used OLS (R command lm()) with 'country' as a dummy variable, to control for unobserved heterogeneity between countries.
Is this exactly the same as using a linear mixed effects model with country random effects for the intercept (lmer() with + (1|country)?
I appear to get exactly the same results when I make predictions with both models (using predict()). I was wondering whether this is normal or whether I have overlooked something?
 A: The difference between the two is that dummy variables do not share information with each other whereas random effects do via a distribution. The dummy variable model estimates 28 different effects. The mixed model estimates a distribution from which the 28 effects are drawn from. They share information because they all contribute to the estimation of the variation. This is my favorite blog post on the topic.
A: The blog post that @statwonk mentioned is superb and worth reading in full. However, for an alternative simple demonstration, consider just the intercepts that we get from these two models, when using the sleepstudy dataset:

*

*Model 0: Dummy variables for each subject (group)

*Model 1: Using random effect for each subject (group)

Here's the distribution of intercepts in each case:

As @statwonk mentioned, this difference is because in the random effects model, there is a "soft constraint" that ensures that the intercepts are all drawn from a common distribution. There is no similar constraint in the dummy variables model.
Here's the code for generating the plot:
library(lme4)
library(magrittr)
library(dplyr)


# ### Model 0: Using dummy variables for group 
m0_dummy <- lm(Reaction ~ Subject, data = sleepstudy)

# Extract coeffiencts: 
df_m0_coef <- m0_dummy$coefficients %>% as.data.frame()

# Convert intercepts from relative to absolute values:  

# By default, one subject is assigned as the baseline, and all other
# subject intercepts are shown as deviations from the baseline. 
# Here, we convert all subject intercepts to absolute values. 

baseline_intercept <- df_m0_coef[1,1]
df_m0_coef <- df_m0_coef %>% 
  mutate(group_intercept = ifelse(. != baseline_intercept,
                                  baseline_intercept + .,
                                  .)) %>% 
  select(group_intercept)


# ### Model 1: Using random effect for group
m1_random_effects <- lmer(Reaction ~ 1 + (1 | Subject), data = sleepstudy)

# Extract coeffiencts:
df_m1_coef <- coef(m1_random_effects)$Subject



# ### Comparison plots 
x_range <- c(150, 400)
y_range <- c(0, 10) 

par(mfrow = c(1, 2))

df_m0_coef$group_intercept %>% 
  hist(xlim = x_range, ylim = y_range, breaks = 10, 
       main = 'Dummy model: \nEstimates of group intercepts')

df_m1_coef$`(Intercept)` %>% 
  hist(xlim = x_range, ylim = y_range, breaks = 10, 
       main = 'Random effects model: \nEstimates of group intercepts')


