Fix matchindex() tie breaking

NEWS.md

@@ -1,3 +1,7 @@
+# mice (development version)
+
+- Fixes miscalibrated donor selection in `matchindex()`, the default PMM matching engine since mice 3.11.8 (2020), for predictive mean matching models where many target cases share the same (or nearly the same) predicted value, e.g. an intercept-only imputation model, models with few categorical predictors, or heavily rounded/discretised outcomes. Previously, ties among donors were broken once per call to `matchindex()` (i.e. once per imputed dataset), so all target cases with a tied predicted value drew from the same small, fixed subset of donors within one completed dataset. This inflated the between-imputation variance and produced overly conservative (too wide) confidence intervals under Rubin's rules. Ties are now broken independently for each target case. **Affected versions:** mice 3.11.8 through 3.19.8 (any analysis using the default `method = "pmm"` with tied or near-tied predicted values). To reproduce results computed with an affected version, reinstall that specific version, e.g. `remotes::install_version("mice", version = "3.19.8")`, rather than relying on a compatibility flag in the current version. See [#757](https://github.com/amices/mice/issues/757) and the original discussion in [#236](https://github.com/amices/mice/issues/236).
+
 # mice 3.19.8
 
 - Implements new behavior in `mice.mids() that inherits the imputation methods from the trained imputation model rather than re-initializing methods anew. This gives consistent behavior between train and deployment imputation models ([#753](https://github.com/amices/mice/issues/753))

R/RcppExports.R

@@ -24,23 +24,34 @@ matcher <- function(obs, mis, k) {
 #' neighbours in \code{d}, randomly draws one of these neighbours, and
 #' returns its position in vector \code{d}.
 #'
-#' Fast predictive mean matching algorithm in seven steps:
+#' Fast predictive mean matching algorithm in eight steps:
 #'
 #' 1. Shuffle records to remove effects of ties
 #'
 #' 2. Obtain sorting order on shuffled data
 #'
 #' 3. Calculate index on input data and sort it
 #'
-#' 4. Pre-sample vector \code{h} with values between 1 and \code{k}
+#' 4. Identify contiguous tie blocks (maximal runs of equal values) in
+#' the sorted donor data
+#'
+#' 5. Pre-sample vector \code{h} with values between 1 and \code{k}
 #'
 #' For each of the \code{n0} elements in \code{t}:
 #'
-#'   5. find the two adjacent neighbours
+#'   6. find the two adjacent neighbours
 #'
-#'   6. find the \code{h_i}'th nearest neighbour
+#'   7. find the \code{h_i}'th nearest neighbour, breaking ties among
+#'   equidistant donors with an independent random draw for this
+#'   target case (not a single order shared by all target cases).
+#'   If the tie block nearest to the target already has at least
+#'   \code{k} members, all \code{k} nearest neighbours necessarily lie
+#'   within that one block, so \code{h_i} is skipped and a donor is
+#'   drawn directly and uniformly from the block; otherwise, the
+#'   neighbours are found by stepping outward from the target,
+#'   drawing without replacement from each tie block encountered
 #'
-#'   7. store the index of that neighbour
+#'   8. store the index of that neighbour
 #'
 #' Return vector of \code{n0} positions in \code{d}.
 #'
@@ -54,6 +65,16 @@ matcher <- function(obs, mis, k) {
 #' See \url{https://github.com/amices/mice/issues/236}.
 #' This function is a replacement for the \code{matcher()} function that has
 #' been in default in \code{mice} since version \code{2.22} (June 2014).
+#'
+#' Prior to the fix in \url{https://github.com/amices/mice/issues/757}, ties
+#' among donors were broken once per call to \code{matchindex()} (i.e. once
+#' per imputed dataset), so all target cases sharing the same predicted
+#' value drew from the same fixed subset of donors within one completed
+#' dataset. This inflated between-imputation variance and produced overly
+#' conservative confidence intervals under Rubin's rules. This was most
+#' visible in intercept-only imputation models, models with few categorical
+#' predictors, or heavily rounded/discretised outcomes. Ties are now broken
+#' independently for each target case.
 #' @examples
 #' set.seed(1)
 #'

R/mice.impute.pmm.R

@@ -43,6 +43,10 @@
 #' \code{2.22} (June 2014) to \code{3.11.7} (Oct 2020). Since version \code{3.12.0}
 #' \code{mice()} uses the much faster \code{matchindex} C function. Use
 #' the deprecated \code{matcher} function only for exact reproduction.
+#' Note that \code{matcher()} itself is miscalibrated when donors are tied
+#' (see \url{https://github.com/amices/mice/issues/757}), so it should not
+#' be used as a "safe" alternative for tied data, only for reproducing
+#' results computed with an older \code{mice} version.
 #' @param \dots Other named arguments.
 #' @return Vector with imputed data, same type as \code{y}, and of length
 #' \code{sum(wy)}

src/matchindex.cpp

@@ -19,23 +19,34 @@ using namespace Rcpp;
 //' neighbours in \code{d}, randomly draws one of these neighbours, and
 //' returns its position in vector \code{d}.
 //'
-//' Fast predictive mean matching algorithm in seven steps:
+//' Fast predictive mean matching algorithm in eight steps:
 //'
 //' 1. Shuffle records to remove effects of ties
 //'
 //' 2. Obtain sorting order on shuffled data
 //'
 //' 3. Calculate index on input data and sort it
 //'
-//' 4. Pre-sample vector \code{h} with values between 1 and \code{k}
+//' 4. Identify contiguous tie blocks (maximal runs of equal values) in
+//' the sorted donor data
+//'
+//' 5. Pre-sample vector \code{h} with values between 1 and \code{k}
 //'
 //' For each of the \code{n0} elements in \code{t}:
 //'
-//'   5. find the two adjacent neighbours
+//'   6. find the two adjacent neighbours
 //'
-//'   6. find the \code{h_i}'th nearest neighbour
+//'   7. find the \code{h_i}'th nearest neighbour, breaking ties among
+//'   equidistant donors with an independent random draw for this
+//'   target case (not a single order shared by all target cases).
+//'   If the tie block nearest to the target already has at least
+//'   \code{k} members, all \code{k} nearest neighbours necessarily lie
+//'   within that one block, so \code{h_i} is skipped and a donor is
+//'   drawn directly and uniformly from the block; otherwise, the
+//'   neighbours are found by stepping outward from the target,
+//'   drawing without replacement from each tie block encountered
 //'
-//'   7. store the index of that neighbour
+//'   8. store the index of that neighbour
 //'
 //' Return vector of \code{n0} positions in \code{d}.
 //'
@@ -49,6 +60,16 @@ using namespace Rcpp;
 //' See \url{https://github.com/amices/mice/issues/236}.
 //' This function is a replacement for the \code{matcher()} function that has
 //' been in default in \code{mice} since version \code{2.22} (June 2014).
+//'
+//' Prior to the fix in \url{https://github.com/amices/mice/issues/757}, ties
+//' among donors were broken once per call to \code{matchindex()} (i.e. once
+//' per imputed dataset), so all target cases sharing the same predicted
+//' value drew from the same fixed subset of donors within one completed
+//' dataset. This inflated between-imputation variance and produced overly
+//' conservative confidence intervals under Rubin's rules. This was most
+//' visible in intercept-only imputation models, models with few categorical
+//' predictors, or heavily rounded/discretised outcomes. Ties are now broken
+//' independently for each target case.
 //' @examples
 //' set.seed(1)
 //'
@@ -118,6 +139,24 @@ IntegerVector matchindex(NumericVector d, NumericVector t, int k = 5) {
     ysort[i] = d(id(i));
   }
 
+  // 3b. Identify tie blocks: blockOf[i] is the id of the maximal run of
+  // equal ysort values that position i belongs to. Donors tied at the same
+  // distance from a target must be selected independently at random for
+  // each target, not in the fixed order left by the single pre-shuffle
+  // above (see https://github.com/amices/mice/issues/236).
+  std::vector<int> blockOf(n1);
+  std::vector<int> blockStart;
+  std::vector<int> blockEnd;
+  for (int i = 0; i < n1; i++) {
+    if (i == 0 || ysort[i] != ysort[i - 1]) {
+      blockStart.push_back(i);
+      blockEnd.push_back(i);
+    } else {
+      blockEnd.back() = i;
+    }
+    blockOf[i] = static_cast<int>(blockStart.size()) - 1;
+  }
+
   // 4. Pre-sample n0 values between 1 and k
   // restrict 1 <= k <= n1
   k = (k <= n1) ? k : n1;
@@ -128,10 +167,47 @@ IntegerVector matchindex(NumericVector d, NumericVector t, int k = 5) {
 
   IntegerVector idx(n0);
 
+  // Scratch space reused across target cases: tracks, within a tie block
+  // touched by the current target, which sorted positions have already
+  // been consumed by that target.
+  std::vector<char> consumed(n1, 0);
+  std::vector<int> touched;
+  touched.reserve(k);
+
+  // Draws a not-yet-consumed donor from the tie block containing sorted
+  // position `pos`, chosen uniformly at random among the remaining
+  // (unconsumed) members of that block for this target only.
+  auto consumeFromBlock = [&](int pos) -> int {
+    int b = blockOf[pos];
+    int lo = blockStart[b];
+    int hi = blockEnd[b];
+    int blockLen = hi - lo + 1;
+    if (blockLen == 1) {
+      consumed[pos] = 1;
+      touched.push_back(pos);
+      return pos;
+    }
+    // Reservoir-style pick among unconsumed members of [lo, hi]
+    int remaining = 0;
+    for (int j = lo; j <= hi; j++) if (!consumed[j]) remaining++;
+    int pick = (remaining > 1) ? (int) std::floor(unif_rand() * remaining) : 0;
+    int seen = 0;
+    int chosen = lo;
+    for (int j = lo; j <= hi; j++) {
+      if (!consumed[j]) {
+        if (seen == pick) { chosen = j; break; }
+        seen++;
+      }
+    }
+    consumed[chosen] = 1;
+    touched.push_back(chosen);
+    return chosen;
+  };
+
   // loop over the target units
   for (int i = 0; i < n0; i++) {
     double val = t(i);
-    int hi = h(i);
+    int hi_k = h(i);
     int count = 0;
 
     // 5. find the two adjacent neighbours
@@ -140,36 +216,65 @@ IntegerVector matchindex(NumericVector d, NumericVector t, int k = 5) {
     int r = iter - ysort.begin();
     int l = r - 1;
 
+    // Fast path: if the nearest side's tie block alone already has at
+    // least k members, all k nearest neighbours necessarily come from
+    // that single block, so the h_i'th draw is uniform over the whole
+    // block regardless of h_i. We can skip the walk and h entirely and
+    // just draw one donor directly from the block.
+    if (l >= 0 || r < n1) {
+      bool takeLeft = (r >= n1) || (l >= 0 && val - ysort[l] < ysort[r] - val);
+      int nearPos = takeLeft ? l : r;
+      int b = blockOf[nearPos];
+      int lo = blockStart[b];
+      int hi = blockEnd[b];
+      int blockLen = hi - lo + 1;
+      // Only safe when the block does not extend past the array bound
+      // being approached (it never does, by construction of blockStart/
+      // blockEnd), and contains at least k members.
+      if (blockLen >= hi_k) {
+        int pick = lo + (int) std::floor(unif_rand() * blockLen);
+        idx(i) = id[pick];
+        continue;
+      }
+    }
+
     // 6. find the h_i'th nearest neighbour
     // 7. store the index of that neighbour
 
     // Compare elements on left and right of crossover
     // point to find the h'th closest match
     // Inspired on Polkas: https://github.com/Polkas/miceFast/issues/10
-    while (count < hi && l >= 0 && r < n1)
+    // Ties (equal ysort values) are resolved by an independent random
+    // draw per target case, see consumeFromBlock() above.
+    while (count < hi_k && l >= 0 && r < n1)
     {
       if (val - ysort[l] < ysort[r] - val)
       {
-        idx(i) = id[l--];
+        idx(i) = id[consumeFromBlock(l--)];
       } else {
-        idx(i) = id[r++];
+        idx(i) = id[consumeFromBlock(r++)];
       }
       count++;
     }
 
     // If right side is exhausted, take left elements
-    while (count < hi && l >= 0)
+    while (count < hi_k && l >= 0)
     {
-      idx(i) = id[l--];
+      idx(i) = id[consumeFromBlock(l--)];
       count++;
     }
 
     // If left side is exhausted, take right elements
-    while (count < hi && r < n1)
+    while (count < hi_k && r < n1)
     {
-      idx(i) = id[r++];
+      idx(i) = id[consumeFromBlock(r++)];
       count++;
     }
+
+    // Reset only the positions touched by this target, so the next
+    // target's tie-breaking is independent again.
+    for (size_t j = 0; j < touched.size(); j++) consumed[touched[j]] = 0;
+    touched.clear();
   }
 
   return idx + 1;