How to calculate mean absolute error for Bayesian modeling based on Leave-One-Out Cross-Validation with Brms package in R

Question

I am fitting a Bayesian model to predict a test score using the Brms package. I would like to know how to calculate the 'Mean Absolute Error' based on 'Leave-One-Out Cross-Validation' (LOO) using the LOO package, but I could not find any information related to how to actualize it by myself.

I would really appreciate if somebody demonstrate me how to calculate MAE based on LOO.

Here is a replicable sample code:

set.seed(123) # for reproducibility
n <- 100 # number of observations
predictor_1 <- rnorm(n)
predictor_2 <- rnorm(n)
test_score <- 5 + 2*predictor_1 + 3*predictor_2 + rnorm(n)
data <- data.frame(test_score, predictor_1, predictor_2)
head(data)

fit <- brm(test_score ~ predictor_1 + predictor_2, data = data)
predicted_test_score <- predict(fit)

# calculate mean absolute error
mae <- mean(abs(predicted_test_score - data$test_score))
mae

Answer 1

The function brms::kfold_predict() can be used for this and there is an easily adaptable example in the help file, see ?brms::kfold_predict().

Perform LOOCV and save the fits (object size grows quickly with nobs!)

fit_LOO <- kfold(fit, save_fits = TRUE, chains = 1, folds = "loo")

Obtain samples (1000 is the default) from the posterior predictive distribution of each of the 100 observations. This yields a 1000x100 matrix.

predict_LOO <- kfold_predict(fit_LOO, method = "predict")

> dim(predict_LOO$yrep)
[1] 1000  100

Define and compute your desired loss (based on posterior means):

MAE <- function(y, yrep) {
  yrep_mean <- colMeans(yrep)
  mean(abs(yrep_mean - y))
}

> MAE(y = predict_LOO$y, yrep = predict_LOO$yrep)
[1] 0.7846185

Answer 2

For you to calculate the MAE based on LOO using the loo package, you need to fit the Bayesian model using the brm function from the brms package, then, perform Leave-One-Out Cross-Validation using the loo function. Also, you will need to extract the pointwise log-likelihood matrix from the brmsfit object and calculate the expected log pointwise predictive density (ELPD) differences using the loo::loo_compare function.Finally, you have to calculate the LOO mean absolute error.

if (!requireNamespace("brms", quietly = TRUE)) {
  install.packages("brms")
}
library(brms)
if (!requireNamespace("loo", quietly = TRUE)) {
install.packages("loo")
}
library(loo)
fit <- brm(test_score ~ predictor_1 + predictor_2, data = data, save_all_pars = TRUE)
loo_fit <- loo(fit, save_psis = TRUE)
log_lik <- log_lik(fit)
elpd_diff <- loo::loo_compare(loo_fit, loo_fit)
loo_mae <- mean(abs(exp(elpd_diff$diff)))
print(loo_mae)