fix(calibration,sensitivity): anchored-curve weights, pooling, newdata support

R/anchor_curve.R

@@ -39,7 +39,45 @@ anchor_draws = function(fit) {
   family = spec$prior$interaction_prior_type %||% "normal"
   list(
     theta = theta, gamma = gamma, indicator = gamma,
-    family = tolower(family)
+    family = tolower(family),
+    diagonal = anchor_diagonal(spec, raw)
+  )
+}
+
+
+# ------------------------------------------------------------------
+# anchor_diagonal
+# ------------------------------------------------------------------
+# The Gamma prior on the precision diagonal, in the frame the weight ratio
+# needs it, or NULL for a model that has no such diagonal.
+#
+# Under vary = "slab-and-diagonal" the refits set the diagonal rate to eta / s
+# (vary_diagonal_rate()), so the weight ratio carries the diagonal prior as
+# well as the slab. The sampler evaluates that prior at K_jj / 2
+# (ggm_model.cpp:320, mixed_mrf_gradient.cpp:559); the raw main draws store
+# K_jj.
+#
+# @param spec  The fit's spec.
+# @param raw   Its raw samples.
+#
+# Returns: list(sum = per-chain vectors of sum_j K_jj / 2, n = number of
+#   diagonal elements, shape = the Gamma shape alpha), or NULL.
+# ------------------------------------------------------------------
+anchor_diagonal = function(spec, raw) {
+  cols = if(identical(spec$model_type, "ggm")) {
+    seq_len(ncol(raw$main[[1]]))
+  } else if(identical(spec$model_type, "mixed_mrf")) {
+    which(endsWith(raw$parameter_names[["main"]], "(precision diag)"))
+  } else {
+    integer(0)
+  }
+  if(length(cols) == 0L) {
+    return(NULL)
+  }
+  list(
+    sum = lapply(raw$main, function(m) 0.5 * rowSums(m[, cols, drop = FALSE])),
+    n = length(cols),
+    shape = spec$prior$scale_shape %||% 1
   )
 }
 
@@ -65,15 +103,20 @@ anchor_draws = function(fit) {
 compare_anchor_draws = function(fit, spec, raw) {
   arguments = extract_arguments(fit)
   names_all = raw$parameter_names
-  num_variables = as.integer(arguments$num_variables %||% arguments$no_variables)
+  num_variables = as.integer(
+    arguments[["num_variables"]] %||% arguments[["no_variables"]]
+  )
 
   num_pairs = length(names_all$pairwise_baseline)
   num_contrasts = length(names_all$pairwise_diff) / num_pairs
   num_main_baseline = length(names_all$main_baseline)
 
   # Main parameters per variable: one per threshold for an ordinal variable,
-  # a linear and a quadratic term for a Blume-Capel one.
-  block = ifelse(arguments$is_ordinal, as.integer(arguments$num_categories), 2L)
+  # a linear and a quadratic term for a Blume-Capel one. The field is named
+  # per model type -- is_ordinal_variable on a compare fit, is_ordinal on the
+  # single-network ones -- and is read exactly, not by partial match.
+  is_ordinal = arguments[["is_ordinal_variable"]] %||% arguments[["is_ordinal"]]
+  block = ifelse(is_ordinal, as.integer(arguments[["num_categories"]]), 2L)
 
   # Indicators run over the upper triangle with the diagonal: (v, v) is that
   # variable's main-effect difference, (i, j) the pair's.
@@ -115,15 +158,41 @@ compare_anchor_draws = function(fit, spec, raw) {
 # scale, for one chain. Only included edges contribute; the weight is
 # the slab-density ratio summed over them.
 #
+# Under vary = "slab-and-diagonal" the refits also move the Gamma prior on the
+# precision diagonal, rate = eta / s, so the target differs from the anchor
+# there too and the weight carries that ratio as well. For Gamma(alpha,
+# eta / s) on each of the n diagonal values x_j = K_jj / 2, the per-draw term is
+#
+#   n * alpha * log((eta / s) / (eta / s_a)) - (eta / s - eta / s_a) * S_t
+#     = -n * alpha * log(s / s_a) - eta * (1 / s - 1 / s_a) * S_t,
+#
+# with S_t = sum_j K_jj / 2 at draw t. The leading term is the same for every
+# draw and cancels in the self-normalized average; the S_t term reweights.
+#
 # @param theta    iter x edges matrix, slab frame.
 # @param gamma    iter x edges indicator matrix.
 # @param family   "normal" or "cauchy".
 # @param s_a      Anchor slab scale.
 # @param s_grid   Target slab scales.
+# @param diagonal The anchor's anchor_diagonal() record, or NULL; `sum` is
+#                 this chain's per-draw S_t.
+# @param eta      Standardized diagonal rate held fixed by the sweep, or NULL
+#                 when the sweep leaves the diagonal prior alone.
 #
 # Returns: iter x length(s_grid) matrix of log weights.
 # ------------------------------------------------------------------
-anchor_log_weights = function(theta, gamma, family, s_a, s_grid) {
+anchor_log_weights = function(theta, gamma, family, s_a, s_grid,
+                              diagonal = NULL, eta = NULL) {
+  # The reweighting identity is written per slab family; a scale-free family
+  # (beta-prime) has no s to sweep and any other one has a different density
+  # ratio.
+  if(!family %in% c("normal", "cauchy")) {
+    stop(
+      "The anchored sensitivity curve supports a normal or Cauchy slab on ",
+      "the swept parameters; this fit uses a '", family, "' slab, for which ",
+      "there is no slab scale to reweight across."
+    )
+  }
   incl = gamma == 1
   m_t = rowSums(incl)
   if(identical(family, "normal")) {
@@ -139,6 +208,16 @@ anchor_log_weights = function(theta, gamma, family, s_a, s_grid) {
     }, numeric(length(m_t)))
   }
   if(is.null(dim(lw))) lw = matrix(lw, nrow = 1L)
+  # Only when the sweep moved the diagonal prior.
+  if(!is.null(eta) && !is.null(diagonal)) {
+    s_t = diagonal$sum
+    dlw = vapply(s_grid, function(s) {
+      -diagonal$n * diagonal$shape * log(s / s_a) -
+        eta * (1 / s - 1 / s_a) * s_t
+    }, numeric(length(s_t)))
+    if(is.null(dim(dlw))) dlw = matrix(dlw, nrow = 1L)
+    lw = lw + dlw
+  }
   lw
 }
 
@@ -152,14 +231,22 @@ anchor_log_weights = function(theta, gamma, family, s_a, s_grid) {
 # @param draws    Output of anchor_draws() for the anchor fit.
 # @param s_a      Anchor slab scale.
 # @param s_grid   Target slab scales (length P).
+# @param eta      Standardized diagonal rate the sweep holds fixed, or NULL
+#                 when the sweep leaves the precision diagonal alone. The
+#                 reported importance ESS follows from the same weights.
 #
 # Returns: list(pip = P x E pooled, ess = length-P pooled importance ESS,
 #   chain_pip = list of P x E per chain, chain_ess = C x P).
 # ------------------------------------------------------------------
-anchor_reweight = function(draws, s_a, s_grid) {
+anchor_reweight = function(draws, s_a, s_grid, eta = NULL) {
   n_chain = length(draws$theta)
+  diag = draws$diagonal
   lw = lapply(seq_len(n_chain), function(c) {
-    anchor_log_weights(draws$theta[[c]], draws$gamma[[c]], draws$family, s_a, s_grid)
+    anchor_log_weights(
+      draws$theta[[c]], draws$gamma[[c]], draws$family, s_a, s_grid,
+      diagonal = if(is.null(diag)) NULL else list(sum = diag$sum[[c]], n = diag$n, shape = diag$shape),
+      eta = eta
+    )
   })
   # The weight sums over gated parameters; the inclusion probability is
   # reported on the indicators, which are the same thing for bgm() and a
@@ -239,9 +326,14 @@ assemble_curve = function(reweights, usable, ess_floor, anchor_index = NULL) {
     cnum = replicate(n_chain, numeric(n_edge), simplify = FALSE)
     for(a in contrib) {
       pa = reweights[[a]]$pip[p, ]
-      # Inverse-variance weight on the PIP scale (var ~ p(1-p)/ESS), clamped
-      # so a saturated edge does not divide by zero.
-      wa = ess_mat[p, a] / pmax(pa * (1 - pa), 1e-6)
+      ea = ess_mat[p, a]
+      # Inverse-variance weight on the PIP scale (var ~ p(1-p)/ESS), with the
+      # variance computed on the add-half smoothed proportion
+      # (ess * pa + 0.5) / (ess + 1), so a saturated anchor (PIP 0 or 1) gets
+      # the variance its ESS supports. The pooled estimate uses the unsmoothed
+      # pa; the 1e-6 floor stays as a backstop.
+      pa_s = (ea * pa + 0.5) / (ea + 1)
+      wa = ea / pmax(pa_s * (1 - pa_s), 1e-6)
       wsum = wsum + wa
       pnum = pnum + wa * pa
       for(cc in seq_len(n_chain)) {

R/calibration_check.R

@@ -89,14 +89,21 @@ fitted_observed_data = function(bgms_object) {
 # 1-based index of each observed value among a discrete variable's categories,
 # i.e. the column of the predicted-probability matrix the case fell in.
 #
+# A value the fit has no category for has nothing to calibrate against. The
+# recode-map lookup returns NA for it and the Blume-Capel branch returns an
+# index outside the category range; both are rejected here.
+#
 # @param values      Observed values on the original scale.
 # @param levels_v    The variable's recode map, or NULL for Blume-Capel.
 # @param shift_v     The variable's additive shift, or NA for regular ordinal.
 # @param variable    Variable name, for the error message.
+# @param num_categories  The variable's category count, or NULL to check only
+#                    that every value was located at all.
 #
 # Returns: integer vector of category indices.
 # ------------------------------------------------------------------
-discrete_category_index = function(values, levels_v, shift_v, variable) {
+discrete_category_index = function(values, levels_v, shift_v, variable,
+                                   num_categories = NULL) {
   index = if(!is.null(levels_v)) {
     if(!is.null(names(levels_v))) {
       # bgmCompare stores a named lookup from the original value to the final
@@ -117,6 +124,31 @@ discrete_category_index = function(values, levels_v, shift_v, variable) {
       "located. Re-fit with the current bgms version."
     )
   }
+  outside = if(is.null(num_categories)) {
+    is.na(index)
+  } else {
+    is.na(index) | index < 1L | index > as.integer(num_categories)
+  }
+  if(any(outside)) {
+    offending = sort(unique(values[outside]))
+    shown = utils::head(offending, 5L)
+    more = if(length(offending) > length(shown)) {
+      sprintf(" (and %d more)", length(offending) - length(shown))
+    } else {
+      ""
+    }
+    stop(
+      "Variable '", variable, "' carries value", if(length(offending) == 1L) {
+        " "
+      } else {
+        "s "
+      },
+      paste(format(shown, trim = TRUE), collapse = ", "), more,
+      " that the fit has no category for, so there is nothing to calibrate ",
+      "them against. Supply 'newdata' on the categories the model was fitted ",
+      "to."
+    )
+  }
   index
 }
 
@@ -310,7 +342,8 @@ pav_panels = function(predicted, observed, levels_list, shifts, grid, nrep, prob
     p = as.vector(cumulative[, -num_categories, drop = FALSE])
 
     observed_index = discrete_category_index(
-      observed[, v], levels_list[[v]], shifts[[v]], v
+      observed[, v], levels_list[[v]], shifts[[v]], v,
+      num_categories = num_categories
     )
     y = as.integer(outer(observed_index, seq_len(num_thresholds), "<="))
     observed_curve = pav_curve(p, y, grid)
@@ -415,8 +448,10 @@ calibration_result = function(panels, nrep, probs, grid, ndraws) {
 #' parameter uncertainty sits inside the distribution the observation is
 #' transformed by; the two panel kinds therefore differ in what they condition
 #' on, the isotonic one on posterior-mean probabilities and this one on the full
-#' predictive distribution. The band is the simultaneous envelope of the
-#' empirical distribution functions of \eqn{n} independent uniforms. It is not
+#' predictive distribution. The band is pointwise: at each grid point it spans
+#' the central interval of the simulated uniform ECDFs given by `probs`, 95%
+#' by default, so a calibrated variable may still exit the band at a few
+#' isolated points. It is not
 #' resampled from the model, because it does not have to be: whatever the
 #' conditional density was, a value drawn from it transforms to an exact
 #' \eqn{\textrm{Uniform}(0, 1)}, so the null is known and one band serves every