Example mixed-effects model integration and optimization in R

lme.R

set.seed(138593)
N <- 100
x <- runif(N, -2, 2)
a <- 1
b <- 1.5
sigma <- 0.3
tau <- 0.4
re_levels <- gl(10, 10)
re <- rnorm(10, 0, tau)
y_true <- a + re[re_levels] + b * x
y <- rnorm(N, mean = y_true, sd = sigma )  
dat <- data.frame(x = x, y = y)

plot(x, y)
for (i in seq_along(re)) {
  abline(a = a + re[i], b = b)
}

ll_fixed_ran <- function(par, re_val, x, y) {
  eta <- re_val + par[1] + par[2] * x
  sum(dnorm(y, mean = eta, sd = exp(par[3]), log = TRUE)) +
    dnorm(re_val, 0, exp(par[4]), log = TRUE)
}
nll_fixed_ran <- function(par) {
  re <- seq(-10 * exp(par[4]), 10 * exp(par[4]), length.out = 100)
  width <- re[2] - re[1]
  nll <- sapply(seq_along(levels(re_levels)), function(j) {
    ll <- sapply(re, function(re_val) {
      ll_fixed_ran(par, re_val, 
        x = dat$x[which(re_levels == j)], y = dat$y[which(re_levels == j)])
    })
    -log(sum(exp(ll) * width)) # integrate
    })
  sum_nll <- sum(nll)
  cat(sum_nll, "\n")
  sum_nll
}
lm_mle <- nlminb(c(0, 0, 0, 0), nll_fixed_ran)

lm_mle
round(lm_mle$par[1:2], 2)
round(exp(lm_mle$par[3]), 2)
round(exp(lm_mle$par[4]), 2)

m <- lme4::lmer(y ~ x + (1|re_levels), REML = FALSE)
arm::display(m)