R/fullmatch.R
In markmfredrickson/optmatch: Functions for Optimal Matching
Defines functions
fullmatch.matrix
fullmatch.numeric
fullmatch.default
fullmatch
setTryRecovery
Documented in
fullmatch
setTryRecovery
#' (Internal) Sets up option to try recovery in \code{fullmatch}.
#'
#' @return NULL
#' @keywords internal
setTryRecovery <- function() {
options("fullmatch_try_recovery" = TRUE)
}
#' Optimal full matching
#'
#' Given two groups, such as a treatment and a control group, and a method of
#' creating a treatment-by-control discrepancy matrix indicating desirability and
#' permissibility of potential matches (or optionally an already created such
#' discrepancy matrix), create optimal full matches of members of the groups.
#' Optionally, incorporate restrictions on matched sets' ratios of treatment to
#' control units.
#'
#' If passing an already created discrepancy matrix, finite entries indicate
#' permissible matches, with smaller discrepancies indicating more desirable
#' matches. The matrix must have row and column names.
#'
#' If it is desirable to create the discrepancies matrix beforehand (for example,
#' if planning on running several different matching schemes), consider using
#' \code{\link{match_on}} to generate the distances. This generic function has
#' several useful methods for handling propensity score models, computing
#' Mahalanobis distances (and other arbitrary distances), and using user supplied
#' functions. These distances can also be combined with those generated by
#' \code{\link{exactMatch}} and \code{\link{caliper}} to create very nuanced
#' matching specifications.
#'
#' The value of \code{tol} can have a substantial effect on computation time;
#' with smaller values, computation takes longer. Not every tolerance can be
#' met, and how small a tolerance is too small varies with the machine and with
#' the details of the problem. If \code{fullmatch} can't guarantee that the
#' tolerance is as small as the given value of argument \code{tol}, then
#' matching proceeds but a warning is issued.
#'
#' By default, \code{fullmatch} will attempt, if the given constraints are
#' infeasible, to find a feasible problem using the same constraints. This
#' will almost surely involve using a more restrictive \code{omit.fraction} or
#' \code{mean.controls}. (This will never automatically omit treatment units.)
#' Note that this does not guarantee that the returned match has the least
#' possible number of omitted subjects, it only gives a match that is feasible
#' within the given constraints. It may often be possible to loosen the
#' \code{omit.fraction} or \code{mean.controls} constraint and still find a
#' feasible match. The auto recovery is controlled by
#' \code{options("fullmatch_try_recovery")}.
#'
#' In full matching problems permitting many-one matches (\code{min.controls}
#' less than 1), the number of controls contributing to matches can exceed
#' what was requested by setting a value of \code{mean.controls} or
#' \code{omit.fraction}. I.e., in this setting \code{mean.controls} sets
#' the minimum ratio of number of controls to number of treatments placed
#' into matched sets.
#'
#' If the program detects that (what it thinks is) a large problem,
#' a warning is issued. Unless you have an older computer, there's a good
#' chance that you can handle larger problems (at the cost of increased
#' computation time). To check the large problem threshold, use
#' \code{\link{getMaxProblemSize}}; to re-set it, use
#' \code{\link{setMaxProblemSize}}.
#'
#' @param x Any valid input to \code{match_on}. \code{fullmatch} will use
#' \code{x} and any optional arguments to generate a distance before performing
#' the matching.
#'
#' If \code{x} is a numeric vector, there must also be passed a vector \code{z}
#' indicating grouping. Both vectors must be named.
#'
#' Alternatively, a precomputed distance may be entered. A matrix of
#' non-negative discrepancies, each indicating the permissibility and
#' desirability of matching the unit corresponding to its row (a 'treatment') to
#' the unit corresponding to its column (a 'control'); or, better, a distance
#' specification as produced by \code{\link{match_on}}.
#'
#' @param min.controls The minimum ratio of controls to treatments that is to
#' be permitted within a matched set: should be non-negative and finite. If
#' \code{min.controls} is not a whole number, the reciprocal of a whole number,
#' or zero, then it is rounded \emph{down} to the nearest whole number or
#' reciprocal of a whole number.
#'
#' When matching within subclasses (such as those created by
#' \code{\link{exactMatch}}), \code{min.controls} may be a named numeric vector
#' separately specifying the minimum permissible ratio of controls to treatments
#' for each subclass. The names of this vector should include names of all
#' subproblems \code{distance}.
#'
#' @param max.controls The maximum ratio of controls to treatments that is
#' to be permitted within a matched set: should be positive and numeric.
#' If \code{max.controls} is not a whole number, the reciprocal of a
#' whole number, or \code{Inf}, then it is rounded \emph{up} to the
#' nearest whole number or reciprocal of a whole number.
#'
#' When matching within subclasses (such as those created by
#' \code{\link{exactMatch}}), \code{max.controls} may be a named numeric vector
#' separately specifying the maximum permissible ratio of controls to treatments
#' in each subclass.
#'
#' @param omit.fraction Optionally, specify what fraction of controls or treated
#' subjects are to be rejected. If \code{omit.fraction} is a positive fraction
#' less than one, then \code{fullmatch} leaves up to that fraction of the control
#' reservoir unmatched. If \code{omit.fraction} is a negative number greater
#' than -1, then \code{fullmatch} leaves up to |\code{omit.fraction}| of the
#' treated group unmatched. Positive values are only accepted if
#' \code{max.controls} >= 1; negative values, only if \code{min.controls} <= 1.
#' If neither \code{omit.fraction} or \code{mean.controls} are specified, then
#' only those treated and control subjects without permissible matches among the
#' control and treated subjects, respectively, are omitted.
#'
#' When matching within subclasses (such as those created by
#' \code{\link{exactMatch}}), \code{omit.fraction} specifies the fraction of
#' controls to be rejected in each subproblem, a parameter that can be made to
#' differ by subclass by setting \code{omit.fraction} equal to a named numeric
#' vector of fractions.
#'
#' At most one of \code{mean.controls} and \code{omit.fraction} can be non-\code{NULL}.
#'
#' @param mean.controls Optionally, specify the average number of controls per
#' treatment to be matched. Must be no less than than \code{min.controls} and no
#' greater than the either \code{max.controls} or the ratio of total number of
#' controls versus total number of treated. Some controls will likely not be
#' matched to ensure meeting this value. If neither \code{omit.fraction} or
#' \code{mean.controls} are specified, then only those treated and control
#' subjects without permissible matches among the control and treated subjects,
#' respectively, are omitted.
#'
#' When matching within subclasses (such as those created by
#' \code{\link{exactMatch}}), \code{mean.controls} specifies the average number of
#' controls per treatment per subproblem, a parameter that can be made to
#' differ by subclass by setting \code{mean.controls} equal to a named numeric
#' vector.
#'
#' At most one of \code{mean.controls} and \code{omit.fraction} can be non-\code{NULL}.
#'
#' @param tol Because of internal rounding, \code{fullmatch} may
#' solve a slightly different matching problem than the one
#' specified, in which the match generated by
#' \code{fullmatch} may not coincide with an optimal solution of
#' the specified problem. \code{tol} times the number of subjects
#' to be matched specifies the extent to
#' which \code{fullmatch}'s output is permitted to differ from an
#' optimal solution to the original problem, as measured by the
#' sum of discrepancies for all treatments and controls placed
#' into the same matched sets.
#'
#' @param data Optional \code{data.frame} or \code{vector} to use to get order
#' of the final matching factor. If a \code{data.frame}, the \code{rownames}
#' are used. If a vector, the \code{names} are first tried, otherwise the contents
#' is considered to be a character vector of names. Useful to pass if you want to
#' combine a match (using, e.g., \code{cbind}) with the data that were used to
#' generate it (for example, in a propensity score matching).
#'
#' @param solver Choose which solver to use. Currently implemented are RELAX-IV
#' and LEMON. Default of \code{""}, a blank string, will use RELAX-IV if the
#' \strong{rrelaxiv} package is installed, otherwise will use LEMON.
#'
#' To explicitly use RELAX-IV, pass string "RELAX-IV".
#'
#' To use LEMON, pass string "LEMON". Optionally, to specify which algorithm
#' LEMON will use, pass the function \link{LEMON} with argument for the
#' algorithm name, "CycleCancelling", "CapacityScaling", "CostScaling", and
#' "NetworkSimplex". See this site for details on their differences:
#' \url{https://lemon.cs.elte.hu/pub/doc/latest/a00606.html}. CycleCancelling is
#' the default.
#'
#' The CycleCancelling algorithm seems to produce results most closely
#' resembling those of optmatch versions prior to 1.0. We have observed the
#' other LEMON algorithms to produce different results when the
#' \code{mean.controls} is unspecified, or specified in such a way as to produce
#' an infeasible matching problem. When using a LEMON algorithm other than
#' CycleCancelling, we recommend setting the \code{fullmatch_try_recovery}
#' option to \code{FALSE}.
#'
#' @param ... Additional arguments, passed to \code{match_on} (e.g. \code{within})
#' or to specific methods.
#'
#' @return A \code{\link{optmatch}} object (\code{factor}) indicating matched groups.
#'
#' @references
#' Hansen, B.B. and Klopfer, S.O. (2006), \sQuote{ Optimal full matching and related designs via network flows},
#' \emph{Journal of Computational and Graphical Statistics}, \bold{15}, 609--627.
#'
#' Hansen, B.B. (2004), \sQuote{Full Matching in an Observational Study
#' of Coaching for the SAT}, \emph{Journal of the American
#' Statistical Association}, \bold{99}, 609--618.
#'
#' Rosenbaum, P. (1991), \sQuote{A Characterization of Optimal Designs for Observational
#' Studies}, \emph{Journal of the Royal Statistical Society, Series B},
#' \bold{53}, 597--610.
#'
#' @example inst/examples/fullmatch.R
#' @keywords nonparametric optimize
#' @export
fullmatch <- function(x,
min.controls = 0,
max.controls = Inf,
omit.fraction = NULL,
mean.controls = NULL,
tol = .001,
data = NULL,
solver = "",
...) {
# if x does not exist then print helpful error msg
x_str <- deparse(substitute(x))
if (length(x_str)>1) x_str <- paste(x_str, collapse="")
data_str <- deparse(substitute(data))
tryCatch(x, error = function(e) {
stop(missing_x_msg(x_str, data_str, ...))})
if (is.null(data)) {
if (is(x, "InfinitySparseMatrix") |
is(x, "matrix") |
is(x, "optmatch.dlist") )
warning("Without 'data' argument the order of the match is not guaranteed
to be the same as your original data.")
}
if (is(x, "optmatch.dlist")) {
warning("The use of 'optmatch.dlist' objects created by 'mdist()' is deprecated.\nPlease use 'match_on()' instead.")
}
UseMethod("fullmatch")
}
#' @export
fullmatch.default <- function(x,
min.controls = 0,
max.controls = Inf,
omit.fraction = NULL,
mean.controls = NULL,
tol = .001,
data = NULL,
solver = "",
within = NULL,
...) {
if (!inherits(x, gsub("match_on.","",methods("match_on")))) {
stop("Invalid input, must be a potential argument to match_on")
}
mfd <- if (!is.null(data)) {
model.frame(data, na.action=na.pass)
} else {
if (inherits(x, "function")) {
stop("A data argument must be given when passing a function")
}
model.frame(x, na.action=na.pass)
}
if (!is(mfd, "data.frame")) {
stop("Please pass data argument")
}
m <- match_on(x, within=within, data=mfd, ...)
out <- fullmatch(m,
min.controls=min.controls,
max.controls=max.controls,
omit.fraction=omit.fraction,
mean.controls=mean.controls,
tol=tol,
data=mfd,
solver=solver,
...)
attr(out, "call") <- match.call()
out
}
#' @export
fullmatch.numeric <- function(x,
min.controls = 0,
max.controls = Inf,
omit.fraction = NULL,
mean.controls = NULL,
tol = .001,
data = NULL,
solver = "",
z,
within = NULL,
...) {
m <- match_on(x, within=within, z=z, ...)
out <- fullmatch(m,
min.controls=min.controls,
max.controls=max.controls,
omit.fraction=omit.fraction,
mean.controls=mean.controls,
tol=tol,
data=data,
solver=solver,
...)
attr(out, "call") <- match.call()
out
}
#' @export
fullmatch.matrix <- function(x,
min.controls = 0,
max.controls = Inf,
omit.fraction = NULL,
mean.controls = NULL,
tol = .001,
data = NULL,
solver = "",
within = NULL,
hint,
...) {
hint <- if (missing(hint)) NULL else nodeinfo(hint)
### Checking Input ###
# this will throw an error if not valid
validDistanceSpecification(x)
# note: we might want to move these checks to validDistSpec
dnms <- dimnames(x)
if (is.null(dnms) | is.null(dnms[[1]]) | is.null(dnms[[2]])) {
stop("argument \'x\' must have dimnames")
}
if (any(duplicated(unlist(dnms)))){
stop("dimnames of argument \'x\' contain duplicates")
}
if (!is.null(within)) warning("Ignoring non-null 'within' argument. When using 'fullmatch' with\n pre-formed distances, please combine them using '+'.")
nmtrt <- dnms[[1]]
nmctl <- dnms[[2]]
# note: this next _should_ be unnecessary, the objects should do this
# but better safe than sorry
if (!isTRUE(all.equal(dim(x), c(length(nmtrt), length(nmctl))))) {
stop("argument \'x\' dimensions do not match row and column names")
}
if (!is.numeric(min.controls)) {
stop("argument \'min.controls\' must be numeric")
}
if (!is.numeric(max.controls)) {
stop("argument \'max.controls\' must be numeric")
}
if (!is.null(omit.fraction)) {
# A vector of all NA's is logical, not numeric, so the first condition is needed.
if (all(is.na(omit.fraction))) {
omit.fraction <- NULL
} else if (any(abs(omit.fraction) > 1, na.rm = TRUE) | !is.numeric(omit.fraction)) {
stop("omit.fraction must be NULL or numeric between -1 and 1")
}
}
if (!is.null(mean.controls)) {
if (all(is.na(mean.controls))) {
mean.controls <- NULL
} else if (any(mean.controls <= 0, na.rm = TRUE) | !is.numeric(mean.controls)) {
stop("mean.controls must be NULL or numeric greater than 0")
}
}
if (!is.null(omit.fraction) & !is.null(mean.controls)) {
stop("omit.fraction and mean.controls cannot both be specified")
}
# Issue #56: Checking for sane input in data
if (!is.null(data)) {
if (!is.vector(data)) {
dnames <- rownames(data)
} else {
dnames <- names(data)
}
if (any(!unlist(dimnames(x)) %in% dnames)) {
stop("Some elements of the distance matrix are not found in the data argument.")
}
}
# problems is guaranteed to be a list of DistanceSpecifictions
# it may only have 1 entry
problems <- findSubproblems(x)
# the number of problems should match the argument lengths for
# min, max, and omit
np <- length(problems)
if (np>1 & is.null(names(problems)))
stop("Subproblems should have names.")
subproblemids <- names(problems)
if (is.null(subproblemids)) subproblemids <- character(1L)
if (is.null(hint)) { hints <- rep(list(NULL), np)
} else {
hints <- split(hint, hint[['groups']],
drop=TRUE # drops levels of hint$groups that aren't represented in hint
)
nohint <- setdiff(subproblemids, names(hints))
hints <- hints[match(subproblemids, names(hints), 0L)]
if (length(hints)>0) for (ii in 1L:length(hints)) hints[[ii]] <- new("NodeInfo", hints[[ii]])
if (length(nohint))
{
nullhint <- rep(list(NULL), length(nohint))
names(nullhint) <- nohint
hints <- c(hints, nullhint)
if (length(nohint)==np) warning("Hint lacks information about subproblems of this problem; ignoring.")
}
hints <- hints[match(subproblemids, names(hints))]
}
if (length(min.controls) > 1 & np != length(min.controls)) {
if (is.null(names(min.controls)))
stop("\'min.controls\' longer than 1 should have names (that match names of subproblems/exact matching categories).")
min.controls <- min.controls[match(subproblemids, names(min.controls), nomatch=0)]
if (length(min.controls)<np)
stop(paste("\'min.controls\' arg of length>1 lacks entries for ",
np-length(min.controls), " of ", np, " subproblems", sep = ""))
}
if (length(max.controls) > 1 & np != length(max.controls)) {
if (is.null(names(max.controls)))
stop("\'max.controls\' longer than 1 should have names (that match names of subproblems/exact matching categories).")
max.controls <- max.controls[match(subproblemids, names(max.controls), nomatch=0)]
if (length(max.controls)<np)
stop(paste("\'max.controls\' arg of length>1 lacks entries for ",
np-length(max.controls), " of ", np, " subproblems", sep = ""))
}
if (!is.null(omit.fraction) & length(omit.fraction) > 1 & np !=
length(omit.fraction)) {
if (is.null(names(omit.fraction)))
stop("\'omit.fraction\' longer than 1 should have names (that match names of subproblems/exact matching categories).")
omit.fraction <- omit.fraction[match(subproblemids, names(omit.fraction), nomatch=0)]
if (length(omit.fraction)<np)
stop(paste("\'omit.fraction\' arg of length>1 lacks entries for ",
np-length(omit.fraction), " of ", np, " subproblems", sep = ""))
}
if (!is.null(mean.controls) & length(mean.controls) > 1 & np !=
length(mean.controls)) {
if (is.null(names(mean.controls)))
stop("\'mean.controls\' longer than 1 should have names (that match names of subproblems/exact matching categories).")
mean.controls <- mean.controls[match(subproblemids, names(mean.controls), nomatch=0)]
if (length(mean.controls)<np)
stop(paste("\'mean.controls\' arg of length>1 lacks entries for ",
np-length(mean.controls), " of ", np, " subproblems", sep = ""))
}
# reset the arguments to be the right length if they are not
if (length(min.controls) == 1) {
min.controls <- rep(min.controls, np)
}
if (length(max.controls) == 1) {
max.controls <- rep(max.controls, np)
}
if (is.null(omit.fraction)) {
omit.fraction <- NA
}
if (length(omit.fraction) == 1) {
omit.fraction <- rep(omit.fraction, np)
}
if (is.null(mean.controls)) {
mean.controls <- NA
}
if (length(mean.controls) == 1) {
mean.controls <- rep(mean.controls, np)
}
if (any(mean.controls + .Machine$double.eps^0.5 < min.controls, na.rm=TRUE)) {
stop("mean.controls cannot be smaller than min.controls")
}
if (any(mean.controls - .Machine$double.eps^0.5> max.controls, na.rm=TRUE)) {
stop("mean.controls cannot be larger than max.controls")
}
if (any(!is.na(mean.controls))) {
if (any(mean.controls > lapply(problems, function(p) {x <- subdim(p)[[1]] ; x[2]/x[1]}), na.rm=TRUE)) {
stop("mean.controls cannot be larger than the ratio of number of controls to treatments")
}
}
if (any(omit.fraction > 0 & max.controls <= .5, na.rm=TRUE)) {
stop("positive \'omit.fraction\' with \'max.controls\' <= 1/2 not permitted")
}
if (any(omit.fraction < 0 & min.controls >= 2, na.rm=TRUE)) {
stop("negative \'omit.fraction\' with \'min.controls\' >= 2 not permitted")
}
# checks solver and evaluates LEMON() if neccessary
solver <- handleSolver(solver)
user.input.mean.controls <- FALSE
if (any(!is.na(mean.controls) & is.na(omit.fraction))) {
user.input.mean.controls <- TRUE
omit.fraction <- 1 - mapply(function(x,y) {z <- subdim(y)[[1]] ; x*z[1]/z[2]}, mean.controls, problems)
}
total.n <- sum(dim(x))
TOL <- tol * total.n
# a helper to handle a single matching problem. all args required.
# input error checking happens in the public fullmatch function.
.fullmatch <- function(d, mnctl, mxctl, omf, hint = NULL, solver) {
# if the subproblem is completely empty, short circuit
if (length(d) == 0 || all(is.infinite(d))) {
x <- dim(d)
cells.a <- rep(NA, x[1])
cells.b <- rep(NA, x[2])
names(cells.a) <- rownames(d)
names(cells.b) <- colnames(d)
tmp <- list(cells = c(cells.a, cells.b), err = -1)
return(tmp)
}
ncol <- dim(d)[2]
nrow <- dim(d)[1]
tol.frac <-
if (total.n > 2 * np) {
(nrow + ncol - 2)/(total.n - 2 * np)
} else 1
# if omf is specified (i.e. not NA), see if is non-negative
# if omf is not specified, check to see if mxctl is > .5
if (switch(1 + is.na(omf), omf >= 0, mxctl > .5)) {
maxc <- min(mxctl, ncol)
minc <- max(mnctl, 1/nrow)
omf.calc <- omf
flipped <- FALSE
} else {
maxc <- min(1/mnctl, ncol)
minc <- max(1/mxctl, 1/nrow)
omf.calc <- -1 * omf
d <- t(d)
flipped <- TRUE
}
## (I'd like to do the following higher in the call stack, and also more
## informatively, obviating subsequent needs for max.cpt, min.cpt,
## omit.fraction etc. Keeping it here for fear of mischief with flipped subproblems.)
if (is.null(hint)) hint <- nodes_shell_fmatch(rownames(d), colnames(d))
temp <- solve_reg_fm_prob(node_info = hint,
distspec = d,
max.cpt = maxc,
min.cpt = minc,
tolerance = TOL * tol.frac,
solver = solver,
omit.fraction = if(!is.na(omf)) { omf.calc }# passes NULL for NA
)
if (!is.null(temp$MCFSolution))
temp$MCFSolution@subproblems[1L,"flipped"] <- flipped
return(temp)
}
# a second helper function, that will attempt graceful recovery in situations where the match
# is infeasible with the given max.controls
.fullmatch.with.recovery <- function(d.r, mnctl.r, mxctl.r, omf.r, hint.r = NULL, solver) {
# if the subproblem isn't clearly infeasible, try to get a match
if (mxctl.r * dim(d.r)[1] >= prod(dim(d.r)[2], 1-omf.r, na.rm=TRUE)) {
tmp <- .fullmatch(d.r, mnctl.r, mxctl.r, omf.r, hint.r, solver)
if (!all(is.na(tmp[1]$cells))) {
# subproblem is feasible with given constraints, no need to recover
new.omit.fraction <<- c(new.omit.fraction, omf.r)
return(tmp)
}
}
# if max.control is in [1, Inf), and we're infeasible
if(is.finite(mxctl.r) & mxctl.r >= 1) {
# Re-solve with no max.control
tmp2 <- list(.fullmatch(d.r, mnctl.r, Inf, omf.r, hint.r, solver))
tmp2.optmatch <- makeOptmatch(d.r, tmp2, match.call(), data)
trial.ss <- stratumStructure(tmp2.optmatch)
treats <- as.numeric(unlist(lapply(strsplit(names(trial.ss), ":"),"[",1)))
ctrls <- as.numeric(unlist(lapply(strsplit(names(trial.ss), ":"),"[",2)))
num.controls <- sum((pmin(ctrls, mxctl.r)*trial.ss)[treats > 0])
if(num.controls == 0) {
# infeasible anyways
if (!exists("tmp")) {
tmp <- .fullmatch(d.r, mnctl.r, mxctl.r, omf.r, hint.r, solver)
}
new.omit.fraction <<- c(new.omit.fraction, omf.r)
return(tmp)
}
new.omf.r <- 1 - num.controls/dim(d.r)[2]
# feasible with the new omit fraction
new.omit.fraction <<- c(new.omit.fraction, new.omf.r)
return(.fullmatch(d.r, mnctl.r, mxctl.r, new.omf.r, hint.r, solver))
} else {
# subproblem is infeasible, but we can't try to fix because no max.controls
if (!exists("tmp")) {
tmp <- .fullmatch(d.r, mnctl.r, mxctl.r, omf.r, hint.r, solver)
}
new.omit.fraction <<- c(new.omit.fraction, omf.r)
return(tmp)
}
}
# In case we need to try and recover from infeasible, save the new.omit.fraction's used for output to user
new.omit.fraction <- numeric(0)
if (is.null(options()$fullmatch_try_recovery)) {
warning("The flag fullmatch_try_recovery is unset, setting to TRUE")
setTryRecovery()
}
if (options()$fullmatch_try_recovery) {
solutions <- mapply(.fullmatch.with.recovery, problems, min.controls, max.controls, omit.fraction, hints, solver, SIMPLIFY = FALSE)
} else {
solutions <- mapply(.fullmatch, problems, min.controls, max.controls, omit.fraction, hints, solver, SIMPLIFY = FALSE)
}
mout <- makeOptmatch(x, solutions, match.call(), data)
names(min.controls) <- names(problems)
names(max.controls) <- names(problems)
attr(mout, "min.controls") <- min.controls
attr(mout, "max.controls") <- max.controls
# length(new.omit.fraction) will be strictly positive if we ever entered .fullmatch.with.recovery
if(length(new.omit.fraction) > 0) {
out.omit.fraction <- new.omit.fraction
} else {
out.omit.fraction <- omit.fraction
}
out.mean.controls <- mapply(function(x,y) (1 - x)*y[2]/y[1], out.omit.fraction, subdim(x))
names(out.mean.controls) <- names(problems)
names(out.omit.fraction) <- names(problems)
if(user.input.mean.controls) {
attr(mout, "mean.controls") <- out.mean.controls
} else {
attr(mout, "omit.fraction") <- out.omit.fraction
}
if(length(new.omit.fraction) > 0 &
!identical(new.omit.fraction, omit.fraction) &
!all(is.na(new.omit.fraction)) &
getOption("optmatch_verbose_messaging", FALSE)) {
if(!any(is.na(new.omit.fraction)) & all(new.omit.fraction == 1)) {
# If we never got a feasible subproblem
warning("The problem appears infeasible with the given constraints.")
} else {
warning("The problem is infeasible with the given constraints; some units were omitted to allow a match.")
}
}
# save hash of distance
attr(mout, "hashed.distance") <- disthash <- hash_dist(x)
attr(mout, "call") <- match.call()
## assemble MCF material
mcfsolutions <- rep(list(NULL), np)
names(mcfsolutions) <- subproblemids
for (ii in 1L:np) {
if (!is.null(solutions[[ii]]$MCFSolution))
{
mcfsolutions[[ii]] <- solutions[[ii]]$MCFSolution
mcfsolutions[[ii]]@subproblems[1L,"hashed_dist"] <- disthash
thesubprob <- subproblemids[ii]
mcfsolutions[[ii]]@subproblems[1L,"groups"] <- thesubprob
mcfsolutions[[ii]]@nodes[,"groups"] <- factor(thesubprob)
if (nrow(mcfsolutions[[ii]]@arcs@matches) > 0) {
mcfsolutions[[ii]]@arcs@matches[,"groups"] <- factor(thesubprob)
}
mcfsolutions[[ii]]@arcs@bookkeeping[,"groups"] <- factor(thesubprob)
bookkeeping_nodes <- c('(_Sink_)', '(_End_)')
for (bn in bookkeeping_nodes) {
nlabs <- node.labels(mcfsolutions[[ii]])
nlabs[nlabs==bn] <- paste0(bn, thesubprob)
node.labels(mcfsolutions[[ii]]) <- nlabs
}
}
}
mcfsolutions <- mcfsolutions[!vapply(mcfsolutions, is.null, logical(1))]
attr(mout, "MCFSolutions") <- if (length(mcfsolutions)==0) {
NULL
} else {
names(mcfsolutions)[1] <- "x"
##b/c in next line `c()` needs to dispatch on an `x` argument
do.call("c", mcfsolutions)
}
# save solver information
attr(mout, "solver") <- solver
if (all(subproblemSuccess(mout) == FALSE)) {
warning(paste("Matching failed.",
"(Restrictions impossible to meet?)\n",
"Enter ?matchfailed for more info."))
} else if (any(subproblemSuccess(mout) == FALSE)) {
warning(paste("At least one subproblem matching failed.\n",
"(Restrictions impossible to meet?)\n",
"Enter ?matchfailed for more info."))
}
return(mout)
}
#' @export
fullmatch.optmatch.dlist <- fullmatch.matrix
#' @export
fullmatch.InfinitySparseMatrix <- fullmatch.matrix
#' @export
fullmatch.BlockedInfinitySparseMatrix <- fullmatch.matrix
#' @aliases fullmatch
#' @rdname fullmatch
#' @export
full <- fullmatch
markmfredrickson/optmatch documentation built on Nov. 24, 2023, 3:38 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions fullmatch.matrix fullmatch.numeric fullmatch.default fullmatch setTryRecovery
Documented in fullmatch setTryRecovery
#' (Internal) Sets up option to try recovery in \code{fullmatch}.
#'
#' @return NULL
#' @keywords internal
setTryRecovery <- function() {
options("fullmatch_try_recovery" = TRUE)
}
#' Optimal full matching
#'
#' Given two groups, such as a treatment and a control group, and a method of
#' creating a treatment-by-control discrepancy matrix indicating desirability and
#' permissibility of potential matches (or optionally an already created such
#' discrepancy matrix), create optimal full matches of members of the groups.
#' Optionally, incorporate restrictions on matched sets' ratios of treatment to
#' control units.
#'
#' If passing an already created discrepancy matrix, finite entries indicate
#' permissible matches, with smaller discrepancies indicating more desirable
#' matches. The matrix must have row and column names.
#'
#' If it is desirable to create the discrepancies matrix beforehand (for example,
#' if planning on running several different matching schemes), consider using
#' \code{\link{match_on}} to generate the distances. This generic function has
#' several useful methods for handling propensity score models, computing
#' Mahalanobis distances (and other arbitrary distances), and using user supplied
#' functions. These distances can also be combined with those generated by
#' \code{\link{exactMatch}} and \code{\link{caliper}} to create very nuanced
#' matching specifications.
#'
#' The value of \code{tol} can have a substantial effect on computation time;
#' with smaller values, computation takes longer. Not every tolerance can be
#' met, and how small a tolerance is too small varies with the machine and with
#' the details of the problem. If \code{fullmatch} can't guarantee that the
#' tolerance is as small as the given value of argument \code{tol}, then
#' matching proceeds but a warning is issued.
#'
#' By default, \code{fullmatch} will attempt, if the given constraints are
#' infeasible, to find a feasible problem using the same constraints. This
#' will almost surely involve using a more restrictive \code{omit.fraction} or
#' \code{mean.controls}. (This will never automatically omit treatment units.)
#' Note that this does not guarantee that the returned match has the least
#' possible number of omitted subjects, it only gives a match that is feasible
#' within the given constraints. It may often be possible to loosen the
#' \code{omit.fraction} or \code{mean.controls} constraint and still find a
#' feasible match. The auto recovery is controlled by
#' \code{options("fullmatch_try_recovery")}.
#'
#' In full matching problems permitting many-one matches (\code{min.controls}
#' less than 1), the number of controls contributing to matches can exceed
#' what was requested by setting a value of \code{mean.controls} or
#' \code{omit.fraction}. I.e., in this setting \code{mean.controls} sets
#' the minimum ratio of number of controls to number of treatments placed
#' into matched sets.
#'
#' If the program detects that (what it thinks is) a large problem,
#' a warning is issued. Unless you have an older computer, there's a good
#' chance that you can handle larger problems (at the cost of increased
#' computation time). To check the large problem threshold, use
#' \code{\link{getMaxProblemSize}}; to re-set it, use
#' \code{\link{setMaxProblemSize}}.
#'
#' @param x Any valid input to \code{match_on}. \code{fullmatch} will use
#' \code{x} and any optional arguments to generate a distance before performing
#' the matching.
#'
#' If \code{x} is a numeric vector, there must also be passed a vector \code{z}
#' indicating grouping. Both vectors must be named.
#'
#' Alternatively, a precomputed distance may be entered. A matrix of
#' non-negative discrepancies, each indicating the permissibility and
#' desirability of matching the unit corresponding to its row (a 'treatment') to
#' the unit corresponding to its column (a 'control'); or, better, a distance
#' specification as produced by \code{\link{match_on}}.
#'
#' @param min.controls The minimum ratio of controls to treatments that is to
#' be permitted within a matched set: should be non-negative and finite. If
#' \code{min.controls} is not a whole number, the reciprocal of a whole number,
#' or zero, then it is rounded \emph{down} to the nearest whole number or
#' reciprocal of a whole number.
#'
#' When matching within subclasses (such as those created by
#' \code{\link{exactMatch}}), \code{min.controls} may be a named numeric vector
#' separately specifying the minimum permissible ratio of controls to treatments
#' for each subclass. The names of this vector should include names of all
#' subproblems \code{distance}.
#'
#' @param max.controls The maximum ratio of controls to treatments that is
#' to be permitted within a matched set: should be positive and numeric.
#' If \code{max.controls} is not a whole number, the reciprocal of a
#' whole number, or \code{Inf}, then it is rounded \emph{up} to the
#' nearest whole number or reciprocal of a whole number.
#'
#' When matching within subclasses (such as those created by
#' \code{\link{exactMatch}}), \code{max.controls} may be a named numeric vector
#' separately specifying the maximum permissible ratio of controls to treatments
#' in each subclass.
#'
#' @param omit.fraction Optionally, specify what fraction of controls or treated
#' subjects are to be rejected. If \code{omit.fraction} is a positive fraction
#' less than one, then \code{fullmatch} leaves up to that fraction of the control
#' reservoir unmatched. If \code{omit.fraction} is a negative number greater
#' than -1, then \code{fullmatch} leaves up to |\code{omit.fraction}| of the
#' treated group unmatched. Positive values are only accepted if
#' \code{max.controls} >= 1; negative values, only if \code{min.controls} <= 1.
#' If neither \code{omit.fraction} or \code{mean.controls} are specified, then
#' only those treated and control subjects without permissible matches among the
#' control and treated subjects, respectively, are omitted.
#'
#' When matching within subclasses (such as those created by
#' \code{\link{exactMatch}}), \code{omit.fraction} specifies the fraction of
#' controls to be rejected in each subproblem, a parameter that can be made to
#' differ by subclass by setting \code{omit.fraction} equal to a named numeric
#' vector of fractions.
#'
#' At most one of \code{mean.controls} and \code{omit.fraction} can be non-\code{NULL}.
#'
#' @param mean.controls Optionally, specify the average number of controls per
#' treatment to be matched. Must be no less than than \code{min.controls} and no
#' greater than the either \code{max.controls} or the ratio of total number of
#' controls versus total number of treated. Some controls will likely not be
#' matched to ensure meeting this value. If neither \code{omit.fraction} or
#' \code{mean.controls} are specified, then only those treated and control
#' subjects without permissible matches among the control and treated subjects,
#' respectively, are omitted.
#'
#' When matching within subclasses (such as those created by
#' \code{\link{exactMatch}}), \code{mean.controls} specifies the average number of
#' controls per treatment per subproblem, a parameter that can be made to
#' differ by subclass by setting \code{mean.controls} equal to a named numeric
#' vector.
#'
#' At most one of \code{mean.controls} and \code{omit.fraction} can be non-\code{NULL}.
#'
#' @param tol Because of internal rounding, \code{fullmatch} may
#' solve a slightly different matching problem than the one
#' specified, in which the match generated by
#' \code{fullmatch} may not coincide with an optimal solution of
#' the specified problem. \code{tol} times the number of subjects
#' to be matched specifies the extent to
#' which \code{fullmatch}'s output is permitted to differ from an
#' optimal solution to the original problem, as measured by the
#' sum of discrepancies for all treatments and controls placed
#' into the same matched sets.
#'
#' @param data Optional \code{data.frame} or \code{vector} to use to get order
#' of the final matching factor. If a \code{data.frame}, the \code{rownames}
#' are used. If a vector, the \code{names} are first tried, otherwise the contents
#' is considered to be a character vector of names. Useful to pass if you want to
#' combine a match (using, e.g., \code{cbind}) with the data that were used to
#' generate it (for example, in a propensity score matching).
#'
#' @param solver Choose which solver to use. Currently implemented are RELAX-IV
#' and LEMON. Default of \code{""}, a blank string, will use RELAX-IV if the
#' \strong{rrelaxiv} package is installed, otherwise will use LEMON.
#'
#' To explicitly use RELAX-IV, pass string "RELAX-IV".
#'
#' To use LEMON, pass string "LEMON". Optionally, to specify which algorithm
#' LEMON will use, pass the function \link{LEMON} with argument for the
#' algorithm name, "CycleCancelling", "CapacityScaling", "CostScaling", and
#' "NetworkSimplex". See this site for details on their differences:
#' \url{https://lemon.cs.elte.hu/pub/doc/latest/a00606.html}. CycleCancelling is
#' the default.
#'
#' The CycleCancelling algorithm seems to produce results most closely
#' resembling those of optmatch versions prior to 1.0. We have observed the
#' other LEMON algorithms to produce different results when the
#' \code{mean.controls} is unspecified, or specified in such a way as to produce
#' an infeasible matching problem. When using a LEMON algorithm other than
#' CycleCancelling, we recommend setting the \code{fullmatch_try_recovery}
#' option to \code{FALSE}.
#'
#' @param ... Additional arguments, passed to \code{match_on} (e.g. \code{within})
#' or to specific methods.
#'
#' @return A \code{\link{optmatch}} object (\code{factor}) indicating matched groups.
#'
#' @references
#' Hansen, B.B. and Klopfer, S.O. (2006), \sQuote{ Optimal full matching and related designs via network flows},
#' \emph{Journal of Computational and Graphical Statistics}, \bold{15}, 609--627.
#'
#' Hansen, B.B. (2004), \sQuote{Full Matching in an Observational Study
#' of Coaching for the SAT}, \emph{Journal of the American
#' Statistical Association}, \bold{99}, 609--618.
#'
#' Rosenbaum, P. (1991), \sQuote{A Characterization of Optimal Designs for Observational
#' Studies}, \emph{Journal of the Royal Statistical Society, Series B},
#' \bold{53}, 597--610.
#'
#' @example inst/examples/fullmatch.R
#' @keywords nonparametric optimize
#' @export
fullmatch <- function(x,
min.controls = 0,
max.controls = Inf,
omit.fraction = NULL,
mean.controls = NULL,
tol = .001,
data = NULL,
solver = "",
...) {
# if x does not exist then print helpful error msg
x_str <- deparse(substitute(x))
if (length(x_str)>1) x_str <- paste(x_str, collapse="")
data_str <- deparse(substitute(data))
tryCatch(x, error = function(e) {
stop(missing_x_msg(x_str, data_str, ...))})
if (is.null(data)) {
if (is(x, "InfinitySparseMatrix") |
is(x, "matrix") |
is(x, "optmatch.dlist") )
warning("Without 'data' argument the order of the match is not guaranteed
to be the same as your original data.")
}
if (is(x, "optmatch.dlist")) {
warning("The use of 'optmatch.dlist' objects created by 'mdist()' is deprecated.\nPlease use 'match_on()' instead.")
}
UseMethod("fullmatch")
}
#' @export
fullmatch.default <- function(x,
min.controls = 0,
max.controls = Inf,
omit.fraction = NULL,
mean.controls = NULL,
tol = .001,
data = NULL,
solver = "",
within = NULL,
...) {
if (!inherits(x, gsub("match_on.","",methods("match_on")))) {
stop("Invalid input, must be a potential argument to match_on")
}
mfd <- if (!is.null(data)) {
model.frame(data, na.action=na.pass)
} else {
if (inherits(x, "function")) {
stop("A data argument must be given when passing a function")
}
model.frame(x, na.action=na.pass)
}
if (!is(mfd, "data.frame")) {
stop("Please pass data argument")
}
m <- match_on(x, within=within, data=mfd, ...)
out <- fullmatch(m,
min.controls=min.controls,
max.controls=max.controls,
omit.fraction=omit.fraction,
mean.controls=mean.controls,
tol=tol,
data=mfd,
solver=solver,
...)
attr(out, "call") <- match.call()
out
}
#' @export
fullmatch.numeric <- function(x,
min.controls = 0,
max.controls = Inf,
omit.fraction = NULL,
mean.controls = NULL,
tol = .001,
data = NULL,
solver = "",
z,
within = NULL,
...) {
m <- match_on(x, within=within, z=z, ...)
out <- fullmatch(m,
min.controls=min.controls,
max.controls=max.controls,
omit.fraction=omit.fraction,
mean.controls=mean.controls,
tol=tol,
data=data,
solver=solver,
...)
attr(out, "call") <- match.call()
out
}
#' @export
fullmatch.matrix <- function(x,
min.controls = 0,
max.controls = Inf,
omit.fraction = NULL,
mean.controls = NULL,
tol = .001,
data = NULL,
solver = "",
within = NULL,
hint,
...) {
hint <- if (missing(hint)) NULL else nodeinfo(hint)
### Checking Input ###
# this will throw an error if not valid
validDistanceSpecification(x)
# note: we might want to move these checks to validDistSpec
dnms <- dimnames(x)
if (is.null(dnms) | is.null(dnms[[1]]) | is.null(dnms[[2]])) {
stop("argument \'x\' must have dimnames")
}
if (any(duplicated(unlist(dnms)))){
stop("dimnames of argument \'x\' contain duplicates")
}
if (!is.null(within)) warning("Ignoring non-null 'within' argument. When using 'fullmatch' with\n pre-formed distances, please combine them using '+'.")
nmtrt <- dnms[[1]]
nmctl <- dnms[[2]]
# note: this next _should_ be unnecessary, the objects should do this
# but better safe than sorry
if (!isTRUE(all.equal(dim(x), c(length(nmtrt), length(nmctl))))) {
stop("argument \'x\' dimensions do not match row and column names")
}
if (!is.numeric(min.controls)) {
stop("argument \'min.controls\' must be numeric")
}
if (!is.numeric(max.controls)) {
stop("argument \'max.controls\' must be numeric")
}
if (!is.null(omit.fraction)) {
# A vector of all NA's is logical, not numeric, so the first condition is needed.
if (all(is.na(omit.fraction))) {
omit.fraction <- NULL
} else if (any(abs(omit.fraction) > 1, na.rm = TRUE) | !is.numeric(omit.fraction)) {
stop("omit.fraction must be NULL or numeric between -1 and 1")
}
}
if (!is.null(mean.controls)) {
if (all(is.na(mean.controls))) {
mean.controls <- NULL
} else if (any(mean.controls <= 0, na.rm = TRUE) | !is.numeric(mean.controls)) {
stop("mean.controls must be NULL or numeric greater than 0")
}
}
if (!is.null(omit.fraction) & !is.null(mean.controls)) {
stop("omit.fraction and mean.controls cannot both be specified")
}
# Issue #56: Checking for sane input in data
if (!is.null(data)) {
if (!is.vector(data)) {
dnames <- rownames(data)
} else {
dnames <- names(data)
}
if (any(!unlist(dimnames(x)) %in% dnames)) {
stop("Some elements of the distance matrix are not found in the data argument.")
}
}
# problems is guaranteed to be a list of DistanceSpecifictions
# it may only have 1 entry
problems <- findSubproblems(x)
# the number of problems should match the argument lengths for
# min, max, and omit
np <- length(problems)
if (np>1 & is.null(names(problems)))
stop("Subproblems should have names.")
subproblemids <- names(problems)
if (is.null(subproblemids)) subproblemids <- character(1L)
if (is.null(hint)) { hints <- rep(list(NULL), np)
} else {
hints <- split(hint, hint[['groups']],
drop=TRUE # drops levels of hint$groups that aren't represented in hint
)
nohint <- setdiff(subproblemids, names(hints))
hints <- hints[match(subproblemids, names(hints), 0L)]
if (length(hints)>0) for (ii in 1L:length(hints)) hints[[ii]] <- new("NodeInfo", hints[[ii]])
if (length(nohint))
{
nullhint <- rep(list(NULL), length(nohint))
names(nullhint) <- nohint
hints <- c(hints, nullhint)
if (length(nohint)==np) warning("Hint lacks information about subproblems of this problem; ignoring.")
}
hints <- hints[match(subproblemids, names(hints))]
}
if (length(min.controls) > 1 & np != length(min.controls)) {
if (is.null(names(min.controls)))
stop("\'min.controls\' longer than 1 should have names (that match names of subproblems/exact matching categories).")
min.controls <- min.controls[match(subproblemids, names(min.controls), nomatch=0)]
if (length(min.controls)<np)
stop(paste("\'min.controls\' arg of length>1 lacks entries for ",
np-length(min.controls), " of ", np, " subproblems", sep = ""))
}
if (length(max.controls) > 1 & np != length(max.controls)) {
if (is.null(names(max.controls)))
stop("\'max.controls\' longer than 1 should have names (that match names of subproblems/exact matching categories).")
max.controls <- max.controls[match(subproblemids, names(max.controls), nomatch=0)]
if (length(max.controls)<np)
stop(paste("\'max.controls\' arg of length>1 lacks entries for ",
np-length(max.controls), " of ", np, " subproblems", sep = ""))
}
if (!is.null(omit.fraction) & length(omit.fraction) > 1 & np !=
length(omit.fraction)) {
if (is.null(names(omit.fraction)))
stop("\'omit.fraction\' longer than 1 should have names (that match names of subproblems/exact matching categories).")
omit.fraction <- omit.fraction[match(subproblemids, names(omit.fraction), nomatch=0)]
if (length(omit.fraction)<np)
stop(paste("\'omit.fraction\' arg of length>1 lacks entries for ",
np-length(omit.fraction), " of ", np, " subproblems", sep = ""))
}
if (!is.null(mean.controls) & length(mean.controls) > 1 & np !=
length(mean.controls)) {
if (is.null(names(mean.controls)))
stop("\'mean.controls\' longer than 1 should have names (that match names of subproblems/exact matching categories).")
mean.controls <- mean.controls[match(subproblemids, names(mean.controls), nomatch=0)]
if (length(mean.controls)<np)
stop(paste("\'mean.controls\' arg of length>1 lacks entries for ",
np-length(mean.controls), " of ", np, " subproblems", sep = ""))
}
# reset the arguments to be the right length if they are not
if (length(min.controls) == 1) {
min.controls <- rep(min.controls, np)
}
if (length(max.controls) == 1) {
max.controls <- rep(max.controls, np)
}
if (is.null(omit.fraction)) {
omit.fraction <- NA
}
if (length(omit.fraction) == 1) {
omit.fraction <- rep(omit.fraction, np)
}
if (is.null(mean.controls)) {
mean.controls <- NA
}
if (length(mean.controls) == 1) {
mean.controls <- rep(mean.controls, np)
}
if (any(mean.controls + .Machine$double.eps^0.5 < min.controls, na.rm=TRUE)) {
stop("mean.controls cannot be smaller than min.controls")
}
if (any(mean.controls - .Machine$double.eps^0.5> max.controls, na.rm=TRUE)) {
stop("mean.controls cannot be larger than max.controls")
}
if (any(!is.na(mean.controls))) {
if (any(mean.controls > lapply(problems, function(p) {x <- subdim(p)[[1]] ; x[2]/x[1]}), na.rm=TRUE)) {
stop("mean.controls cannot be larger than the ratio of number of controls to treatments")
}
}
if (any(omit.fraction > 0 & max.controls <= .5, na.rm=TRUE)) {
stop("positive \'omit.fraction\' with \'max.controls\' <= 1/2 not permitted")
}
if (any(omit.fraction < 0 & min.controls >= 2, na.rm=TRUE)) {
stop("negative \'omit.fraction\' with \'min.controls\' >= 2 not permitted")
}
# checks solver and evaluates LEMON() if neccessary
solver <- handleSolver(solver)
user.input.mean.controls <- FALSE
if (any(!is.na(mean.controls) & is.na(omit.fraction))) {
user.input.mean.controls <- TRUE
omit.fraction <- 1 - mapply(function(x,y) {z <- subdim(y)[[1]] ; x*z[1]/z[2]}, mean.controls, problems)
}
total.n <- sum(dim(x))
TOL <- tol * total.n
# a helper to handle a single matching problem. all args required.
# input error checking happens in the public fullmatch function.
.fullmatch <- function(d, mnctl, mxctl, omf, hint = NULL, solver) {
# if the subproblem is completely empty, short circuit
if (length(d) == 0 || all(is.infinite(d))) {
x <- dim(d)
cells.a <- rep(NA, x[1])
cells.b <- rep(NA, x[2])
names(cells.a) <- rownames(d)
names(cells.b) <- colnames(d)
tmp <- list(cells = c(cells.a, cells.b), err = -1)
return(tmp)
}
ncol <- dim(d)[2]
nrow <- dim(d)[1]
tol.frac <-
if (total.n > 2 * np) {
(nrow + ncol - 2)/(total.n - 2 * np)
} else 1
# if omf is specified (i.e. not NA), see if is non-negative
# if omf is not specified, check to see if mxctl is > .5
if (switch(1 + is.na(omf), omf >= 0, mxctl > .5)) {
maxc <- min(mxctl, ncol)
minc <- max(mnctl, 1/nrow)
omf.calc <- omf
flipped <- FALSE
} else {
maxc <- min(1/mnctl, ncol)
minc <- max(1/mxctl, 1/nrow)
omf.calc <- -1 * omf
d <- t(d)
flipped <- TRUE
}
## (I'd like to do the following higher in the call stack, and also more
## informatively, obviating subsequent needs for max.cpt, min.cpt,
## omit.fraction etc. Keeping it here for fear of mischief with flipped subproblems.)
if (is.null(hint)) hint <- nodes_shell_fmatch(rownames(d), colnames(d))
temp <- solve_reg_fm_prob(node_info = hint,
distspec = d,
max.cpt = maxc,
min.cpt = minc,
tolerance = TOL * tol.frac,
solver = solver,
omit.fraction = if(!is.na(omf)) { omf.calc }# passes NULL for NA
)
if (!is.null(temp$MCFSolution))
temp$MCFSolution@subproblems[1L,"flipped"] <- flipped
return(temp)
}
# a second helper function, that will attempt graceful recovery in situations where the match
# is infeasible with the given max.controls
.fullmatch.with.recovery <- function(d.r, mnctl.r, mxctl.r, omf.r, hint.r = NULL, solver) {
# if the subproblem isn't clearly infeasible, try to get a match
if (mxctl.r * dim(d.r)[1] >= prod(dim(d.r)[2], 1-omf.r, na.rm=TRUE)) {
tmp <- .fullmatch(d.r, mnctl.r, mxctl.r, omf.r, hint.r, solver)
if (!all(is.na(tmp[1]$cells))) {
# subproblem is feasible with given constraints, no need to recover
new.omit.fraction <<- c(new.omit.fraction, omf.r)
return(tmp)
}
}
# if max.control is in [1, Inf), and we're infeasible
if(is.finite(mxctl.r) & mxctl.r >= 1) {
# Re-solve with no max.control
tmp2 <- list(.fullmatch(d.r, mnctl.r, Inf, omf.r, hint.r, solver))
tmp2.optmatch <- makeOptmatch(d.r, tmp2, match.call(), data)
trial.ss <- stratumStructure(tmp2.optmatch)
treats <- as.numeric(unlist(lapply(strsplit(names(trial.ss), ":"),"[",1)))
ctrls <- as.numeric(unlist(lapply(strsplit(names(trial.ss), ":"),"[",2)))
num.controls <- sum((pmin(ctrls, mxctl.r)*trial.ss)[treats > 0])
if(num.controls == 0) {
# infeasible anyways
if (!exists("tmp")) {
tmp <- .fullmatch(d.r, mnctl.r, mxctl.r, omf.r, hint.r, solver)
}
new.omit.fraction <<- c(new.omit.fraction, omf.r)
return(tmp)
}
new.omf.r <- 1 - num.controls/dim(d.r)[2]
# feasible with the new omit fraction
new.omit.fraction <<- c(new.omit.fraction, new.omf.r)
return(.fullmatch(d.r, mnctl.r, mxctl.r, new.omf.r, hint.r, solver))
} else {
# subproblem is infeasible, but we can't try to fix because no max.controls
if (!exists("tmp")) {
tmp <- .fullmatch(d.r, mnctl.r, mxctl.r, omf.r, hint.r, solver)
}
new.omit.fraction <<- c(new.omit.fraction, omf.r)
return(tmp)
}
}
# In case we need to try and recover from infeasible, save the new.omit.fraction's used for output to user
new.omit.fraction <- numeric(0)
if (is.null(options()$fullmatch_try_recovery)) {
warning("The flag fullmatch_try_recovery is unset, setting to TRUE")
setTryRecovery()
}
if (options()$fullmatch_try_recovery) {
solutions <- mapply(.fullmatch.with.recovery, problems, min.controls, max.controls, omit.fraction, hints, solver, SIMPLIFY = FALSE)
} else {
solutions <- mapply(.fullmatch, problems, min.controls, max.controls, omit.fraction, hints, solver, SIMPLIFY = FALSE)
}
mout <- makeOptmatch(x, solutions, match.call(), data)
names(min.controls) <- names(problems)
names(max.controls) <- names(problems)
attr(mout, "min.controls") <- min.controls
attr(mout, "max.controls") <- max.controls
# length(new.omit.fraction) will be strictly positive if we ever entered .fullmatch.with.recovery
if(length(new.omit.fraction) > 0) {
out.omit.fraction <- new.omit.fraction
} else {
out.omit.fraction <- omit.fraction
}
out.mean.controls <- mapply(function(x,y) (1 - x)*y[2]/y[1], out.omit.fraction, subdim(x))
names(out.mean.controls) <- names(problems)
names(out.omit.fraction) <- names(problems)
if(user.input.mean.controls) {
attr(mout, "mean.controls") <- out.mean.controls
} else {
attr(mout, "omit.fraction") <- out.omit.fraction
}
if(length(new.omit.fraction) > 0 &
!identical(new.omit.fraction, omit.fraction) &
!all(is.na(new.omit.fraction)) &
getOption("optmatch_verbose_messaging", FALSE)) {
if(!any(is.na(new.omit.fraction)) & all(new.omit.fraction == 1)) {
# If we never got a feasible subproblem
warning("The problem appears infeasible with the given constraints.")
} else {
warning("The problem is infeasible with the given constraints; some units were omitted to allow a match.")
}
}
# save hash of distance
attr(mout, "hashed.distance") <- disthash <- hash_dist(x)
attr(mout, "call") <- match.call()
## assemble MCF material
mcfsolutions <- rep(list(NULL), np)
names(mcfsolutions) <- subproblemids
for (ii in 1L:np) {
if (!is.null(solutions[[ii]]$MCFSolution))
{
mcfsolutions[[ii]] <- solutions[[ii]]$MCFSolution
mcfsolutions[[ii]]@subproblems[1L,"hashed_dist"] <- disthash
thesubprob <- subproblemids[ii]
mcfsolutions[[ii]]@subproblems[1L,"groups"] <- thesubprob
mcfsolutions[[ii]]@nodes[,"groups"] <- factor(thesubprob)
if (nrow(mcfsolutions[[ii]]@arcs@matches) > 0) {
mcfsolutions[[ii]]@arcs@matches[,"groups"] <- factor(thesubprob)
}
mcfsolutions[[ii]]@arcs@bookkeeping[,"groups"] <- factor(thesubprob)
bookkeeping_nodes <- c('(_Sink_)', '(_End_)')
for (bn in bookkeeping_nodes) {
nlabs <- node.labels(mcfsolutions[[ii]])
nlabs[nlabs==bn] <- paste0(bn, thesubprob)
node.labels(mcfsolutions[[ii]]) <- nlabs
}
}
}
mcfsolutions <- mcfsolutions[!vapply(mcfsolutions, is.null, logical(1))]
attr(mout, "MCFSolutions") <- if (length(mcfsolutions)==0) {
NULL
} else {
names(mcfsolutions)[1] <- "x"
##b/c in next line `c()` needs to dispatch on an `x` argument
do.call("c", mcfsolutions)
}
# save solver information
attr(mout, "solver") <- solver
if (all(subproblemSuccess(mout) == FALSE)) {
warning(paste("Matching failed.",
"(Restrictions impossible to meet?)\n",
"Enter ?matchfailed for more info."))
} else if (any(subproblemSuccess(mout) == FALSE)) {
warning(paste("At least one subproblem matching failed.\n",
"(Restrictions impossible to meet?)\n",
"Enter ?matchfailed for more info."))
}
return(mout)
}
#' @export
fullmatch.optmatch.dlist <- fullmatch.matrix
#' @export
fullmatch.InfinitySparseMatrix <- fullmatch.matrix
#' @export
fullmatch.BlockedInfinitySparseMatrix <- fullmatch.matrix
#' @aliases fullmatch
#' @rdname fullmatch
#' @export
full <- fullmatch
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.