Skip to contents

Compute Pseudo-Correlation Using Beta-Binomial Model

Source: R/general_tools.R
get_pseudo_correlation.Rd

This function calculates a pseudo R²-like correlation metric using a beta-binomial model implemented in C++. It takes in a data matrix ZDB_matrix and two model matrices for inclusion and exclusion, respectively. The function now supports both sparse and dense matrices for m1 and m2, and allows selection between Cox-Snell and Nagelkerke R² metrics.

Usage

get_pseudo_correlation(
  ZDB_matrix,
  m1_inclusion = NULL,
  m2_exclusion = NULL,
  metric = "CoxSnell",
  suppress_warnings = TRUE,
  verbose = FALSE,
  permutation_count = 100L,
  permutation_seed = NULL
)

Arguments

ZDB_matrix

A numeric dense matrix of shape (events x samples). Should have rownames representing events.

m1_inclusion

A numeric matrix (dense or sparse) of the same number of rows as ZDB_matrix, representing inclusion features.

m2_exclusion

A numeric matrix (dense or sparse) of the same number of rows as ZDB_matrix, representing exclusion features.

metric

Character string specifying which R² metric to compute. Options are "CoxSnell" (default) or "Nagelkerke".

suppress_warnings

Logical. If TRUE (default), suppresses warnings during computation (e.g., due to ill-conditioned inputs).

verbose

Logical. If TRUE, prints progress and informational messages. Default is FALSE.

permutation_count

Integer. Number of times the null model is generated by permuting the cells (columns) of ZDB_matrix and recomputing the pseudo correlation. Each permutation yields one null value per event, so the n draws form a per-event empirical null. Defaults to 100L. Use 1L to reproduce the original single-permutation behaviour. Values of 1000 or more are recommended for reliable empirical FDR estimation (the smallest achievable empirical p-value is approximately 1 / (permutation_count + 1)).

permutation_seed

Optional single numeric value. If supplied, the random number generator is seeded with it before the permutations are drawn, making the null distribution (and therefore the empirical p-values) reproducible. The caller's global RNG stream is saved and restored, so passing a seed does not perturb downstream randomness. Defaults to NULL (use the ambient RNG state).

Value

An S3 object of class splikit_pseudo_correlation_result with three components:

statistics

A per-event data.table with the columns described below.

null_distribution

A long data.table containing one row per retained event and permutation, with columns event, permutation, and null_pseudo_correlation. This table is suitable for direct export or pooled descriptive analyses; event-level p-values remain based on each event's own permutation draws.

metadata

A list recording the permutation settings and null-draw counts.

The statistics component contains:

event

The event names from ZDB_matrix rownames.

pseudo_correlation

The computed pseudo R² correlation values using the specified metric.

null_distribution

Mean of the per-event null correlation values across the permutation_count permutations (equal to the single null draw when permutation_count = 1).

null_sd

Standard deviation of the per-event null draws (NA when fewer than two valid draws are available).

n_perm_valid

Number of permutations that produced a non-NA null value for the event.

emp_pvalue

Two-sided empirical p-value, (b + 1) / (n_perm_valid + 1), where b is the number of valid null draws whose absolute value is greater than or equal to the absolute observed correlation.

emp_padj

Benjamini-Hochberg adjusted emp_pvalue across the retained events.

Examples

# \donttest{
set.seed(42)
# get the m1 object
junction_abundance_object <- load_toy_SJ_object()
m1_obj <- make_m1(junction_ab_object = junction_abundance_object)

# obtaining the m1 and eventdata
m1_inclusion <- m1_obj$m1_inclusion_matrix
eventdata <- m1_obj$event_data
m2_exclusion <- make_m2(m1_inclusion, eventdata)

# subset to a few hundred cells to keep the example quick
cells <- seq_len(min(250, ncol(m1_inclusion)))
m1_inclusion <- m1_inclusion[, cells]
m2_exclusion <- m2_exclusion[, cells]

# creating a dummy ZDB
ZDB_matrix <- matrix(rnorm(n = (nrow(m1_inclusion) * ncol(m1_inclusion)), sd = 7),
nrow = nrow(m1_inclusion),
ncol = ncol(m1_inclusion))
rownames(ZDB_matrix) <- rownames(m1_inclusion)

# m1 and m2 can now be either sparse or dense matrices
# Example with dense matrices (backward compatible)
m1_dense <- as.matrix(m1_inclusion)
m2_dense <- as.matrix(m2_exclusion)
pseudo_r_square_cox <- get_pseudo_correlation(ZDB_matrix, m1_dense, m2_dense,
                                              permutation_count = 20,
                                              permutation_seed = 1)
print(pseudo_r_square_cox)
#> Pseudo-correlation result
#> Statistics: 2 retained event(s)
#>                        event pseudo_correlation null_distribution   null_sd
#>                       <char>              <num>             <num>     <num>
#> 1: chr19:27376501-27377078_E        -0.04460012        0.01023085 0.1169676
#> 2: chr19:27376505-27377078_E         0.12066213        0.04112645 0.1044667
#>    n_perm_valid emp_pvalue  emp_padj
#>           <int>      <num>     <num>
#> 1:           20  0.7619048 0.7619048
#> 2:           20  0.2857143 0.5714286
#> Null distribution: 40 row(s); 40 valid draw(s)

# Example with sparse matrices (more memory efficient)
pseudo_r_square_sparse <- get_pseudo_correlation(ZDB_matrix, m1_inclusion, m2_exclusion,
                                                 permutation_count = 20,
                                                 permutation_seed = 1)

# Per-event empirical p-values and the exportable long null distribution
head(pseudo_r_square_sparse$statistics[,
  c("event", "pseudo_correlation", "emp_pvalue")])
#>                        event pseudo_correlation emp_pvalue
#>                       <char>              <num>      <num>
#> 1: chr19:27376501-27377078_E        -0.04460012  0.7619048
#> 2: chr19:27376505-27377078_E         0.12066213  0.2857143
dim(pseudo_r_square_sparse$null_distribution)
#> [1] 40  3

# Example using Nagelkerke R-squared instead of Cox-Snell
pseudo_r_square_nagel <- get_pseudo_correlation(ZDB_matrix, m1_inclusion, m2_exclusion,
                                                metric = "Nagelkerke",
                                                permutation_count = 20)
print(pseudo_r_square_nagel)
#> Pseudo-correlation result
#> Statistics: 2 retained event(s)
#>                        event pseudo_correlation null_distribution   null_sd
#>                       <char>              <num>             <num>     <num>
#> 1: chr19:27376501-27377078_E         -0.1349213       -0.02549244 0.2603528
#> 2: chr19:27376505-27377078_E          0.3650189        0.08476141 0.2470716
#>    n_perm_valid emp_pvalue  emp_padj
#>           <int>      <num>     <num>
#> 1:           20  0.8095238 0.8095238
#> 2:           20  0.1904762 0.3809524
#> Null distribution: 40 row(s); 40 valid draw(s)
# }