diff --git a/NEWS.md b/NEWS.md index 0591962a..ab8d00b9 100644 --- a/NEWS.md +++ b/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 diff --git a/R/InfinitySparseMatrix.R b/R/InfinitySparseMatrix.R index 149b4fff..e6c674ca 100644 --- a/R/InfinitySparseMatrix.R +++ b/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)) diff --git a/R/match_on.R b/R/match_on.R index bf93840c..7064f0f8 100644 --- a/R/match_on.R +++ b/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) diff --git a/man/scoreCaliper.Rd b/man/scoreCaliper.Rd index 23e7c342..3074e9f2 100644 --- a/man/scoreCaliper.Rd +++ b/man/scoreCaliper.Rd @@ -5,7 +5,7 @@ \title{(Internal) Helper function to create an InfinitySparseMatrix from a set of scores, a treatment indicator, and a caliper width.} \usage{ -scoreCaliper(x, z, caliper) +scoreCaliper(x, z, caliper, within = NULL) } \arguments{ \item{x}{The scores, a vector indicating the 1-D location of each @@ -15,6 +15,13 @@ unit.} \item{caliper}{The width of the caliper with respect to the scores \code{x}.} + +\item{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.} } \value{ An \code{InfinitySparseMatrix} object, suitable to be diff --git a/tests/testthat/test.match_on.R b/tests/testthat/test.match_on.R index d7b43280..12c71393 100644 --- a/tests/testthat/test.match_on.R +++ b/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 +})