Advanced diagnostics

This vignette is for when you already know which function to call (see vignette("RASbound-intro")) and want to understand its result. It shows how to:

A running example

x0 <- matrix(c(2, 4, 9,
               6, 8, 1,
               3, 5, 7), nrow = 3, byrow = TRUE)
row_targets <- c(25, 20, 25)
col_targets <- c(20, 25, 25)
upper <- matrix(c(Inf, Inf, 6,
                  Inf, Inf, Inf,
                  Inf, Inf, Inf), nrow = 3, byrow = TRUE)  # cap cell [1, 3] at 6

fit <- bras_freeze_hybrid(x0, row_targets, col_targets, upper = upper)
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
fit
#> <bras_fit> method: freeze_hybrid 
#>   converged: TRUE  | feasible: TRUE  | outer passes: 2 
#>   cells fixed at a bound: 1

Reading `fixed_cells`

fit$fixed_cells
#>     row col bound value outer_pass
#> row   1   3 upper     6          1

Each row is a cell that was pinned to a bound during the fit. The column outer_pass says when. Cells fixed in pass 1 were the clearest violations of the unconstrained solution. Cells fixed in later passes became binding only after earlier fixations changed the problem that was left.

This column is your audit trail. It shows that a cell reached its bound because of the equilibrium of the free problem at that point, not because of an arbitrary early choice.

If a cell binds when you did not expect it to, look at the data first. Usually the seed value, or the target of its row or column, is not consistent with the bounds you set.

Diagnosing infeasibility

Sometimes no matrix satisfies both the targets and the bounds. By default (on_infeasible = "error") the function stops with a message. The same information is also available in the result:

too_tight <- matrix(1, nrow = 3, ncol = 3)
fit_bad <- suppressWarnings(  # x0 is also above `upper` here; that is not the point
  bras_freeze_sim(x0, row_targets, col_targets, upper = too_tight,
                  on_infeasible = "warning")
)
fit_bad$feasible
#> [1] FALSE
fit_bad$diagnostics$infeasible_at
#> $rows
#> [1] 1 2 3
#> 
#> $cols
#> [1] 1 2 3

infeasible_at$rows and infeasible_at$cols list the margins that could not be met. Start there instead of scanning the whole matrix. There are two common causes.

The bounds are too tight for a row or column. The sum of upper over its cells is smaller than its target, or the sum of lower is larger. Check it directly:
```
rowSums(too_tight)  # every row can reach at most 3, but the targets are 25, 20, 25
#> [1] 3 3 3
```
A cell that has to contribute is seeded at zero. The next section explains this case. The fix is different: change x0, not the bounds.

`feasible = TRUE` but `converged = FALSE`

This is not the same as proven infeasibility. It means the outer loop reached max_outer before it found a solution or ruled one out:

fit_slow <- bras_freeze_sim(x0, row_targets, col_targets, upper = upper,
                            max_outer = 1, on_infeasible = "warning")
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
#> Warning: Stopped after 1 outer pass(es) without fully resolving bound
#> violations. The returned matrix may still contain cells outside their bounds;
#> consider increasing `max_outer`.
fit_slow$feasible    # TRUE: not proven infeasible
#> [1] TRUE
fit_slow$converged   # FALSE: it just did not finish
#> [1] FALSE

Here a larger max_outer may solve it. See ?bras_fit for the full distinction.

Cells that must be exactly zero

You do not need lower = upper = 0 for a cell that must be zero. Every variant scales cells multiplicatively, so a cell with x0[i, j] = 0 stays exactly 0:

x0_with_zero <- x0
x0_with_zero[1, 1] <- 0
fit_zero <- bras_freeze_sim(x0_with_zero, row_targets, col_targets)
fit_zero$matrix[1, 1]  # exactly 0
#> [1] 0

Two points to keep in mind.

A seed outside its own [lower, upper] range is moved to the nearest bound, with a warning. For most cells this only saves some work. For a cell seeded at 0 with a positive lower, it is necessary, because no iteration could move the cell away from zero:

bad_lower <- matrix(0, 3, 3); bad_lower[1, 1] <- 5
fit_imputed <- bras_freeze_sim(x0_with_zero, row_targets, col_targets,
                               lower = bad_lower)
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,1)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. For cell(s) below `lower` (including
#> any seeded at exactly 0 with a positive `lower`), this substitution is what
#> makes the cell solvable at all, not just a speed optimization. If this was
#> unintentional -- a mismatch between `x0` and the bounds you meant to set --
#> double-check your inputs; the final result may otherwise reflect a seed value
#> you did not intend.
fit_imputed$matrix[1, 1]  # at least 5, not 0
#> [1] 7.015922

Read the warning as a prompt to check your inputs. It usually means that x0 and the bounds do not agree.

A row or column that depends on zero-seeded cells can be infeasible even with the loosest bounds (lower = 0, upper = Inf), because those cells cannot contribute anything. The feasibility check of the bras_freeze_*() functions knows this: it treats the achievable range of a zero-seeded cell as {0}. Such a case is reported like any other infeasibility.

If a cell is zero in your seed only because of a data problem (rounding, suppression, a missing value), remember that RAS can never make it positive. This holds for every form of RAS, bounded or not. If the cell should be able to become positive, seed it with a small positive value instead of 0.

Comparing variants on the same problem

All functions return the same bras_fit structure, so comparing them is direct. This is useful as a sanity check, and it is what a simulation study would automate:

fits <- list(
  freeze_sim  = bras_freeze_sim(x0, row_targets, col_targets, upper = upper),
  freeze_seq  = bras_freeze_seq(x0, row_targets, col_targets, upper = upper),
  naive       = bras_naive(x0, row_targets, col_targets, upper = upper),
  soft        = bras_soft(x0, row_targets, col_targets, upper = upper),
  cap_iter    = bras_cap_iter(x0, row_targets, col_targets, upper = upper)
)
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.

fits$calibrate_raking <- bras_calibrate_raking(x0, row_targets, col_targets, upper = upper)
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
fits$calibrate_linear <- bras_calibrate_linear(x0, row_targets, col_targets, upper = upper)
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
# bras_calibrate_logit() is left out on purpose. The section
# "When bras_calibrate_logit() fails" below uses this same problem.

data.frame(
  method    = names(fits),
  converged = vapply(fits, function(f) f$converged, logical(1)),
  max_row_err = vapply(fits, function(f) max(f$diagnostics$row_margin_error), numeric(1)),
  max_col_err = vapply(fits, function(f) max(f$diagnostics$col_margin_error), numeric(1)),
  dist_from_seed = vapply(fits, function(f) sum(abs(f$matrix - x0)), numeric(1))
)
#>                            method converged  max_row_err  max_col_err
#> freeze_sim             freeze_sim      TRUE 7.531138e-09 3.552714e-15
#> freeze_seq             freeze_seq      TRUE 7.531138e-09 3.552714e-15
#> naive                       naive      TRUE 6.700361e-09 4.934506e-09
#> soft                         soft      TRUE 6.912888e-09 3.552714e-15
#> cap_iter                 cap_iter      TRUE 7.531138e-09 0.000000e+00
#> calibrate_raking calibrate_raking      TRUE 6.883027e-11 3.698375e-11
#> calibrate_linear calibrate_linear      TRUE 7.105427e-15 2.131628e-14
#>                  dist_from_seed
#> freeze_sim             31.35153
#> freeze_seq             31.35153
#> naive                  31.35153
#> soft                   31.00000
#> cap_iter               31.35153
#> calibrate_raking       31.35153
#> calibrate_linear       32.03954

dist_from_seed is the L1 distance to the seed. It is a quick measure of how much each method changed the data. It is not a substitute for the divergence measures of a proper simulation study.

Do not read this column as a ranking of all methods. The freeze family, bras_naive(), bras_soft(), bras_cap_iter() and bras_calibrate_raking() all change the data by multiplicative scaling, like classic RAS. bras_calibrate_linear() minimizes a least-squares distance. A smaller dist_from_seed for it does not mean it is closer to the right answer. It means that distance function chose a different matrix. Which distance to minimize is a modelling choice. Use dist_from_seed only to compare methods of the same family.

`bras_calibrate_linear()` and `bras_calibrate_raking()` give different matrices

Both wrappers meet the same margins, but they do not return the same matrix. This is a reproducible case, with a fixed seed. It is the case used in the package’s own cross-validation test (tests/testthat/test-crossvalidation-and-properties.R):

set.seed(20260919)
x0_cmp <- matrix(round(runif(16, 1, 20)), nrow = 4)
row_targets_cmp <- round(runif(4, 80, 150))
props <- runif(4); props <- props / sum(props)
col_targets_cmp <- round(sum(row_targets_cmp) * props)
col_targets_cmp[4] <- sum(row_targets_cmp) - sum(col_targets_cmp[1:3])

fit_ras_cmp    <- ras(x0_cmp, row_targets_cmp, col_targets_cmp, tol = 1e-10, max_iter = 500)
fit_raking_cmp <- bras_calibrate_raking(x0_cmp, row_targets_cmp, col_targets_cmp)
fit_linear_cmp <- bras_calibrate_linear(x0_cmp, row_targets_cmp, col_targets_cmp)

# raking is very close to our own engine (same distance function):
max(abs(fit_ras_cmp$matrix - fit_raking_cmp$matrix))
#> [1] 5.457679e-11

# linear is not, although it meets the same margins to the same tolerance:
max(abs(fit_ras_cmp$matrix - fit_linear_cmp$matrix))
#> [1] 12.12134
rbind(row_ras = rowSums(fit_ras_cmp$matrix), row_linear = rowSums(fit_linear_cmp$matrix))
#>            [,1] [,2] [,3] [,4]
#> row_ras      92  121  139  135
#> row_linear   92  121  139  135

bras_calibrate_raking() agrees with ras() up to numerical tolerance. This is expected, because both minimize the same multiplicative distance. bras_calibrate_linear() meets the same margins, but some of its cells differ a lot from the ras() result. This is not numerical noise: its least-squares distance prefers another matrix.

Both matrices are correct solutions of the constraints. They differ in which notion of “closest to the seed” selects among all valid matrices. Which one you need depends on what the adjustment is supposed to represent. RASbound cannot decide that for you.

When `bras_calibrate_logit()` fails

bras_calibrate_logit() was left out of the comparison above on purpose. On this problem it fails:

# Self-contained: the running example again, so this chunk gives the same
# result in whatever order the document is run.
x0 <- matrix(c(2, 4, 9,
               6, 8, 1,
               3, 5, 7), nrow = 3, byrow = TRUE)
row_targets <- c(25, 20, 25)
col_targets <- c(20, 25, 25)
upper <- matrix(c(Inf, Inf, 6,
                  Inf, Inf, Inf,
                  Inf, Inf, Inf), nrow = 3, byrow = TRUE)
bras_calibrate_logit(x0, row_targets, col_targets, upper = upper)
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
#> Error:
#> ! survey::calibrate(calfun = "logit") failed: Valores infinitos o faltantes en 'x'. This calibration function requires finite bounds and can be numerically fragile when most cells are only loosely bounded (see the "Logit calibration cannot represent an unbounded cell" section of ?bras_calibrate_logit). Try increasing `unbounded_mult`, switching to bras_calibrate_linear() instead -- it does not require finite bounds and does not need this substitution -- or supplying a tighter, meaningful `upper` yourself.

The reason is the following. The bound upper[1, 3] = 6 is tight, and the seed x0[1, 3] = 9 is above it. The seed is therefore moved to exactly 6 before the calibration starts. The logit calibration needs finite bounds, so every Inf is replaced by unbounded_mult * max(row_targets, col_targets). By default that is 20 * 25 = 500. One cell sits on a real bound of 6 next to cells with a bound of 500, about 80 times larger. With such different scales, the Newton-Raphson solve of the logit distance can diverge, which gives the error above.

In practice, the automatic replacement of Inf works well when the bounds are on a similar scale, or when only a few cells really matter. It does not replace real finite bounds when you know that some cells need a tight limit. If you meet this failure on your own data, use bras_calibrate_raking() or bras_calibrate_linear(). Both accept Inf and need no replacement. Or use one of the bras_freeze_*() functions.

When `bras_freeze_hybrid()`’s per-pass fallback is not enough

bras_freeze_hybrid() fixes all violating cells at once. If that would make a single pass infeasible, it falls back to sequential fixation for that pass only. This per-pass fallback does not make it as robust as bras_freeze_seq() in every case, as this real input shows:

x0 <- matrix(c(29, 16, 22,
               14,  4, 20,
               27, 20, 23,
               29, 23, 13,
               26, 20, 12), nrow = 5, byrow = TRUE)
row_targets <- c(91, 45, 67, 71, 70)
col_targets <- c(133, 89, 122)
upper <- matrix(c(40, 47, 27,
                  17, 19, 13,
                  44, 16, 17,
                  29, 19, 41,
                  25, 16, 33), nrow = 5, byrow = TRUE)

fit_sim <- suppressWarnings(
  bras_freeze_sim(x0, row_targets, col_targets, upper = upper, on_infeasible = "warning")
)
fit_seq <- suppressWarnings(
  bras_freeze_seq(x0, row_targets, col_targets, upper = upper, on_infeasible = "warning")
)
fit_hybrid <- bras_freeze_hybrid(x0, row_targets, col_targets, upper = upper)
#> Warning: `x0` had 6 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (5,1), (3,2), (4,2), (5,2), (2,3), and 1 more); imputed to the
#> nearest bound before solving, as if that value had been supplied as the seed.
#> If this was unintentional -- a mismatch between `x0` and the bounds you meant
#> to set -- double-check your inputs; the final result may otherwise reflect a
#> seed value you did not intend.
#> Warning: `x0` had 6 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (5,1), (3,2), (4,2), (5,2), (2,3), and 1 more); imputed to the
#> nearest bound before solving, as if that value had been supplied as the seed.
#> If this was unintentional -- a mismatch between `x0` and the bounds you meant
#> to set -- double-check your inputs; the final result may otherwise reflect a
#> seed value you did not intend.

data.frame(
  method    = c("sim", "seq", "hybrid"),
  feasible  = c(fit_sim$feasible, fit_seq$feasible, fit_hybrid$feasible),
  converged = c(fit_sim$converged, fit_seq$converged, fit_hybrid$converged),
  n_fixed   = c(nrow(fit_sim$fixed_cells), nrow(fit_seq$fixed_cells), nrow(fit_hybrid$fixed_cells))
)
#>   method feasible converged n_fixed
#> 1    sim    FALSE     FALSE      10
#> 2    seq     TRUE      TRUE       5
#> 3 hybrid     TRUE      TRUE       5

sim ends infeasible after fixing 10 cells, while seq solves the problem with only 5. hybrid also solves it, with 5 fixed cells. The per-pass fallback alone is not enough here. It looks at one pass at a time, and fixing a few more cells than necessary in several different passes can add up to an infeasible run, even if each pass looked fine on its own.

For that reason bras_freeze_hybrid() has a second, global safety net. If the per-pass approach ends without a solution, it retries the whole problem from the start with sequential fixation, as bras_freeze_seq() does, before giving up. That is how hybrid succeeds here, and fit$method shows that this retry was used:

fit_hybrid$method
#> [1] "freeze_hybrid_retry_sequential"

A stall: no violated bound, but the margins do not match

A stall is a third outcome, different from infeasibility and from running out of max_outer. It is a known limitation of alternating-projection methods such as RAS. After enough cells are fixed, the free cells that remain in a row, and the columns they belong to, can reach a true fixed point that does not satisfy the margins. The margin error stops changing at all across iterations. It is not slow convergence.

x0_stall <- matrix(c(4, 20,  5, 49,
                     10, 41, 44, 35,
                     29,  9, 24, 18,
                     12, 30, 16, 46,
                     46, 14,  6, 10), nrow = 5, byrow = TRUE)
row_targets_stall <- c(95, 120, 63, 67, 134)
col_targets_stall <- c(146, 215, 82, 36)
upper_stall <- matrix(c(33, 31, 13, 39,
                        35, 41, 13, 46,
                        17, 34, 48,  6,
                        23, 53, 51, 35,
                        53, 57, 50, 52), nrow = 5, byrow = TRUE)

fit_stall <- suppressWarnings(
  bras_freeze_sim(x0_stall, row_targets_stall, col_targets_stall,
                  upper = upper_stall, on_infeasible = "warning")
)

fit_stall$feasible             # TRUE: not proven infeasible
#> [1] TRUE
fit_stall$converged            # FALSE: the margins are not matched
#> [1] FALSE
fit_stall$diagnostics$stalled  # TRUE: a fixed point, not slow progress
#> [1] TRUE

The feasibility check done before solving cannot detect this. For each row and column separately, it checks whether the bounds of the remaining free cells can reach the residual target. That condition is necessary but not sufficient. A sufficient condition has to hold jointly for every subset of rows and columns (it is related to network-flow feasibility). That check is harder and is not implemented.

What to do. A larger max_iter or max_outer will not help, and the warning says so. A stall often means that no matrix satisfies the targets and the bounds together. The last section of this vignette shows such a case on real data, with a simple test that proves it. You can loosen the bounds of the rows and columns involved. The largest values of fit_stall$diagnostics$row_margin_error and col_margin_error point to them. Or you can accept some residual mismatch, if that is tolerable for your use.

A real infeasible case: a stall with no solution behind it

A real 55 x 55 input-output-style table (rounded to integers, and shipped in inst/extdata as the infeasible_*.csv files) shows a stall on real data, and how to confirm that the problem has no solution.

read_mat <- function(f) {
  as.matrix(utils::read.csv(system.file("extdata", f, package = "RASbound"),
                            header = FALSE))
}
x0i <- unname(read_mat("infeasible_x0.csv"))
upperi <- unname(read_mat("infeasible_upper.csv"))
rti <- read_mat("infeasible_row_targets.csv")[, 1]
cti <- read_mat("infeasible_col_targets.csv")[, 1]

fit_i <- suppressWarnings(
  bras_freeze_seq(x0i, rti, cti, upper = upperi,
                  tol = 1e-6, max_iter = 1000, max_outer = 500,
                  on_infeasible = "warning")
)
fit_i$converged
#> [1] FALSE
fit_i$diagnostics$stalled
#> [1] TRUE

The simple feasibility check compares the target of each row and column with the total capacity of its cells. It passes here:

cap <- ifelse(x0i > 0, upperi, 0)   # effective capacity under RAS
any(rti > rowSums(cap))
#> [1] FALSE
any(cti > colSums(cap))
#> [1] FALSE

That check is necessary but not sufficient, because the cells of a row also have to respect the column targets. A set of cells can put at most min(col_target, capacity) into each column. For row 13:

i <- 13
rti[i]                         # what the row must receive
#> [1] 11737980
sum(cap[i, ])                  # its raw capacity: it looks sufficient
#> [1] 13042200
sum(pmin(cti, cap[i, ]))       # what it can really absorb
#> [1] 11270963

Row 13 must receive more than it can absorb once the column totals are taken into account. So no matrix satisfies this problem, and the stall is the correct diagnosis. (Solving the transportation problem as a linear program or a maximum flow gives the same answer.) If every upper limit is relaxed by 25%, the problem has a solution:

fit_r <- suppressWarnings(
  bras_freeze_seq(x0i, rti, cti, upper = round(upperi * 1.25),
                  tol = 1e-6, max_iter = 1000, max_outer = 500,
                  on_infeasible = "warning")
)
fit_r$converged
#> [1] TRUE