sensitivityIxJ: Exact Sensitivity Analysis

• sensitivityIxJ is an exact, nonparametric sensitivity analysis framework for multi-level contingency tables that quantifies how unmeasured confounding alters association tests.
• It extends Rosenbaum’s model to multi-level treatments by deriving the worst-case null distribution for any permutation-invariant test with computational shortcuts for specialized tests.
• Implemented in R, the framework enables robust inference without collapsing treatment or outcome levels, preserving power in observational studies.

sensitivityIxJ is an exact, nonparametric sensitivity-analysis framework for observational studies with $I\times J$ or $I\times J\times K$ contingency tables, designed to assess how inference on association changes under unmeasured confounding when both treatment and outcome may be non-binary. It extends Rosenbaum’s sensitivity model for generic bias, derives the exact worst-case null distribution for any permutation-invariant test, and provides computationally specialized procedures for sign-score, ordinal, and general permutation-invariant statistics. The framework is implemented in the R package sensitivityIxJ and is motivated by the observation that existing sensitivity analyses for contingency tables often assume a binary treatment variable or impose strong parametric assumptions on non-binary treatment variables (Chiu et al., 23 Jul 2025).

1. Statistical setting and inferential target

The framework treats observational data as an $I\times J$ contingency table, or as an $I\times J\times K$ array when stratification is present. For subject $s$, the treatment is denoted $Z_s\in\{1,\ldots,I\}$ and the observed outcome is $r_s\in\{1,\ldots,J\}$. Under no unmeasured confounding and Fisher’s sharp null, treatment assignment is uniform over all assignments with fixed row totals $N_{i\cdot}$:

$$\mathbb{P}(Z=z\mid\mathcal{F},\,\mathcal{Z}(\mathbf{N}_{I\cdot})) =\frac1{|\mathcal{Z}(\mathbf{N}_{I\cdot})|},$$

where $\mathcal{Z}(\mathbf{N}_{I\cdot})$ is the set of all assignments with fixed row totals $I\times J\times K$0.

Under Fisher’s sharp null, the column totals $I\times J\times K$1 are also fixed. This allows the inferential problem to be phrased as a conditional randomization problem over the set of treatment assignments compatible with the observed margins. Any permutation-invariant test statistic can then be written as $I\times J\times K$2, where $I\times J\times K$3 denotes the contingency table.

The central object is the worst-case null $I\times J\times K$4-value at a specified hidden-bias level $I\times J\times K$5. Rather than testing association only under ignorability, sensitivityIxJ asks how large the $I\times J\times K$6-value can become once treatment assignment is allowed to depend on an unobserved confounder subject to Rosenbaum-type bounds. This reformulation is what makes the framework a sensitivity analysis for contingency-table association tests rather than a standard exact test.

2. Extended Rosenbaum model for multi-level treatments

To allow hidden bias, sensitivityIxJ extends Rosenbaum’s “generic-bias” model from binary treatment settings to $I\times J\times K$7 treatment levels. An unobserved confounder $I\times J\times K$8 affects treatment assignment through binary indicators $I\times J\times K$9 that specify which treatments are pooled by the confounder, and a nonnegative bias parameter $I\times J$0 that controls the strength of the distortion:

$I\times J$1

Equivalently, for any $I\times J$2,

$I\times J$3

When two subjects share the same observed covariates, this yields the odds-ratio bound

$I\times J$4

whenever $I\times J$5, and the ratio equals $I\times J$6 if $I\times J$7.

This parameterization is flexible because $I\times J$8 determines which treatment levels are regarded as similarly affected by the hidden confounder. In the example emphasized in the package documentation and application, $I\times J$9 encodes bias affecting any care versus no care. A plausible implication is that the model is especially natural for observational studies in which treatment levels can be grouped into substantively meaningful exposure classes without collapsing the outcome table itself.

3. Exact worst-case null distribution

For a permutation-invariant statistic $I\times J\times K$0, the worst-case null $I\times J\times K$1-value at bias level $I\times J\times K$2 is

$I\times J\times K$3

A key result is that the maximizer $I\times J\times K$4 lies at a corner of the unit cube, $I\times J\times K$5; this is stated as Lemma 2. The resulting combinatorial reduction is central to the exact method (Chiu et al., 23 Jul 2025).

For general permutation-invariant tests, the search for the worst-case $I\times J\times K$6 can be restricted to

$I\times J\times K$7

whose size is at most $I\times J\times K$8 by Theorem 1. For ordinal tests with additive scores $I\times J\times K$9, Theorem 2 further reduces the candidate set to

$s$0

which has size $s$1. For sign-score tests with $s$2, Theorem 3 shows that the search collapses to a single corner,

$s$3

Given a corner $s$4, the exact null $s$5-value is

$s$6

Theorem 4 rewrites this quantity in terms of kernels that count tables sharing specified margins and bias counts:

$s$7

The closed-form recursion formulas in Theorem 4 and Lemmas 1–2 accelerate enumeration by several orders of magnitude.

4. Specialized algorithms and computational structure

The framework provides specialized procedures for three sub-families of tests. For sign-score tests with binary outcome ($s$8), the null distribution reduces to a single extended hypergeometric distribution for $s$9 with weights $Z_s\in\{1,\ldots,I\}$0, and fast multivariate hypergeometric routines compute $Z_s\in\{1,\ldots,I\}$1-values in $Z_s\in\{1,\ldots,I\}$2. For ordinal tests, only $Z_s\in\{1,\ldots,I\}$3 corners need be examined, and dynamic programming over $Z_s\in\{1,\ldots,I\}$4 has size $Z_s\in\{1,\ldots,I\}$5 rather than $Z_s\in\{1,\ldots,I\}$6. For general permutation-invariant tests, the method still searches up to $Z_s\in\{1,\ldots,I\}$7 candidate corners, but kernel recursion remains available for exact null computation.

The package implements both exact enumeration and a fast table-sampling method, SIS-G, which uses the same kernels for importance-sampling weights. This division between exact and sampling-based computation is important because the exact method is designed for size control, whereas the sampling approach addresses computational burden in larger problems.

A useful way to interpret the computational design is as a hierarchy. The most general formulation accommodates any permutation-invariant test; structural constraints such as ordinality or binary outcomes are then exploited to shrink the search space for hidden-bias configurations. This suggests that sensitivityIxJ is not a single test but a computational framework that adapts to the algebraic form of the statistic under study.

5. Power properties, approximations, and package interface

Simulation studies reported for the method show that tests using all $Z_s\in\{1,\ldots,I\}$8 levels, described as “Opt” ordinal tests, uniformly dominate tests that collapse either rows or columns. In all four data-generating scenarios, the $Z_s\in\{1,\ldots,I\}$9 (Opt) test has highest power against bias (Chiu et al., 23 Jul 2025). Power is preserved under collapsing only if two conditions both hold: the test scores $r_s\in\{1,\ldots,J\}$0 are equal on the combined levels, and the bias indicators $r_s\in\{1,\ldots,J\}$1 are identical on those levels. Otherwise power is lost. Cross-cut tests based only on extreme levels can perform well when the true association is strongly “U-shaped” or “monotone” on the extremes, but they lose power in moderate associations.

The paper also distinguishes exact inference from asymptotic approximation. Normal approximations based on exact first and second moments from Proposition 3 and Corollary 2 are much faster, but they can exceed nominal size for small $r_s\in\{1,\ldots,J\}$2 or imbalanced margins. Exact methods are therefore recommended for size control.

The R package exposes these ideas through a compact interface.

Function	Role	Key arguments or output
sensixj_table	Worst-case $r_s\in\{1,\ldots,J\}$3-value at a specified $r_s\in\{1,\ldots,J\}$4	N, delta, gamma, stat, method
sensixj_senstest	Sensitivity curve over a grid of $r_s\in\{1,\ldots,J\}$5 values	returns data frame of $r_s\in\{1,\ldots,J\}$6 vs. $r_s\in\{1,\ldots,J\}$7
sensixj_plot	Visualization of sensitivity curve	plots $r_s\in\{1,\ldots,J\}$8 vs. $r_s\in\{1,\ldots,J\}$9
sensixj_summary	Reporting robustness threshold	prints largest $N_{i\cdot}$0 at which $N_{i\cdot}$1

The documented defaults are method="exact", stat="chi2", and $N_{i\cdot}$2. The package also supports ordinal analyses through a specification such as stat=list(type="ordinal",alpha,beta). In practical terms, this software layer turns the theoretical worst-case randomization calculation into a routine sensitivity workflow for observed contingency tables.

6. Early Childhood Longitudinal Study application

The principal empirical application re-analyzes the effect of pre-kindergarten care on math achievement using data from the Early Childhood Longitudinal Study, Kindergarten Class of 1998–1999. The analysis focuses on Black girls and Black boys from low-income, single-parent families. Treatment has $N_{i\cdot}$3 levels: no non-parental care, relative care, and center-based care. Outcome has $N_{i\cdot}$4 levels: number and shape, relative size, and ordinality and sequence. Gender defines $N_{i\cdot}$5 strata, so the observed structure is two $N_{i\cdot}$6 tables (Chiu et al., 23 Jul 2025).

The chosen test is ordinal, evaluated with two score sets. The “Prior” scores are

$N_{i\cdot}$7

and the “Profile” scores are estimated by alternating MLE on $N_{i\cdot}$8. The sensitivity model uses $N_{i\cdot}$9, corresponding to hidden bias that distinguishes no care from any care.

The reported results show that under no hidden bias ($$\mathbb{P}(Z=z\mid\mathcal{F},\,\mathcal{Z}(\mathbf{N}_{I\cdot})) =\frac1{|\mathcal{Z}(\mathbf{N}_{I\cdot})|},$$0), the $$\mathbb{P}(Z=z\mid\mathcal{F},\,\mathcal{Z}(\mathbf{N}_{I\cdot})) =\frac1{|\mathcal{Z}(\mathbf{N}_{I\cdot})|},$$1-values are approximately $$\mathbb{P}(Z=z\mid\mathcal{F},\,\mathcal{Z}(\mathbf{N}_{I\cdot})) =\frac1{|\mathcal{Z}(\mathbf{N}_{I\cdot})|},$$2 for girls and $$\mathbb{P}(Z=z\mid\mathcal{F},\,\mathcal{Z}(\mathbf{N}_{I\cdot})) =\frac1{|\mathcal{Z}(\mathbf{N}_{I\cdot})|},$$3 for boys. Girls remain significant up to $$\mathbb{P}(Z=z\mid\mathcal{F},\,\mathcal{Z}(\mathbf{N}_{I\cdot})) =\frac1{|\mathcal{Z}(\mathbf{N}_{I\cdot})|},$$4, whereas boys lose significance by $$\mathbb{P}(Z=z\mid\mathcal{F},\,\mathcal{Z}(\mathbf{N}_{I\cdot})) =\frac1{|\mathcal{Z}(\mathbf{N}_{I\cdot})|},$$5. A joint test via the truncated-product method also rejects at $$\mathbb{P}(Z=z\mid\mathcal{F},\,\mathcal{Z}(\mathbf{N}_{I\cdot})) =\frac1{|\mathcal{Z}(\mathbf{N}_{I\cdot})|},$$6. The substantive interpretation given in the study is that girls’ pre-K benefit is robust to moderate hidden bias, whereas boys’ benefit is less robust.

This application illustrates the intended use of sensitivityIxJ: not merely testing association in a multi-level contingency table, but quantifying how much unmeasured confounding would be required to overturn that association. The broader methodological significance is that this can be done without dichotomizing treatments or outcomes and without abandoning exact, non-asymptotic control of the null distribution.