| Title: | Exact Nonparametric Sensitivity Analysis for I by J Contingency Tables |
|---|---|
| Description: | Implements exact, normally approximated, and sampling-based sensitivity analysis for observational studies with contingency tables. Includes exact (kernel-based), normal approximation, and sequential importance sampling (SIS) methods using 'Rcpp' for computational efficiency. The methods build upon the framework introduced in Rosenbaum (2002) <doi:10.1007/978-1-4757-3692-2> and the generalized design sensitivity framework developed by Chiu (2025) <doi:10.48550/arXiv.2507.17207>. |
| Authors: | Elaine Chiu [aut, cre] |
| Maintainer: | Elaine Chiu <[email protected]> |
| License: | GPL-3 |
| Version: | 0.1.5 |
| Built: | 2026-05-14 07:33:10 UTC |
| Source: | https://github.com/cran/sensitivityIxJ |
This function estimates the total number of tables satisfying the margin constraints. We recommend the readers to use this before determining the number of Monte Carlo simulations for the sampling-based p-value calculation.
estimate_tables(treatment.margins, outcome.margins)estimate_tables(treatment.margins, outcome.margins)
treatment.margins |
A vector specifying the total number of subjects receiving each treatment level. The first entry of this vector represents the number of subjects receiving the first treatment level in the study. |
outcome.margins |
A vector specifying the total number of subjects showing each outcome level. The first entry of this vector represents the number of subjects showing the first outcome level in the study. |
this function returns the estimated total number of tables satisfying the margin constraints.
This function computes exact p-values for sensitivity analysis in I by J contingency tables under the generic bias model. It enumerates all tables in the reference set and calculates the maximum p-value over the sensitivity parameter space (u allocations). The function is designed for general permutation-invariant test statistics.
exact.general.sen.IxJ( u_space = NULL, obs.table, table_space, gamma, delta, row = "treatment", verbose = FALSE )exact.general.sen.IxJ( u_space = NULL, obs.table, table_space, gamma, delta, row = "treatment", verbose = FALSE )
u_space |
A matrix where each row represents a unique allocation of |
obs.table |
A matrix or table object representing the observed contingency table. |
table_space |
A list of matrices or table objects representing the space of contingency tables to consider (typically all tables with test statistic >= observed). |
gamma |
A scalar |
delta |
A binary vector with no more than two unique values, corresponding to treatment levels.
The length must match the number of treatments (rows of |
row |
A string indicating whether rows represent "outcome" or "treatment". Must be either "outcome" or "treatment". Default is "treatment". |
verbose |
A logical flag indicating whether to print progress messages showing the
current maximizer and probability at each step. Default is |
The function performs exact sensitivity analysis by:
Enumerating all possible u allocations (or using provided u_space)
Computing the p-value for each allocation by summing probabilities over table_space
Finding the allocation that maximizes the p-value
For computational efficiency, the function only supports certain table dimensions. If u_space is not provided, default corner allocations are generated for J <= 5.
A list containing:
Probability under Randomized Controlled Trial (RCT) with u_allocation set to zero.
Maximum probability found across all allocations in u_space.
The u_allocation vector that yields max.prob.
Extracted gamma value from the generic bias model.
remind the practitioners of their delta
## Example 1: 3 by 3 table with custom test statistic # Create an observed table (example data) obs.table <- matrix(c(5, 3, 2, 6, 11, 7, 3, 0, 3), ncol = 3, byrow = TRUE) # Define a test statistic emphasizing certain cells transform.fun <- function(tb){ test.stat <- 4 * tb[3, 3] + 3 * tb[2, 3] return(test.stat) } obs.stat <- transform.fun(obs.table) # Find the reference set (tables with test statistic >= observed) table.set <- possible.table( threshold = obs.stat, table = obs.table, direction = "greater than", transform.fun = transform.fun ) # Perform sensitivity analysis sen.result <- exact.general.sen.IxJ( obs.table = obs.table, table_space = table.set, gamma = 0.5, delta = c(0, 1, 1) ) sen.result## Example 1: 3 by 3 table with custom test statistic # Create an observed table (example data) obs.table <- matrix(c(5, 3, 2, 6, 11, 7, 3, 0, 3), ncol = 3, byrow = TRUE) # Define a test statistic emphasizing certain cells transform.fun <- function(tb){ test.stat <- 4 * tb[3, 3] + 3 * tb[2, 3] return(test.stat) } obs.stat <- transform.fun(obs.table) # Find the reference set (tables with test statistic >= observed) table.set <- possible.table( threshold = obs.stat, table = obs.table, direction = "greater than", transform.fun = transform.fun ) # Perform sensitivity analysis sen.result <- exact.general.sen.IxJ( obs.table = obs.table, table_space = table.set, gamma = 0.5, delta = c(0, 1, 1) ) sen.result
This function computes exact p-values for score-based sensitivity analysis in I by J contingency tables under the generic bias model. It is specifically designed for ordinal data where both treatment levels and outcomes have a natural ordering, and tests for trend using assigned scores.
exact.score.sen.IxJ( obs.table, gamma, delta, row = "treatment", verbose = FALSE, treatment.scores, outcome.scores )exact.score.sen.IxJ( obs.table, gamma, delta, row = "treatment", verbose = FALSE, treatment.scores, outcome.scores )
obs.table |
A matrix or table object representing the observed contingency table. |
gamma |
A nonnegative scalar. |
delta |
A binary vector with no more than two unique values, corresponding
to treatment levels. The length must match the number of treatments
(rows of |
row |
A string indicating whether rows represent "outcome" or "treatment". Must be either "outcome" or "treatment". Default is "treatment". |
verbose |
A logical flag indicating whether to print progress messages.
Default is |
treatment.scores |
A numeric vector specifying the scores for each treatment level. Must have the same length as the number of treatments and exhibit a monotone trend. |
outcome.scores |
A numeric vector specifying the scores for each outcome level. Must have the same length as the number of outcomes and exhibit a monotone trend. |
The score test assumes both treatments and outcomes are ordinal with monotone trends. The test statistic is computed as the sum of products of cell counts with their corresponding treatment and outcome scores. The function automatically generates an appropriate u-space for score tests, focusing on allocations that respect the ordinal nature of the data.
A list containing:
Probability under Randomized Controlled Trial (RCT) with u_allocation set to zero.
Maximum probability found across all allocations in u_space.
The u_allocation vector that yields max.prob.
Extracted gamma value from the generic bias model.
The delta vector
The observed test statistic based on the observed table and the assigned scores.
The observed table.
exact.general.sen.IxJ for general test statistics,
possible.table for generating reference sets
## Example 1: Binary outcome table (2 by 2) obs.table <- matrix(c(12, 18, 17, 3), ncol = 2, byrow = TRUE, dimnames = list(treatment = c("control", "treated"), outcome = c("failure", "success"))) # Perform score-based sensitivity analysis result_2x2 <- exact.score.sen.IxJ(obs.table = obs.table, gamma = 0.5, delta = c(0, 1), treatment.scores = c(0, 1), outcome.scores = c(0, 1)) result_2x2 ## Example 2: Three-level ordinal outcome (3 by 3) obs.table <- matrix(c(12, 18, 17, 3, 12, 25, 0, 3, 4), ncol = 3, byrow = FALSE, dimnames = list(treatment = c("low", "medium", "high"), outcome = c("poor", "fair", "good"))) # Test for trend with ordinal scores result_3x3 <- exact.score.sen.IxJ(obs.table = obs.table, gamma = 0.5, delta = c(0, 1, 1), treatment.scores = c(0, 1, 2), outcome.scores = c(1, 2, 3)) result_3x3## Example 1: Binary outcome table (2 by 2) obs.table <- matrix(c(12, 18, 17, 3), ncol = 2, byrow = TRUE, dimnames = list(treatment = c("control", "treated"), outcome = c("failure", "success"))) # Perform score-based sensitivity analysis result_2x2 <- exact.score.sen.IxJ(obs.table = obs.table, gamma = 0.5, delta = c(0, 1), treatment.scores = c(0, 1), outcome.scores = c(0, 1)) result_2x2 ## Example 2: Three-level ordinal outcome (3 by 3) obs.table <- matrix(c(12, 18, 17, 3, 12, 25, 0, 3, 4), ncol = 3, byrow = FALSE, dimnames = list(treatment = c("low", "medium", "high"), outcome = c("poor", "fair", "good"))) # Test for trend with ordinal scores result_3x3 <- exact.score.sen.IxJ(obs.table = obs.table, gamma = 0.5, delta = c(0, 1, 1), treatment.scores = c(0, 1, 2), outcome.scores = c(1, 2, 3)) result_3x3
This function computes the probability of a single contingency table under the generic bias model given an unmeasured confounder. This is an auxiliary function for the p-value computation.
generic.I.by.J.sensitivity.point.probability( table, row = "treatment", u_allocation, gamma, delta, shared_divisor = 1e+06 )generic.I.by.J.sensitivity.point.probability( table, row = "treatment", u_allocation, gamma, delta, shared_divisor = 1e+06 )
table |
A matrix or table object representing the observed contingency table. |
row |
A string indicating whether rows represent "outcome" or "treatment". Must be either "outcome" or "treatment". Default is "treatment". |
u_allocation |
A vector where each entry represents the number of |
gamma |
A scalar |
delta |
A binary vector with no more than two unique values, corresponding to treatment levels.
The length must match the number of treatments (rows of |
shared_divisor |
A scalar to rescale the numerator and the denominator of the probability mass function to prevent overflow. Defaulted to 1000000. |
This function returns the probability mass of this table given a unmeasured confounder.
This function performs a score test for ordinal association in contingency tables
with fixed margins under a sensitivity analysis framework. It uses the test statistic
where (treatment score X outcome score)
and computes the mean and variance under the null hypothesis using the Poisson-binomial
approximation for the multivariate Fisher's noncentral hypergeometric distribution.
norm_single_u_allocation_p_value( obs.table, gamma, delta, u_allocation, row = "treatment", treatment.scores, outcome.scores, shared_divisor = 1e+06 )norm_single_u_allocation_p_value( obs.table, gamma, delta, u_allocation, row = "treatment", treatment.scores, outcome.scores, shared_divisor = 1e+06 )
obs.table |
A matrix or table object representing the contingency table. Rows represent treatments, columns represent outcomes. |
gamma |
A nonnegative scalar |
delta |
A binary vector. Must have length equal to the number of treatments. Can contain at most two unique values. |
u_allocation |
A numeric vector of unmeasured confounder allocations for each outcome category. Must have length equal to the number of outcomes. Each entry must be non-negative and not exceed the corresponding outcome margin. |
row |
Character indicating whether rows represent "treatment" or "outcome". If "outcome", the table will be transposed. Default is "treatment". |
treatment.scores |
A numeric vector of scores for treatments. Must be monotone (either increasing or decreasing). Higher scores typically indicate more intense treatments. |
outcome.scores |
A numeric vector of scores for outcomes. Must be monotone (either increasing or decreasing). Higher scores typically indicate better outcomes. |
shared_divisor |
A numeric value used for numerical stability in the Rcpp computations. Default is 1e6. |
The function implements a score test for ordinal association where the test statistic is a weighted sum of the cell counts, with weights given by the product of treatment and outcome scores. The function uses specialized Rcpp functions to compute the mean and variance-covariance structure of the free cells (those with degrees of freedom), then recovers the full mean vector and covariance matrix using the marginal constraints. The test statistic and its moments are then computed using matrix operations.
A list containing:
T_obs |
The observed test statistic |
mu_T |
The expected value of the test statistic under the null hypothesis |
var_T |
The variance of the test statistic under the null hypothesis |
z_score |
The standardized test statistic (T_obs - mu_T) / sqrt(var_T) |
p_value |
The one sided upper tailed p-value based on the normal approximation |
treatment.scores |
The treatment scores used |
outcome.scores |
The outcome scores used |
# 2x3 contingency table obs.table <- matrix(c(10, 20, 30, 15, 25, 10), nrow = 2, byrow = TRUE) # Sensitivity parameters and u_allocation gamma = 0.5 delta <- c(0, 1) u_allocation <- c(5, 10, 8) # Ordinal scores treatment.scores <- c(0, 1) outcome.scores <- c(0, 1, 2) # Perform test result <- norm_single_u_allocation_p_value( obs.table = obs.table, gamma = gamma, delta = delta, u_allocation = u_allocation, treatment.scores = treatment.scores, outcome.scores = outcome.scores )# 2x3 contingency table obs.table <- matrix(c(10, 20, 30, 15, 25, 10), nrow = 2, byrow = TRUE) # Sensitivity parameters and u_allocation gamma = 0.5 delta <- c(0, 1) u_allocation <- c(5, 10, 8) # Ordinal scores treatment.scores <- c(0, 1) outcome.scores <- c(0, 1, 2) # Perform test result <- norm_single_u_allocation_p_value( obs.table = obs.table, gamma = gamma, delta = delta, u_allocation = u_allocation, treatment.scores = treatment.scores, outcome.scores = outcome.scores )
This function implements normal approximation methods for sensitivity analysis in I by J contingency tables under the generic bias model. It computes asymptotically valid p-values for score test statistics based on the product of treatment and outcome scores, providing rapid analysis for large tables.
norm.score.sen.IxJ( obs.table, gamma, delta, row = "treatment", treatment.scores, outcome.scores, shared_divisor = 1e+06, u_space = NULL, verbose = FALSE )norm.score.sen.IxJ( obs.table, gamma, delta, row = "treatment", treatment.scores, outcome.scores, shared_divisor = 1e+06, u_space = NULL, verbose = FALSE )
obs.table |
A matrix or table object representing the observed contingency table. |
gamma |
a nonnegative scalar. |
delta |
a binary vector
to treatment levels. Its length must match the number of treatments
(rows of |
row |
A string indicating whether rows represent "outcome" or "treatment". Must be either "outcome" or "treatment". Default is "treatment". |
treatment.scores |
A numeric vector of scores for treatments. Must be monotone (either increasing or decreasing). Higher scores typically indicate more intense treatments. Length must equal the number of treatments. |
outcome.scores |
A numeric vector of scores for outcomes. Must be monotone (either increasing or decreasing). Higher scores typically indicate better outcomes. Length must equal the number of outcomes. |
shared_divisor |
Numeric value used for numerical stability in calculations. Default is 1e6. |
u_space |
A numeric matrix where each row is a candidate |
verbose |
Logical; if |
For an I by J table, the test statistic is a weighted
sum of cell counts where weights are products of treatment and outcome scores:
across all cells.
The method computes:
Mean and variance of the test statistic under the generic bias model
Standardized z-scores assuming asymptotic normality
the one-sided upper tailed probability
When u_space is not provided, the function automatically generates
corner u allocations which often
contain the worst-case scenarios for sensitivity analysis.
A list containing:
The observed test statistic value.
Mean of the test statistic under RCT (gamma = 0) or (Gamma = 1).
Mean of the test statistic under the sensitivity model at maximizer.
Variance of the test statistic under RCT.
Variance of the test statistic under the sensitivity model at maximizer.
P-value under RCT (no unmeasured confounding).
Maximum p-value across all u-allocations (sensitivity bound).
The u-allocation vector that yields max.prob.
The treatment scores used in the analysis.
The outcome scores used in the analysis.
exact.general.sen.IxJ for exact methods,
sampling.general.sen.IxJ for Monte Carlo methods
# 2 by 3 table example with ordinal scores obs.table <- matrix(c(10, 20, 30, 15, 25, 10), nrow = 2, byrow = TRUE) treatment.scores <- c(0, 1) # Control vs Treatment outcome.scores <- c(0, 1, 2) # Ordinal outcomes result <- norm.score.sen.IxJ(obs.table = obs.table, gamma = 0.5, delta = c(0, 1), treatment.scores = treatment.scores, outcome.scores = outcome.scores ) # 3 by 3 table with customized scores obs.table <- matrix(data=c(10,30,10,14,15,24,4,5,15), nrow = 3) treatment.scores <- c(0, 0.5, 1) # Three treatment levels outcome.scores <- c(0, 1, 2) # three outcome levels result <- norm.score.sen.IxJ(obs.table = obs.table, gamma = 0.5, delta = c(0, 0, 1), treatment.scores = treatment.scores, outcome.scores = outcome.scores )# 2 by 3 table example with ordinal scores obs.table <- matrix(c(10, 20, 30, 15, 25, 10), nrow = 2, byrow = TRUE) treatment.scores <- c(0, 1) # Control vs Treatment outcome.scores <- c(0, 1, 2) # Ordinal outcomes result <- norm.score.sen.IxJ(obs.table = obs.table, gamma = 0.5, delta = c(0, 1), treatment.scores = treatment.scores, outcome.scores = outcome.scores ) # 3 by 3 table with customized scores obs.table <- matrix(data=c(10,30,10,14,15,24,4,5,15), nrow = 3) treatment.scores <- c(0, 0.5, 1) # Three treatment levels outcome.scores <- c(0, 1, 2) # three outcome levels result <- norm.score.sen.IxJ(obs.table = obs.table, gamma = 0.5, delta = c(0, 0, 1), treatment.scores = treatment.scores, outcome.scores = outcome.scores )
Given an observed contingency table (with fixed margins) and a user-defined
function that computes a test statistic, this function enumerates all
valid contingency tables that meet a particular criterion (i.e., those
with T(table) >= threshold or T(table) <= threshold).
possible.table( threshold, table, direction = c("greater than", "less than"), transform.fun )possible.table( threshold, table, direction = c("greater than", "less than"), transform.fun )
threshold |
A numeric value representing the cutoff for the test statistic. |
table |
A matrix or table object specifying the observed contingency table. The function will preserve the row/column sums (margins) of this table when generating possible tables. |
direction |
A character string, either |
transform.fun |
A user-defined function of the form |
The function systematically iterates over all valid ways to fill a contingency
table of dimension I x J such that row and column sums match the original
table's margins. Only those tables that satisfy transform.fun(tbl) >= threshold
(for direction = "greater than") or transform.fun(tbl) <= threshold
(for direction = "less than") are returned.
Currently, this function supports certain table dimensions (I up to 5
when J=2, etc.). If the dimension is not supported, it returns
NULL with a warning.
A list of all contingency tables (each table is returned as a
table object) that meet the specified threshold criterion.
If no tables match or the dimension is not supported, NULL is returned.
# Suppose we have a 3x3 table: obs_table <- matrix(c(5, 3, 2, 6, 11, 7, 3, 0,3), ncol = 3, byrow = TRUE) obs_table # Define a simple test statistic function (e.g., sum of diagonal) diag_sum <- function(tbl) sum(diag(tbl)) # Find all 3x3 tables with the same margins where the diagonal sum >= 19 result <- possible.table(threshold = 19, table = obs_table, direction = "greater than", transform.fun = diag_sum) result[[1]] # Inspect the first matching table# Suppose we have a 3x3 table: obs_table <- matrix(c(5, 3, 2, 6, 11, 7, 3, 0,3), ncol = 3, byrow = TRUE) obs_table # Define a simple test statistic function (e.g., sum of diagonal) diag_sum <- function(tbl) sum(diag(tbl)) # Find all 3x3 tables with the same margins where the diagonal sum >= 19 result <- possible.table(threshold = 19, table = obs_table, direction = "greater than", transform.fun = diag_sum) result[[1]] # Inspect the first matching table
This function implements a sampling-based approach to sensitivity analysis for general (permutation-invariant) test statistics in I by J contingency tables under the generic bias model. It uses Sequential Importance Sampling (SIS) from Eisinger and Chen (2017) to approximate p-values when exact enumeration is computationally infeasible.
sampling.general.sen.IxJ( obs.table, gamma, delta, row = "treatment", mc.iteration = 5000, transform.fun, verbose = FALSE, u_space = NULL )sampling.general.sen.IxJ( obs.table, gamma, delta, row = "treatment", mc.iteration = 5000, transform.fun, verbose = FALSE, u_space = NULL )
obs.table |
A matrix or table object representing the observed contingency table. |
gamma |
a nonnegative scalar |
delta |
A binary vector with no more than two unique values,
corresponding to treatment levels. The length must match the number of
treatments (rows if |
row |
A string indicating whether rows represent "outcome" or "treatment". Must be either "outcome" or "treatment". Default is "treatment". |
mc.iteration |
Integer specifying the number of Monte Carlo iterations for each u-allocation. Higher values increase accuracy but require more computation time. Default is 5000. |
transform.fun |
A user-defined function that takes a contingency table and returns a numeric test statistic. This function should be permutation-invariant within treatment groups. |
verbose |
Logical flag indicating whether to print progress messages showing current u-allocation and estimated probabilities. Default is FALSE. |
u_space |
Optional matrix where each row represents a candidate u-allocation. If NULL, a default set of corner allocations is generated. Each row should have length equal to the number of outcomes. |
This function performs Monte Carlo sensitivity analysis for arbitrary test
statistics that are permutation-invariant. Unlike
the score test version, this function can handle any user-defined test
statistic through the transform.fun parameter.
The algorithm:
Generates tables using Sequential Importance Sampling (SIS)
Evaluates the test statistic on each sampled table
Estimates p-values using importance weights
Searches over u-allocations to find the maximum p-value
When u_space is not provided, the function generates default corner
allocations that explore extreme cases in the unmeasured confounder space.
A list containing:
rct.prob |
Estimated probability under RCT (all u-allocations zero) |
max.prob |
Maximum estimated probability across all u-allocations |
maximizer |
The u-allocation vector yielding max.prob |
obs.stat |
Observed test statistic value from transform.fun |
obs.table |
The input observed table |
gamma |
Extracted gamma value from the generic bias model |
delta |
The delta vector |
Eisinger, R. D., & Chen, Y. (2017). Sampling for Conditional Inference on Contingency Tables. Journal of Computational and Graphical Statistics, 26(1), 79–87.
exact.general.sen.IxJ for exact computation when feasible,
sampling.score.sen.IxJ for score tests with ordinal data,
# Example with custom test statistic emphasizing corner cells obs.table <- matrix(c(10, 5, 8, 12), ncol = 2, byrow = TRUE) # Define test statistic: sum of diagonal elements diag_stat <- function(tb) { sum(diag(tb)) } # Run sensitivity analysis result <- sampling.general.sen.IxJ( obs.table = obs.table, gamma = 0.5, delta = c(0, 1), transform.fun = diag_stat, mc.iteration = 5000, verbose = TRUE ) # Custom u-space example custom_u <- matrix(c(0, 0, 5, 0, 0, 8, 5, 8), ncol = 2, byrow = TRUE) result_custom <- sampling.general.sen.IxJ( obs.table = obs.table, gamma = 0.5, delta = c(0, 1), transform.fun = diag_stat, u_space = custom_u )# Example with custom test statistic emphasizing corner cells obs.table <- matrix(c(10, 5, 8, 12), ncol = 2, byrow = TRUE) # Define test statistic: sum of diagonal elements diag_stat <- function(tb) { sum(diag(tb)) } # Run sensitivity analysis result <- sampling.general.sen.IxJ( obs.table = obs.table, gamma = 0.5, delta = c(0, 1), transform.fun = diag_stat, mc.iteration = 5000, verbose = TRUE ) # Custom u-space example custom_u <- matrix(c(0, 0, 5, 0, 0, 8, 5, 8), ncol = 2, byrow = TRUE) result_custom <- sampling.general.sen.IxJ( obs.table = obs.table, gamma = 0.5, delta = c(0, 1), transform.fun = diag_stat, u_space = custom_u )
This function implements a sampling-based approach to sensitivity analysis for score tests in I by J contingency tables under the generic bias model. It uses Sequential Importance Sampling (SIS) from Eisinger and Chen (2017) to approximate p-values when exact computation is infeasible due to large reference set
sampling.score.sen.IxJ( obs.table, gamma, delta, row = "treatment", mc.iteration = 5000, treatment.scores, outcome.scores, verbose = FALSE )sampling.score.sen.IxJ( obs.table, gamma, delta, row = "treatment", mc.iteration = 5000, treatment.scores, outcome.scores, verbose = FALSE )
obs.table |
A matrix or table object representing the observed contingency table. |
gamma |
a nonnegative scalar |
delta |
A binary vector with no more than two unique values, corresponding to treatment levels. The length must match the number of treatments. Must have a monotone trend for ordinal tests. |
row |
A string indicating whether rows represent "outcome" or "treatment". Must be either "outcome" or "treatment". Default is "treatment". |
mc.iteration |
Integer specifying the number of Monte Carlo iterations for each u-allocation. Higher values increase accuracy but require more computation time. Default is 5000. |
treatment.scores |
A numeric vector specifying scores for each treatment level. Must have the same length as the number of treatments and exhibit a monotone trend. |
outcome.scores |
A numeric vector specifying scores for each outcome level. Must have the same length as the number of outcomes and exhibit a monotone trend. |
verbose |
Logical flag indicating whether to print progress messages during computation. Default is FALSE. |
The function uses importance sampling to estimate p-values for score tests when both treatments and outcomes are ordinal. The score test statistic is computed as the sum of products of cell counts with their corresponding treatment and outcome scores.
The u-space is automatically constructed to respect the ordinal nature of the data, focusing on allocations that assign bias to higher outcome levels first (assuming an increasing trend in outcome scores).
Unlike exact methods, this function provides Monte Carlo estimates that converge to true values as mc.iteration increases. Results may vary slightly between runs unless the random seed is fixed.
A list containing:
rct.prob |
Estimated probability under RCT (all u-allocations zero) |
max.prob |
Maximum estimated probability across all u-allocations |
maximizer |
The u-allocation vector yielding max.prob |
obs.stat |
Observed test statistic value |
obs.table |
The input observed table |
gamma |
Extracted gamma value from the generic bias model |
delta |
delta vector |
Eisinger, R. D., & Chen, Y. (2017). Sampling for Conditional Inference on Contingency Tables. Journal of Computational and Graphical Statistics, 26(1), 79–87.
exact.score.sen.IxJ for exact computation when feasible,
sampling.general.sen.IxJ for general test statistics.
# Binary outcome example obs.table <- matrix(c(15, 10, 5, 20), ncol = 2, byrow = TRUE) result <- sampling.score.sen.IxJ(obs.table = obs.table, gamma = 0.5, delta = c(0, 1), treatment.scores = c(0, 1), outcome.scores = c(0, 1), mc.iteration = 5000)# Binary outcome example obs.table <- matrix(c(15, 10, 5, 20), ncol = 2, byrow = TRUE) result <- sampling.score.sen.IxJ(obs.table = obs.table, gamma = 0.5, delta = c(0, 1), treatment.scores = c(0, 1), outcome.scores = c(0, 1), mc.iteration = 5000)