expect_identical(). Under an optimized BLAS such as Intel
MKL the two differ in the last digit or two, and CRAN’s MKL check failed
three of them. They now compare to a relative tolerance of 1e-10.
Nothing the package computes changed.breakdown_point(),
breakdown_frontier(), identified_set() and
sensitivity_params(), with coef() and
nobs() methods, return the reported quantities without
reaching into the object. breakdown_point() reports the
magnitude, as print() does; signed = TRUE
keeps the sign an Oster delta carries. breakdown_frontier()
names its first column after the parameter swept rather than
index. The vignettes and replication scripts now use
them.print() on a regsensitivity object is
now the headline: the model, the medium-regression coefficient, the
hypothesis, the breakdown point and a thinned results table.
summary() used to print exactly the same thing. It now adds
the call, the control and comparison sets, every sensitivity parameter,
the summary statistics and the full results table; for an identified-set
sweep that table gains a holds column, TRUE
where the set lies wholly on the hypothesis’ side.
regsen_boot(), regsen_multi() and
regsen_summary() results gain the methods they lacked. Each
now has print(), summary() and
plot(): the bootstrap summary adds the bootstrap standard
error, the bias, both intervals and the replicate quantiles, and its
plot draws the replicates with the interval shaded;
confint() returns the interval and recomputes it at another
level or type from the stored replicates. The
multi-treatment summary ranks treatments from most to least fragile and
counts those whose breakdown point is below 1. plot() on a
regsen_summary() result draws the DMP bounds and the Oster
frontier side by side. summary() on these used to fall
through to summary.default().
regsensitivity(formula, data, subcommand = "bounds", ...),
following the R convention for formula interfaces; the subcommand used
to come first, as in Stata. The old order still runs, with a deprecation
warning.u' Var(W1)^{-1} v for
several u, v. They were computed through an
explicit inverse, from chol2inv() or, when the Cholesky
factorization failed, solve() on a ridge-perturbed matrix.
Each vector is now whitened once by a triangular solve on the Cholesky
factor and the forms are crossprod()s of the results. The
fallback for a matrix singular to working precision is the symmetric
eigendecomposition, which gives the pseudo-inverse rather than an
inverse whose entries grow without bound. Results are unchanged to the
digits the numeric snapshots record.plot() on an analysis = "oster",
delta_type = "eq" sweep now groups the roots by branch
rather than by sort order. Up to three roots come back per
delta, sorted within each delta, so the root
continuing a given branch moves between the
beta1/beta2/beta3 columns
wherever a fold adds two roots or one escapes to infinity. Drawing a
column as a line cut continuous branches into pieces and joined
unrelated ones; a heuristic that broke the line at an asymptote and at
delta = 1 suppressed the false connectors but left visible
gaps in curves that are not broken.
Branch identity does not need a heuristic. delta(beta)
is a ratio of cubics; between two consecutive poles of it the map is
continuous and single valued, so the roots in one pole interval are one
branch, and walking that interval in increasing beta
traverses it end to end. A branch that folds back is not a function of
delta, so the branches are drawn with
geom_path().
A zero of the denominator of delta(beta) that is
also a zero of its numerator is a hole, not an asymptote, and no longer
splits a branch. With a single control covariate it is the only real
root of the denominator, which made that case a single continuous curve
reported as two.
The spurious root – the one solving the cubic at every
delta, since it is that common root – no longer survives at
scattered delta and draws a stray spur along the asymptote.
Its residual is set by how accurately polyroot() locates
it, which is well above the few digits of machine epsilon the filter
allowed.
plot() warns when a delta grid is too
coarse to draw anything, and errors when Oster’s cubic has no real
solution anywhere on it. Passing delta = c(-3, 3) – the
endpoints, as Stata’s delta() takes them, rather than the
grid – silently produced an empty panel.
The region rxbar > rmax(cbar) > rybar no
longer raises an error. regsen_bounds() and
regsen_breakdown() compute it like any other, and the
numbers agree with an independent search over the constraint set of DMP
(2026) Theorem 5 to within 0.01, continuously across
rmax.
What stood in the way was the solver for the constraint Assumption A3
places on z. That constraint is a quadratic inequality
whose leading coefficient turns negative exactly when
rxbar * ||c|| > 1; its solution set is then the
complement of an interval, not an interval, and the solver returned the
hull of the two pieces. Sampling that hull hands the optimizer selection
equations that violate A3, so the region was refused rather than
reported. The feasible set is now represented as the list of intervals
it is, and sampled in proportion to their lengths.
The breakdown point in rxbar is no longer capped at
rmax(cbar). That ceiling binds only when rybar
is unrestricted; with rybar finite the breakdown point can
lie well past it. Where no rxbar overturns the conclusion –
the horizontal arm of the breakdown frontier – the reported value is now
+Inf instead of the ceiling.
regsen_breakdown() gains
direction = "rybar", giving the frontier
rybar_bf(rxbar) of DMP (2026) Theorem 4 over a grid of
rxbar. This is the form the paper’s Figure 1 plots, and the
only one that can describe the frontier’s horizontal arm.
regsen_bounds() and regsen_breakdown()
gain clow, the lower end of DMP Assumption A6,
R(W2 ~ W1 . W0) in [clow, cbar]. Asserting
that the controls are at least somewhat endogenous tightens the
identified set and raises the breakdown point. clow = 0,
the default, is the previous behaviour.regsen_bounds(), regsen_breakdown() and
regsen_multi() gain an ncores argument,
joining regsen_boot(), and all four default to the session
setting of the new regsen_cores(). The grid points of an
identified set, the values of a breakdown frontier and the treatments of
a multi-treatment sweep are independent and are spread across the cores;
on the bundled data a finite-rybar frontier runs about
three times faster on six cores. Results are identical for any number of
cores: nothing but the bootstrap draws random numbers, and it seeds each
replicate itself.regsen_cores("auto") sets the session
to all but two of the machine’s cores; regsen_cores(n) sets
a number. Work inside a bootstrap replicate or a treatment always runs
serially, so cores are not oversubscribed by nesting.options(regsensitivity.backend = "psock"), and the
test suite runs it everywhere, so a change that works under fork but not
on Windows is caught before it ships. Arguments a worker closure touches
are forced before the closure is built: a socket worker is a fresh
session and would otherwise receive an unevaluated promise into a frame
it does not have.mclapply() preschedules jobs onto cores and marks
every job on a core as failed when one errors, so a single bad bootstrap
resample used to take a whole core’s worth of replicates down with
it.regsen_boot() now reports a bias-corrected and
accelerated (BCa) interval alongside the percentile interval, and
returns it by default; type = "perc" restores the old
choice. The acceleration comes from a delete-one jackknife over the
units the bootstrap resamples – rows, or clusters under a cluster
bootstrap. In the simulation design of the package’s paper at
n = 500, coverage improves from 92.0% to 94.0% against a
nominal 95% at the same interval width; the two agree by
n = 1000. z0 and acceleration are
returned so the correction can be inspected, and an endpoint that falls
on an extreme replicate is reported as such.model.frame() call per draw, which is what makes the BCa
jackknife – one computation per row – affordable.rmax, the rxbar at which the identified
set becomes unbounded, was wrong for cbar < 1. It solved
the cbar branch of DMP appendix equation (S18)
unconditionally, when the maximising ||c|| is
rxbar itself wherever A6 allows it. With
cbar = 0.5 on the bundled data, calibrating against the ten
geographic and climate covariates, it reported 1.48 where the set is in
fact unbounded from 1.19, so regsen_bounds() printed a
large finite number in place of an infinite one across that range, and
the default rxbar grid ran past the point where the
analysis says anything.beta = bnd_ub(v)) was computed from the lower end of the
identified set. With rybar finite the set is not symmetric
around beta_med, so the reported value was wrong; the
tested end is now chosen by the direction of the hypothesis.beta = bnd_eq(v)) now breaks
down at whichever end of the identified set can reach the hypothesised
value, rather than always at the upper one.compare or nocompare that is not
a control on the right-hand side of the formula – a typo, or the
treatment itself – is now an error. It used to be dropped silently, so a
misspelt covariate quietly changed which variables calibrated the
analysis.cbar must lie
in [0, 1], rxbar and rybar must
be non-negative, r2long non-negative and
maxovb non-negative or NA, and none may be
missing. Out-of-range values used to reach the closed forms and the
optimizer and come back as NA without explanation.analysis = "oster" with
delta_type = "bound", a maxovb cap is now
applied at delta >= 1 too. The raw set is the whole real
line there, and the cap was skipped, so the reported bounds were
-Inf and +Inf where the constraint says
beta_med - maxovb and beta_med + maxovb.W0 were. Var(W1) is singular
otherwise, and the analysis went through a regularised inverse of it.
The test for collinearity with W0 is now relative to the
covariate’s own variance, so a covariate measured in small units is no
longer mistaken for a constant.regsen_boot() keeps replicates on which the breakdown
point is +Inf – the hypothesis survived every value of the
sensitivity parameter – in both intervals, where they count as
+Inf. They were dropped along with failed replicates, so an
interval whose upper tail is unbounded was reported with a finite upper
endpoint. The count is returned as $infinite and
printed.print() of a regsen_boot() result reports
the interval for the magnitude of a signed (Oster) breakdown point
correctly. Taking abs() of each endpoint reversed a
negative interval and, for one straddling zero, hid that the magnitude
may be as small as zero.regsen_table(label = ...) warns when
caption is not given, since knitr::kable()
emits a float only with a caption and the label had nowhere to attach;
it was dropped silently. The r2long column head uses
\mathrm rather than \text, so it no longer
needs amsmath.scale_colour_regsen() and
scale_fill_regsen() recycle the palette (with a warning)
when a sweep has more than eight groups, instead of erroring inside
plot().beta_med
(finite rybar, or Oster).plot() of a result whose sensitivity parameters are all
single values says so, rather than failing with a subscript error.subset containing NA treats
NA as FALSE, as subset() does.
The NA used to select a row of missing values and leave an
NA in the row bookkeeping that regsen_boot()
uses to line up a cluster column.ngrid argument of regsen_bounds() and
regsen_breakdown() is removed: it was documented but never
used.rybar box is
left at its default. A numeric input cannot carry Inf, so
blank now means unrestricted, and a cleared box greys the panel instead
of raising.New vignette, Reading the output: interpretation, edge cases
and decisions (vignette("interpreting-results")). It
says what each number means, gives a decision path from breakdown point
and rho_k to a verdict, and walks through the cases that
confuse readers – a robust and a fragile conclusion on worlds where the
omitted variable exists, a breakdown point of zero, a set that is
unbounded early because the calibration set is weak, what goes in
compare, when cbar matters, an infinite
breakdown point under finite rybar, the two Oster breakdown
points, and how to read the bootstrap print – each with the sentence to
write in the paper. It ends with a list of common mistakes and a
reporting checklist.
The two Oster breakdown points now carry the names Masten and Poirier
(2026) give them, in ?regsen_breakdown, the
interpretation vignette and the paper: the explain away point
(beta = bnd_eq(0); the psacalc number) and the
sign change point (the default). The documentation also states
that the sign change value is capped at one by their Theorem 2, not by
convention, and that only maxovb lifts the cap. On the
frontier data the two are -23.3 and 0.974, and a reader who took one for
the other reached the opposite verdict.
W2,
reads the true beta_long, r_X,
r_Y and c off it by the projections DMP (2026)
define them through, and requires the package’s set at exactly those
parameters to contain beta_long – with rybar
unrestricted, with rybar finite (the optimizer), and with
two-sided A6. A search over constructed omitted variables on fixed
observed data stays inside the rybar = Inf set and reaches
within 15% of its ends, so the set is neither invalid nor loose.beta_med at
rxbar = 0; the breakdown point is where the tested bound
crosses the hypothesis; rmax is where the set becomes
unbounded; rescaling Y or X scales the bounds
as the model says, and rescaling, shifting or reordering the covariates
leaves bounds and breakdown points unchanged, for both DMP and
Oster.rybar identified
set now polishes its result with a local search. Bounds agree with a run
twenty times as long to within 3e-4, where the previous setting could be
off by 0.02.rybar several times faster despite the more accurate
optimizer.The Oster equality identified-set plot no longer draws a false
vertical connector where the numerically ordered solution branches swap
arms across an asymptote. Branches are now split exactly where a step
would cross a vertical asymptote of delta(beta), computed
from the model rather than guessed from the drawing. With several
r2long values the curves are now drawn (and coloured) per
r2long instead of being interleaved into one line.
plot() gains xline, drawing vertical
reference lines. This closes a gap against the Stata package, whose
vignette uses xline() to mark values such as
rmax or a breakdown point; those figures could not be
reproduced in R before.
New theme_regsen() and
scale_colour_regsen() (with
scale_color_regsen() and scale_fill_regsen()
aliases), exported so a plot can be rebuilt with the same look. Defaults
are aimed at print: a faint y-grid, a thin panel border, legend on top,
and the colourblind-safe Okabe-Ito palette, which also survives
greyscale printing.
Axis and legend titles now use plotmath, so they read as the symbols in the papers rather than as parameter names.
plot() gains base_size; use
9 for a two-column journal figure.
Passing NA to title,
subtitle, xtitle or ytitle drops
the annotation. The auto-generated subtitle is no longer applied by
default, since a journal figure carries its description in the
caption.
calibrate_rho() and
calibrate_partial_r2(), computing the point-identified
reference values of DMP (2026) Section 3.4 (Table 4) and Table 3, so the
sensitivity parameters can be read against the observed covariates they
are calibrated to.results$breakdown is documented as
signed. The sign carries the direction of selection;
passing abs() of it back into regsen_bounds()
lands on a different branch of the Oster cubic. The scalar
$breakdown reports the magnitude and exists on
regsen_bounds() output only – a
regsen_breakdown() result carries the value in
results$breakdown alone.rybar = rxbar and the
Table 4 calibration tick marks; both zero crossings are reported and
agree with the published 0.804 and 0.959. The right panel previously
plotted the breakdown point against cbar, a different
object from the paper’s, which is the breakdown frontier in (rxbar,
rybar) space; it is now that frontier, one curve per
cbar.plot() no longer paints a line along the panel edge
when a stretch of the sweep is unbounded. Only the first point of each
unbounded run is routed off-panel and the rest are blank, so the curve
leaves the plot once instead of running along its boundary. Runs are
detected per series, so one beginning at a group boundary still gets its
own exit.regsen_bounds(analysis = "oster", delta_type = "bound")
reported a bounded identified set where the set is in fact unbounded.
The solver correctly returns -Inf/Inf for
|delta| <= dbar with dbar >= 1, but the
routine accumulating bounds across the grid discarded infinite values
alongside NA and carried the last finite bound forward. On
the bundled data, delta = c(0.5, 1, 2) reported
[1.31, 2.98] at every point instead of
(-Inf, Inf) from dbar = 1 on.
NA still carries over, since it means the solver failed
and carries no information; an infinity is a bound and now propagates.
This is the one direction a sensitivity analysis must not err in, since
it made results look more robust than they are. Previously validated
numbers are unchanged – the snapshot tests pass untouched.
regsen_explore(), a shiny application for moving
the sensitivity parameters by hand and watching the identified set
respond. is a Suggests, so the rest of the package works without it.
Inputs are validated and the model is fitted once before the app opens,
so a bad call reports at the console rather than as a banner in a
browser tab.regsen_multi(), running the analysis once per
candidate treatment and collecting the breakdown points into one table,
with a plot() method ordering treatments by fragility.
Every treatment is analysed in the same specification – the others stay
as controls – so the rows are comparable. A treatment whose analysis
fails yields an NA row with the reason rather than aborting
the sweep.regsen_boot() gains ncores, running
replications in parallel: forking on macOS and Linux, a PSOCK cluster on
Windows. A progress bar is not shown in parallel, since it cannot report
meaningfully from several workers.seed are now identical regardless
of ncores. Each replicate seeds itself from a vector drawn
once up front, rather than relying on parallel RNG substreams, which
would have made the numbers depend on how the work happened to be
divided.regsen_tidy() and regsen_glance(),
broom-style tidiers. They are registered on
generics::tidy()/glance() when is installed,
so results flow into and anything else speaking that vocabulary;
generics stays a Suggests.autoplot() method, forwarding to
plot().stopifnot()s the comparison, so a
regression stops the build instead of quietly shipping a wrong
replication.as.data.frame() method returning the results as a
plain data frame, so they can be passed to any table package
(kableExtra, gt, modelsummary, huxtable, tinytable, xtable).regsen_table() rendering results as LaTeX, HTML or
Markdown, so the same call works in an Rmd/Qmd knitting to PDF, HTML or
Markdown. LaTeX output uses booktabs and puts notes in a
threeparttable, with notes = TRUE generating a
note recording the analysis, hypothesis and sample size. Infinite bounds
render as the correct symbol per format rather than being blanked
out.rybar < Inf and cbar > 0) explored
parameter points where varx_bounds() returned NA endpoints.
The downstream quad_ineq_bounds() did unguarded
if (b1 <= r1 && b2 <= r1) comparisons which
evaluated to NA and triggered “missing value where
TRUE/FALSE needed”. Symptoms: BFG2020 vignette failed to build on
Windows (Status: 1 ERROR on win-builder), but built clean on macOS due
to floating-point determinism differences in nloptr’s DIRECT-L.if (...) numeric comparisons in
dmp.R and guard against NA propagation. Every comparison
now uses isTRUE() / isFALSE() so a single NA
cannot poison the dispatch logic.tests/testthat/test-na-robustness.R).rybar = Inf analytic,
rybar < Inf, cbar = 0 analytic,
rybar < Inf, cbar > 0 via DIRECT global optimization
(nloptr)).eq and bound
modes) and breakdown points.maxovb extension.regsen_bounds(),
regsen_breakdown(), regsen_summary(),
regsensitivity() dispatcher, hypothesis helpers
(bnd_lb(), bnd_ub(),
bnd_eq()).plot.regsensitivity() (ggplot2).regsensitivity: end-to-end tour with the BFG2020
application.dmp2022-replication: reproduces every table and figure
of Diegert, Masten & Poirier (2026) that uses the bundled data.mp2022-stylized: reproduces the stylized Oster examples
from Masten & Poirier (2026).bfg2020 (Bazzi, Fiszbein, Gebresilasse
2020, subset).