A worked example: a 57 x 57 input-output-style table

library(RASbound)

This vignette applies the package to the kind of problem it was written for. A square, non-negative, input-output-style table must be adjusted to new row and column totals, and every cell must stay inside a plausible range. The data come from a real application. They are rounded to integers and ship with the package in inst/extdata.

read_mat <- function(f) {
  as.matrix(utils::read.csv(system.file("extdata", f, package = "RASbound"),
                            header = FALSE))
}
x0 <- unname(read_mat("example_x0.csv"))
lower <- unname(read_mat("example_lower.csv"))
upper <- unname(read_mat("example_upper.csv"))
row_targets <- read_mat("example_row_targets.csv")[, 1]
col_targets <- read_mat("example_col_targets.csv")[, 1]

dim(x0)
#> [1] 57 57
c(sum(row_targets), sum(col_targets))
#> [1] 115090168 115090168
sum(x0 == 0)                  # exact zeros in the seed
#> [1] 458
sum(lower > 0)                # cells with a positive lower bound
#> [1] 2075
tail(row_targets, 3)          # the last two rows (and columns) have target 0
#> [1] 80197     0     0

This input is demanding for three reasons. Many cells have a positive lower bound as well as an upper bound. The seed contains exact zeros. And the last two rows and columns have a target of exactly 0.

Hard-freeze variants

args <- list(x0 = x0, row_targets = row_targets, col_targets = col_targets,
             lower = lower, upper = upper,
             tol = 1e-6, max_iter = 1000, max_outer = 500)
fit_sim <- suppressWarnings(do.call(bras_freeze_sim, args))
fit_seq <- suppressWarnings(do.call(bras_freeze_seq, args))
fit_hyb <- suppressWarnings(do.call(bras_freeze_hybrid, args))

sapply(list(sim = fit_sim, seq = fit_seq, hybrid = fit_hyb),
       function(f) c(converged = f$converged,
                     feasible = f$feasible,
                     n_fixed = nrow(f$fixed_cells)))
#>           sim seq hybrid
#> converged   1   1      1
#> feasible    1   1      1
#> n_fixed    35  34     35

Every cell should respect its range, and the margins should be matched:

ok <- function(f) c(
  within_bounds = all(f$matrix >= lower - 1e-6 & f$matrix <= upper + 1e-6),
  max_margin_error = max(f$diagnostics$row_margin_error,
                         f$diagnostics$col_margin_error)
)
sapply(list(sim = fit_sim, seq = fit_seq, hybrid = fit_hyb), ok)
#>                           sim          seq       hybrid
#> within_bounds    1.000000e+00 1.000000e+00 1.000000e+00
#> max_margin_error 7.674098e-07 7.674098e-07 7.674098e-07

The variants fix cells in a different order. When several bounds bind at once, they can return slightly different matrices. This is how far apart they are:

max(abs(fit_sim$matrix - fit_seq$matrix))
#> [1] 2803.843
max(abs(fit_sim$matrix - fit_hyb$matrix))
#> [1] 0

Where the bounds bind

plot() shows the fitted matrix and marks the cells that were fixed at a bound. An upward triangle is a cell at its upper bound, and a downward triangle is a cell at its lower bound:

plot(fit_seq)

The exact cells and values are in fit_seq$fixed_cells. If your matrix has row and column names, plot() shows them as axis labels. See ?plot.bras_fit for how to read the colours.

The survey wrappers on the same problem

The bras_calibrate_*() wrappers remove the rows and columns with a zero target before calling survey, and put them back as zeros afterwards. On the original, unrounded version of this table, all three wrappers still failed inside survey::calibrate(). The cause is not known. This is what they return on the rounded data:

try_wrapper <- function(f) {
  tryCatch({
    fit <- suppressWarnings(f(x0, row_targets, col_targets,
                              lower = lower, upper = upper))
    paste("returned; converged =", fit$converged)
  }, error = function(e) paste("failed:", substr(conditionMessage(e), 1, 70)))
}
c(raking = try_wrapper(bras_calibrate_raking),
  linear = try_wrapper(bras_calibrate_linear))
#>                                                                             raking 
#> "failed: survey::calibrate(calfun = \"raking\") failed: Calibration failed. See ?" 
#>                                                                             linear 
#> "failed: survey::calibrate(calfun = \"linear\") failed: Calibration failed. See ?"

If a wrapper fails on a real table like this one, use a bras_freeze_*() function instead.