draft of ppc_calibration plots

Author: TeemuSailynoja

R/ppc-calibration.R

@@ -1,4 +1,4 @@
-#' PPC calibration
+# x' PPC calibration
 #'
 #' Assess the calibration of the predictive distributions `yrep` in relation to
 #' the data `y'.
@@ -14,49 +14,119 @@
 #'
 #' @section Plot Descriptions:
 #' \describe{
+#'   \item{`ppc_calibration()`,`ppc_calibration_grouped()`}{
+#'
+#'   },
 #'   \item{`ppc_calibration_overlay()`,`ppc_calibration_overlay_grouped()`}{
 #'
+#'   },
+#'   \item{`ppc_loo_calibration()`,`ppc_loo_calibration_grouped()`}{
+#'
 #'   }
 #' }
 #'
 NULL
 
+
 #' @rdname PPC-calibration
 #' @export
-.ppc_calibration_data <- function(y, prep, group = NULL) {
-  y <- validate_y(y)
-  n_obs <- length(y)
-  prep <- validate_predictions(prep, n_obs)
-  if (any(prep > 1 | prep < 0)) {
-    stop("Values of ´prep´ should be predictive probabilities between 0 and 1.")
-  }
-  if (!is.null(group)) {
-    group <- validate_group(group, n_obs)
-  } else {
-    group <- rep(1, n_obs * nrow(prep))
-  }
+ppc_calibration_overlay <- function(
+    y, prep, ..., linewidth = 0.25, alpha = 0.5) {
+  check_ignored_arguments(...)
+  data <- .ppc_calibration_data(y, prep)
+  ggplot(data) +
+    geom_abline(color = "black", linetype = 2) +
+    geom_line(
+      aes(value, cep, group = rep_id, color = "yrep"),
+      linewidth = linewidth, alpha = alpha
+    ) +
+    scale_color_ppc() +
+    bayesplot_theme_get() +
+    legend_none() +
+    coord_equal(xlim = c(0, 1), ylim = c(0, 1), expand = FALSE) +
+    xlab("Predicted probability") +
+    ylab("Conditional event probability") +
+    NULL
+}
 
-  if (requireNamespace("monotone", quietly = TRUE)) {
-    monotone <- monotone::monotone
-  } else {
-    monotone <- function(y) {
-      stats::isoreg(y)$yf
-    }
-  }
-  .ppd_data(prep, group = group) %>%
-    group_by(group, rep_id) %>%
-    mutate(
-      ord = order(value),
-      value = value[ord],
-      cep = monotone(y[ord])
-    ) |>
-    ungroup()
+#' @rdname PPC-calibration
+#' @export
+ppc_calibration_overlay_grouped <- function(
+    y, prep, group, ..., linewidth = 0.25, alpha = 0.7) {
+  check_ignored_arguments(...)
+  data <- .ppc_calibration_data(y, prep, group)
+  ggplot(data) +
+    geom_abline(color = "black", linetype = 2) +
+    geom_line(aes(value, cep, group = rep_id, color = "yrep"),
+      linewidth = linewidth, alpha = alpha
+    ) +
+    facet_wrap(vars(group)) +
+    scale_color_ppc() +
+    bayesplot_theme_get() +
+    legend_none() +
+    coord_equal(xlim = c(0, 1), ylim = c(0, 1), expand = FALSE) +
+    xlab("Predicted probability") +
+    ylab("Conditional event probability") +
+    NULL
 }
 
 #' @rdname PPC-calibration
 #' @export
-ppc_calibration_overlay <- function(
-    y, prep, ..., linewidth = 0.25, alpha = 0.7) {
+ppc_calibration <- function(
+    y, prep, prob = .95, show_mean = TRUE, ..., linewidth = 0.25, alpha = 0.7) {
+  check_ignored_arguments(...)
+  data <- .ppc_calibration_data(y, prep) %>%
+    group_by(y_id) %>%
+    summarise(
+      value = median(value),
+      lb = quantile(cep, .5 - .5 * prob),
+      ub = quantile(cep, .5 + .5 * prob),
+      cep = median(cep)
+    )
+
+  ggplot(data) +
+    aes(value, cep) +
+    geom_abline(color = "black", linetype = 2) +
+    geom_ribbon(aes(ymin = lb, ymax = ub, fill = "yrep"), alpha = alpha) +
+    geom_line(
+      aes(color = "y"),
+      linewidth = linewidth
+    ) +
+    scale_color_ppc() +
+    scale_fill_ppc() +
+    bayesplot_theme_get() +
+    # legend_none() +
+    coord_equal(xlim = c(0, 1), ylim = c(0, 1), expand = FALSE) +
+    xlab("Predicted probability") +
+    ylab("Conditional event probability") +
+    NULL
+}
+
+#' @rdname PPC-calibration
+#' @export
+ppc_calibration_grouped <- function(
+    y, prep, show_mean, ..., linewidth = 0.25, alpha = 0.7) {
+  check_ignored_arguments(...)
+  data <- .ppc_calibration_data(y, prep, group)
+  ggplot(data) +
+    geom_abline(color = "black", linetype = 2) +
+    geom_line(aes(value, cep, group = rep_id, color = "yrep"),
+      linewidth = linewidth, alpha = alpha
+    ) +
+    facet_wrap(vars(group)) +
+    scale_color_ppc() +
+    bayesplot_theme_get() +
+    legend_none() +
+    coord_equal(xlim = c(0, 1), ylim = c(0, 1), expand = FALSE) +
+    xlab("Predicted probability") +
+    ylab("Conditional event probability") +
+    NULL
+}
+
+#' @rdname PPC-calibration
+#' @export
+ppc_loo_calibration <- function(
+    y, prep, lw, ..., linewidth = 0.25, alpha = 0.7) {
   check_ignored_arguments(...)
   data <- .ppc_calibration_data(y, prep)
   ggplot(data) +
@@ -76,8 +146,8 @@ ppc_calibration_overlay <- function(
 
 #' @rdname PPC-calibration
 #' @export
-ppc_calibration_overlay_grouped <- function(
-    y, prep, group, ..., linewidth = 0.25, alpha = 0.7) {
+ppc_loo_calibration_grouped <- function(
+    y, prep, lw, ..., linewidth = 0.25, alpha = 0.7) {
   check_ignored_arguments(...)
   data <- .ppc_calibration_data(y, prep, group)
   ggplot(data) +
@@ -94,3 +164,39 @@ ppc_calibration_overlay_grouped <- function(
     ylab("Conditional event probability") +
     NULL
 }
+
+.ppc_calibration_data <- function(y, prep, group = NULL) {
+  y <- validate_y(y)
+  n_obs <- length(y)
+  prep <- validate_predictions(prep, n_obs)
+  if (any(prep > 1 | prep < 0)) {
+    stop("Values of ´prep´ should be predictive probabilities between 0 and 1.")
+  }
+  if (!is.null(group)) {
+    group <- validate_group(group, n_obs)
+  } else {
+    group <- rep(1, n_obs * nrow(prep))
+  }
+
+  if (requireNamespace("monotone", quietly = TRUE)) {
+    monotone <- monotone::monotone
+  } else {
+    monotone <- function(y) {
+      stats::isoreg(y)$yf
+    }
+  }
+  .ppd_data(prep, group = group) %>%
+    group_by(group, rep_id) %>%
+    mutate(
+      ord = order(value),
+      value = value[ord],
+      cep = monotone(y[ord])
+    ) |>
+    ungroup()
+}
+
+.loo_resample_data <- function(prep, lw, psis_object) {
+  lw <- .get_lw(lw, psis_object)
+  stopifnot(identical(dim(prep), dim(lw)))
+  # Work in progress here...
+}