suppressPackageStartupMessages(library(lme4))
set.seed(20261006)
J <- 80L
m <- 50L
rho <- 0.02
cluster <- factor(rep(seq_len(J), each = m))
u <- rnorm(J, sd = sqrt(rho))
y <- 10 + u[cluster] + rnorm(J * m, sd = sqrt(1 - rho))
dat <- data.frame(y, cluster)
naive <- lm(y ~ 1, data = dat)
clustered <- lmer(y ~ 1 + (1 | cluster), data = dat, REML = TRUE)
tau2 <- as.numeric(VarCorr(clustered)$cluster[1, 1])
sigma2 <- sigma(clustered)^2
icc <- tau2 / (tau2 + sigma2)
print(c(J = J, m = m, n = nrow(dat), generating_ICC = rho,
        fitted_ICC = icc, fitted_between_variance = tau2,
        fitted_within_variance = sigma2,
        naive_mean = coef(naive)[1], mixed_mean = fixef(clustered)[1],
        naive_SE = sqrt(vcov(naive)[1, 1]),
        mixed_SE = sqrt(vcov(clustered)[1, 1]),
        generating_design_effect = 1 + (m - 1) * rho,
        fitted_design_effect = 1 + (m - 1) * icc))
print(isSingular(clustered))
print(packageVersion('lme4'))
sessionInfo()

# Numerical and structural checks for the controlled example.
stopifnot(nrow(dat) == 4000L, nlevels(dat$cluster) == 80L,
          abs(coef(naive)[1] - fixef(clustered)[1]) < 1e-8,
          sqrt(vcov(clustered)[1, 1]) > sqrt(vcov(naive)[1, 1]),
          !isSingular(clustered), clustered@optinfo$conv$opt == 0)