Implementing Rosenbaum Bounds in Statistical Software

Why Software Inputs and Test Statistics Matter

A sensitivity analysis can produce a precise-looking p-value while answering the wrong question if the matched sets are encoded incorrectly or the test statistic is misunderstood. In a Rosenbaum sensitivity analysis, \(\Gamma\) limits how different the treatment assignment odds may be for two units who are matched on measured covariates. The analysis then asks how large a p-value could be under the most unfavorable assignment probabilities allowed by that limit.

This tutorial uses matched pairs and a one-sided sign statistic to make the software calculations transparent. The same principles of checking set structure, direction, and the meaning of the reported bound apply when software uses other statistics. The calculations here are sensitivity analyses of a sharp null, not estimates of a treatment effect and not probabilities that the null is true.

What the Software Needs to Represent

For pair-based calculations, the minimum useful input is a pair identifier, an indicator of which unit was treated, and the outcome for each unit. Keep one record per unit rather than collapsing the data before you have checked the matching structure. Each pair should have exactly two records and exactly one treated unit. A pair identifier must identify the matched set—not the person, treatment group, or row number.

A matched set with more than two units is not automatically a matched pair. Many sensitivity procedures can handle sets with one treated unit and multiple controls, but their assignment probabilities and test-statistic calculations differ from the pairwise sign calculation below. Do not split a larger set into artificial pairs unless the design and analysis specifically justify that operation.

The sign statistic also requires a clear direction. For a hypothesized increase, count treated-minus-control differences above zero. If the software instead forms control-minus-treated differences, the signs reverse. Tied pairs have \(d_i=0\): they contribute neither a positive nor a negative sign and should be handled according to the selected test. For the sign test, omit them and report the resulting number of non-tied pairs.

Assignment Model and Worst-Case Sign-Test Bound

The sensitivity calculation is based on the sharp null hypothesis that treatment has no effect on any unit’s outcome. Under this null, the outcomes are fixed while the treatment label within each pair is varied. The assignment model assumes one treated unit per pair, independent assignment across pairs conditional on the matched design, and treatment-odds ratios within each pair bounded by \(\Gamma\). Under the sharp null and \(\Gamma=1\), each of the two assignments in a pair has probability one-half.

For \(\Gamma>1\), the assignment can favor the member of a pair whose outcome makes the observed direction look more consistent with a treatment effect. In any pair, the probability of that favorable orientation is at most \[ \frac{\Gamma}{1+\Gamma}. \] For the one-sided sign statistic, the worst-case upper-tail probability is therefore the upper tail of a binomial random variable with \(n\) trials and success probability \(\Gamma/(1+\Gamma)\):

This is a worst-case p-value bound for the specified statistic, alternative, matched-pair structure, and \(\Gamma\). It is not a confidence interval for an effect. To see why the bound has this form, let \(I_i\) indicate whether pair \(i\) has the favorable sign under assignment. Each \(I_i\) is a Bernoulli variable whose success probability is at most \(q=\Gamma/(1+\Gamma)\). With independent assignment across pairs, the sum of these indicators is stochastically no larger than a binomial sum with all success probabilities equal to \(q\). Consequently, its probability of reaching or exceeding the observed \(S\) is at most the binomial upper tail shown above. Setting every probability to its allowed maximum gives the worst case for this upper-tail test.

The derivation depends on the sign statistic and the stated assignment model. A signed-rank statistic, a rank-sum statistic, or a statistic using the magnitudes of differences does not generally have this same binomial formula. Use the formula only when the analysis is actually a one-sided sign test on independent matched pairs.

Implementing and Verifying a Gamma Grid

A software workflow should make the analysis reproducible: define the matched sets and outcome direction, choose the statistic, specify a grid of \(\Gamma\) values, and verify selected results independently. For the sign statistic, the software calculation can be checked with a binomial survival function. In R, for example, the upper-tail probability \(P(X\geq S)\) can be computed as `pbinom(S - 1, size = n, prob = Gamma/(1 + Gamma), lower.tail = FALSE)`. Subtracting one is essential because the function’s upper tail is strictly greater than its first argument.

The following grid uses seven non-tied pairs with five positive differences. Every value is calculated from the same statistic and data; only \(\Gamma\) changes.

Practical Checks and Common Failure Modes

A disciplined implementation separates data checks from sensitivity calculations. First confirm the number of matched sets, the number of treated units per set, and the outcome difference direction. Then compute the statistic independently. Finally, compare at least one software output with a hand calculation or an exact distribution calculation. If the software accepts a general matched-set object, verify that its documented test statistic and assignment model match the intended analysis.

• Wrong unit or set identifier: If records from different matched sets share an identifier, the program may treat them as one larger set. Check that every pair has exactly two units and one treated unit.
• Reversed alternative: Counting control-minus-treated positives answers the opposite directional question from counting treated-minus-control positives. Confirm the alternative and the software’s sign convention.
• Including ties as negative signs: A zero difference is not a negative difference. For the sign test, exclude ties from both \(n\) and \(S\).
• Using the wrong tail: A lower-tail or two-sided calculation is not the one-sided upper-tail bound defined here. Check the tail convention and whether the survival-function argument is \(S-1\).
• Changing the statistic without changing the calculation: The binomial tail applies to the sign count. A software option for signed ranks or another statistic requires its own sensitivity calculation.

A common mistake is to see a p-value rise with \(\Gamma\) and report it as the probability that unmeasured confounding exists or that the treatment effect is absent. Neither interpretation is correct. The bound describes the largest one-sided p-value for a specified sharp-null test under the assignment probabilities allowed by that \(\Gamma\); it does not estimate the prevalence of confounding or assign probability to a hypothesis.

The exact binomial sign test is attractive because its inputs and worst-case tail are easy to audit. Other test statistics can use more information in the outcomes, such as their ranks or magnitudes, but may require different software routines and more involved calculations. E-values are another sensitivity metric, but they answer a different question and should not be substituted for a Rosenbaum p-value bound. The choice among approaches should follow the estimand, design, and inferential question, not whichever software output is easiest to obtain.

Finally, \(\Gamma\) is a sensitivity parameter, not something estimated from the matched outcomes. A value such as \(\Gamma=2\) permits treatment assignment odds to differ by as much as a factor of two within a matched pair due to hidden bias, under the model. The interpretation is conditional on the sensitivity model being a useful description of the remaining assignment bias; it does not establish that any particular hidden confounder has that strength.