# =====================================================================
# Toy demo: Bayesian random-intercept (hierarchical) model with R-INLA
# ---------------------------------------------------------------------
# Model:
#   y_ij = beta0 + beta1 * x_ij + u_j + eps_ij
#   u_j   ~ Normal(0, sigma_u^2)   (group-level random intercept)
#   eps_ij ~ Normal(0, sigma_e^2)  (residual noise)
# =====================================================================

# install.packages("INLA",
#   repos = c(getOption("repos"), INLA = "https://inla.r-inla-download.org/R/stable"),
#   dep = TRUE)
library(INLA)

# ---------------------------------------------------------------------
# 1. Simulate a toy dataset with KNOWN ("true") parameters
# ---------------------------------------------------------------------
set.seed(42)

n_groups   <- 30      # number of clusters (e.g., hospitals / batches / subjects)
n_per_grp  <- 20      # observations per group
N          <- n_groups * n_per_grp

# True parameter values (the targets we hope to recover)
beta0_true   <- 1.0   # intercept
beta1_true   <- 2.0   # slope for x
sigma_u_true <- 1.5   # SD of the group random intercepts
sigma_e_true <- 1.0   # residual SD

# Group IDs and a continuous covariate
group <- rep(1:n_groups, each = n_per_grp)
x     <- rnorm(N, mean = 0, sd = 1)

# Draw one random intercept per group, then expand to all rows
u_group <- rnorm(n_groups, mean = 0, sd = sigma_u_true)
u       <- u_group[group]

# Linear predictor + Gaussian noise
eps <- rnorm(N, mean = 0, sd = sigma_e_true)
y   <- beta0_true + beta1_true * x + u + eps

mydata <- data.frame(y = y, x = x, group = group)
head(mydata)

# ---------------------------------------------------------------------
# 2. Fit the random-intercept model with INLA
# ---------------------------------------------------------------------
# f(group, model = "iid") = exchangeable group-level random intercept
formula <- y ~ x + f(group, model = "iid")

result <- inla(
  formula,
  family            = "gaussian",
  data              = mydata,
  control.predictor = list(compute = TRUE),
  control.compute   = list(dic = TRUE, waic = TRUE)
)

summary(result)

# ---------------------------------------------------------------------
# 3. Compare INLA posteriors with the TRUE values
# ---------------------------------------------------------------------
# Fixed effects (beta0, beta1)
cat("\n--- Fixed effects ---\n")
print(result$summary.fixed[, c("mean", "sd", "0.025quant", "0.975quant")])
cat("True beta0 =", beta0_true, " | True beta1 =", beta1_true, "\n")

# INLA reports PRECISIONS (1/variance) for hyperparameters.
# Convert posterior precision summaries to SDs for interpretation.
prec2sd <- function(prec) 1 / sqrt(prec)

hyper <- result$summary.hyperpar
cat("\n--- Hyperparameters (precision scale) ---\n")
print(hyper[, c("mean", "0.025quant", "0.975quant")])

# Point estimates of SDs from posterior-mean precisions
sigma_e_hat <- prec2sd(hyper["Precision for the Gaussian observations", "mean"])
sigma_u_hat <- prec2sd(hyper["Precision for group", "mean"])

cat("\n--- Recovered SDs vs truth ---\n")
cat(sprintf("Residual SD : estimate = %.3f  (true = %.3f)\n", sigma_e_hat, sigma_e_true))
cat(sprintf("Group   SD  : estimate = %.3f  (true = %.3f)\n", sigma_u_hat, sigma_u_true))

# ---------------------------------------------------------------------
# 4. Recover the group-level random intercepts u_j
# ---------------------------------------------------------------------
re <- result$summary.random$group   # one row per group
plot(u_group, re$mean,
     xlab = "True random intercept u_j",
     ylab = "INLA posterior mean",
     main = "Recovery of group random intercepts",
     pch = 19)
abline(0, 1, col = "red", lwd = 2)   # y = x reference line

cat("\nCorrelation(true u_j, estimated u_j) =",
    round(cor(u_group, re$mean), 3), "\n")