## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  # Evaluated wherever INLA and the packages used here are installed.
  eval     = all(vapply(
    c("INLA", "scoringRules", "dplyr", "ggplot2", "broom"),
    requireNamespace, logical(1), quietly = TRUE
  ))
)


## ----message = FALSE----------------------------------------------------------
library(inlabru)
library(inlaops)
library(scoringRules)
library(broom)
library(dplyr)
library(ggplot2)


## -----------------------------------------------------------------------------
set.seed(42)
n   <- 500
dat <- data.frame(
  x1 = rnorm(n),
  x2 = rnorm(n),
  x3 = rnorm(n)
)
dat$y <- 1.5 + 0.9 * dat$x1 + 0.3 * dat$x2 + rnorm(n, sd = 0.8)

# Hold out 20% for evaluation
idx_holdout <- sample(n, size = 100L)
train_dat   <- dat[-idx_holdout, ]
test_dat    <- dat[ idx_holdout, ]


## -----------------------------------------------------------------------------
fit_full <- bru(
  y ~ Intercept(1) + x1 + x2,
  family  = "gaussian",
  data    = train_dat,
  options = bru_options(control.compute = list(config = TRUE, dic = TRUE, waic = TRUE))
)

fit_reduced <- bru(
  y ~ Intercept(1) + x1,
  family  = "gaussian",
  data    = train_dat,
  options = bru_options(control.compute = list(config = TRUE, dic = TRUE, waic = TRUE))
)


## -----------------------------------------------------------------------------
scores_full <- bru_score(
  fit_full,
  newdata      = test_dat,
  pred_formula = ~ rnorm(
    length(x1), mean = Intercept + x1 + x2,
    sd   = 1 / sqrt(Precision_for_the_Gaussian_observations)
  ),
  y_col        = "y",
  n_samples    = 2000L,
  seed         = 1L
)

scores_reduced <- bru_score(
  fit_reduced,
  newdata      = test_dat,
  pred_formula = ~ rnorm(
    length(x1), mean = Intercept + x1,
    sd   = 1 / sqrt(Precision_for_the_Gaussian_observations)
  ),
  y_col        = "y",
  n_samples    = 2000L,
  seed         = 1L
)


## -----------------------------------------------------------------------------
bind_rows(
  mutate(scores_full,    model = "full"),
  mutate(scores_reduced, model = "reduced")
) |>
  group_by(model) |>
  summarise(
    mean_crps = mean(crps),
    mean_logs = mean(logs),
    mean_rmse = sqrt(mean(se))
  )


## -----------------------------------------------------------------------------
bind_rows(
  mutate(scores_full,    model = "full"),
  mutate(scores_reduced, model = "reduced")
) |>
  ggplot(aes(x = crps, fill = model)) +
  geom_density(alpha = 0.5) +
  labs(x = "CRPS (lower is better)", y = "Density") +
  theme_minimal()


## -----------------------------------------------------------------------------
glance(fit_full)