munklak simulation

munklak_sim.R

library(lfe)
library(lme4)
library(MASS)
library(tidyverse)
library(modelsummary)

set.seed(290653)
sim <- function(Nobs = 100, Ngrp = 50, cor_x_u = .6, x_sd = .2, y_sd = 1, ...) {
    group <- mvrnorm(Ngrp, c(0, 0), matrix(c(1, cor_x_u, cor_x_u, 1), ncol = 2)) %>%
             data.frame %>% setNames(c('U', 'Ucor')) %>%
             mutate(group = 1:Ngrp,
					Z = rnorm(n()))
    dat <- expand_grid(group = 1:Ngrp, observation = 1:Nobs) %>%
           left_join(group, by = 'group') %>%
           mutate(X = rnorm(Ngrp * Nobs, mean = Ucor, sd = x_sd),
                  e = rnorm(Ngrp * Nobs, mean = 0, sd = 1),
                  Y = X + U + Z + e) %>%
 		   group_by(group) %>%
		   mutate(Xbar = mean(X)) 
    return(dat)
}

dat <- sim()
mod <- list()
mod[['FE']] <- felm(Y ~ X | group, data = dat)
mod[['MLM']] <- lmer(Y ~ X + Z + (1 | group), data = dat)
mod[['Mundlak']] <- lmer(Y ~ I(X - Xbar) + Xbar + Z + (1 | group), data = dat)
msummary(mod, coef_omit = 'group', fmt = '%.5f')