dont_fit_random_effects_with_optimization.R

# This illustrates why Neal's Funnel really messes with us when we're trying 
# to estimate random effects models using maximum likelihood. Ie. why we need
# approximations with lme4, INLA, or need to go full Bayes. The mode is an 
# awful one. 

# Let's simulate some data from a really simple model with 10 
# random intercepts whose variance you want to estimate

N <- 1000
individual <- rep(1:10, each = 100)
theta_sigma <- rnorm(1, 0, 2)
theta <- rnorm(10, 0, exp(theta_sigma/2))
sigma <- rnorm(1, 0, 1)
y <- rnorm(N, theta[individual], exp(sigma/2))


# Now specify the negative log likelihood to pass to the optimizer

negative_log_likelihood_of_this_model <- function(pars) {
  sigma_est <- pars[1]
  sigma_theta_est <- pars[2]
  theta_est <- pars[3:12]
  target <- 0
  # add log of prior for sigmas
  target <- target + log(dnorm(sigma_est, 0, 1)) + log(dnorm(sigma_theta_est, 0, 2))
  
  # Add log of random effects' contribution to likelihood
  target <- target + sum(log(dnorm(theta_est, 0, exp(sigma_theta_est/2))))
  
  # Add log likelihood
  target <- target + sum(log(dnorm(y, theta[individual], exp(sigma_est/2))))
  
  # Return negative log likelihood
  -target
}

# What is the negative penalized likelihood at the correct parameters? 
negative_log_likelihood_of_this_model(c(sigma, theta_sigma, theta))

# Let's try to get a penalized likelihood estimate
penalized_likelihood_fit <- optim(par = rnorm(12), fn = negative_log_likelihood_of_this_model, method = "SANN", control = list(maxit = 3e5, trace = 1))


# Plot the actuals and estimates
library(tidyverse)
data_frame(fit = penalized_likelihood_fit$par, actual = c(sigma,  theta_sigma, theta), parameter = c("sigma", "theta_sigma", rep("theta", 10))) %>% 
  ggplot(aes(x = fit, y = actual, colour = parameter)) +
  geom_point() +
  ggthemes::theme_hc()