Within score caliper

NEWS.md

@@ -1,6 +1,8 @@
 ## Changes in **optmatch** Version 0.10.8
 
-- Updates to internal C++ code.
+- `optmatch:::scoreCaliper()` gains an optional
+  `within=` argument (#245)
+- Updates to internal C++ code
 
 ## Changes in **optmatch** Version 0.10.7
 

R/InfinitySparseMatrix.R

@@ -1053,7 +1053,7 @@ dbind <- function(..., force_unique_names = FALSE) {
                          sum(vapply(lapply(mats, methods::slot, "dimension"), "[", 1, 2))))
 
   # This needs to be much smarter, especially if any element is already a BISM
-  groups <- as.factor(rep(0:(length(mats)-1), times =
+  groups <- as.factor(rep(seq_along(mats), times =
                                       vapply(lapply(mats, slot, "colnames"), length, 1) +
                                         vapply(lapply(mats, slot, "rownames"), length, 1)))
   names(groups) <- do.call(c, Map(c, cnameslist, rnameslist))

R/match_on.R

@@ -716,7 +716,7 @@ match_on.numeric <- function(x, within = NULL, caliper = NULL, exclude = NULL, d
             colnames = c(allowed@colnames, n0_names))
 
     } else {
-        allowed <- scoreCaliper(x, z, caliper)
+        allowed <- scoreCaliper(x, z, caliper, within)
     }
 
     if (!is.null(within)) {
@@ -750,9 +750,79 @@ match_on.numeric <- function(x, within = NULL, caliper = NULL, exclude = NULL, d
 #' @param z The treatment assignment vector (same length as \code{x})
 #' @param caliper The width of the caliper with respect to the scores
 #'   \code{x}.
+#' @param within A valid distance specification, such as the result of
+#'   \code{\link{exactMatch}} or \code{\link{caliper}}. Finite entries indicate
+#'   which distances to create. Including this argument can significantly speed
+#'   up computation for sparse matching problems. Specify this filter either via
+#'   \code{within} or via \code{strata} elements of a formula; mixing these
+#'   methods will fail.
 #' @return An \code{InfinitySparseMatrix} object, suitable to be
 #'   passed to \code{\link{match_on}} as an \code{within} argument.
-scoreCaliper <- function(x, z, caliper) {
+scoreCaliper <- function(x, z, caliper, within=NULL) {
+  if (!is.null(within)) {
+    allowed_list <- list()
+    for (blk in levels(within@groups)) {
+      ix <- within@groups == blk
+        if (sum(ix & z) > 0 & sum(ix & !z) > 0) {
+          allowed_list[[blk]] <- scoreCaliperBlock(x[ix], z[ix], caliper)
+        }
+    }
+    allowed <- dbind(allowed_list)
+    allowed <- addEligibleTreatments(allowed,
+                  within@rownames[!within@rownames %in% allowed@rownames], Inf)
+    allowed <- addIneligibleControls(allowed,
+                  within@colnames[!within@colnames %in% allowed@colnames], Inf)
+  } else {
+    allowed <- scoreCaliperBlock(x, z, caliper)
+  }
+  return(allowed)
+}
+
+##' Helper function to add rows to an existing ISM
+##'
+##' Adds rows of value \code{val} for each entry in \code{names}.
+##' @param ism An ISM.
+##' @param names A vector of names to be added to the rows.
+##' @param val The value to be added.
+##' @return The ISM with any padded rows.
+##' @keywords internal
+addEligibleTreatments <- function(ism, names, val=0) {
+  numnames <- length(names)
+  if (numnames == 0) return(ism) # Short circuit
+  ismdim <- ism@dimension
+  ism@.Data <- c(ism@.Data, rep(val, numnames*ismdim[2]))
+  ism@rows <- c(ism@rows,
+                rep(seq(ismdim[1] + 1, ismdim[1] + numnames),
+                    each=ismdim[2]))
+  ism@cols <- c(ism@cols, rep(1:ismdim[2], times=numnames))
+  ism@dimension <- c(ismdim[1] + as.integer(numnames), ismdim[2])
+  ism@rownames <- c(ism@rownames, names)
+  ism
+}
+
+##' Helper function to add cols to an existing ISM
+##'
+##' Adds cols of value \code{val} for each entry in \code{names}.
+##' @param ism An ISM.
+##' @param names A vector of names to be added to the cols
+##' @param val The value to be added.
+##' @return The ISM with any padded cols.
+##' @keywords internal
+addIneligibleControls <- function(ism, names, val=Inf) {
+  numnames <- length(names)
+  if (numnames == 0) return(ism) # Short circuit
+  ismdim <- ism@dimension
+  ism@.Data <- c(ism@.Data, rep(val, numnames*ismdim[1]))
+  ism@cols <- c(ism@cols,
+                rep(seq(ismdim[2] + 1, ismdim[2] + numnames),
+                    each=ismdim[1]))
+  ism@rows <- c(ism@rows, rep(1:ismdim[1], times=numnames))
+  ism@dimension <- c(ismdim[1], ismdim[2] + as.integer(numnames))
+  ism@colnames <- c(ism@colnames, names)
+  ism
+}
+
+scoreCaliperBlock <- function(x, z, caliper) {
   z <- toZ(z)
 
   treated <- x[z]
@@ -770,7 +840,7 @@ scoreCaliper <- function(x, z, caliper) {
 
   stops <- findInterval(treated + caliper + .Machine$double.eps, control)
   starts <- length(control) - findInterval(-(treated - caliper -
-                                             .Machine$double.eps), rev(-control))
+                                               .Machine$double.eps), rev(-control))
 
   edges <- pmax(0, (stops - starts))
   n <- sum(edges)

tests/testthat/test.match_on.R

@@ -1067,8 +1067,45 @@ test_that("Exclude argument for match_on with caliper arg", {
     expect_equivalent(mm[, 'z'], 7:9)
 })
 
-
 test_that("No longer support user-defined distances in match_on.formula", {
   data(nuclearplants)
   expect_warning(match_on(pr ~ cost, data = nuclearplants, method = optmatch:::compute_euclidean), "not supported")
 })
+
+test_that("add rows to an ISM", {
+  w <- matrix(c(1, Inf, 2, 3, Inf, 4), nrow = 3)
+  rownames(w) <- c("a", "b", "c")
+  B <- as.InfinitySparseMatrix(w[1:2,1:2])
+  BPlus <- addEligibleTreatments(B, c("c"), c(2, 4))
+  expect_equivalent(as.matrix(BPlus), w)
+})
+
+test_that("add cols to an ISM", {
+  w <- matrix(c(1, Inf, 2, 3, Inf, 4, 5, 6, 7), nrow = 3)
+  colnames(w) <- c("a", "b", "c")
+  B <- as.InfinitySparseMatrix(w[,1:2])
+  BPlus <- addIneligibleControls(B, c("c"), c(5, 6, 7))
+  expect_equivalent(as.matrix(BPlus), w)
+})
+
+test_that("scoreCaliper returns BISM when given a BISM within arg", {
+  scores <- rep(1:3, each = 4)
+  z <- rep(c(0,1), 6)
+  names(z) <- names(scores) <- letters[1:12]
+  b <- rep(1:3, 4)
+
+  ez <- exactMatch(z ~ b)
+  a <- scoreCaliper(scores, z, caliper=1, within=ez)
+  expect_true(is(a, "BlockedInfinitySparseMatrix"))
+  expect_equal(ez@groups[sort(names(ez@groups))],
+               a@groups[sort(names(a@groups))])
+  ## currently getting:
+  ## > ez@groups
+  ## a b c d e f g h i j k l
+  ## 1 2 3 1 2 3 1 2 3 1 2 3
+  ## Levels: 1 2 3
+  ## > a@groups
+  ## a g d j e k b h c i f l
+  ## 0 0 0 0 1 1 1 1 2 2 2 2
+  ## Levels: 0 1 2
+})