Score Shift Sweep: Every Element vs Rest (One-vs-Rest Wilcoxon)
Source:R/05_inferential.R
av_calc_shift.RdSweeps every element of group_by, testing each element's score
distribution against the pooled rest (the all-element generalization of
av_compare_target). Reports a tidy table ranked by shift, with
BH correction across elements and, for >= 3 elements, a Kruskal-Wallis omnibus.
One-vs-rest is a marginal screen ("which element stands out vs the baseline?");
for a specific pairwise contrast use av_compare_groups.
Usage
av_calc_shift(
obj,
group_by = "output_type",
abs_score = TRUE,
omnibus = TRUE,
min_n = 5L,
balance_by = NULL,
per = "variant_id"
)Arguments
- obj
An
AlphaVarSetobject.- group_by
Character. Column whose elements are swept (default
"output_type").- abs_score
Logical. If TRUE (default), uses absolute scores.
- omnibus
Logical. If TRUE (default), adds a Kruskal-Wallis omnibus across all elements when there are >= 3 of them.
- min_n
Minimum n per arm (default 5). Elements below it are dropped.
- balance_by
Character or NULL. When set, collapses each arm to one balanced value per unit
per(two-step weighted mean overbalance_by) before each element's test. Must differ fromperandgroup_by.- per
Character. Unit (replication) unit for
balance_by; default"variant_id". When aperunit straddles several elements ofgroup_by(a per-track column), the sweep uses paired Wilcoxon signed-rank tests plus a Friedman omnibus blocked byper; otherwise Mann-Whitney plus a balanced Kruskal-Wallis omnibus. Ignored whenbalance_byis NULL.
Value
The updated AlphaVarSet with results at
obj$stats$comparisons[[paste0(group_by, "_vs_rest")]]: a list with
method, global_p (Kruskal-Wallis or Friedman p, or
NA_real_ when no omnibus), global_effect /
global_effect_type (epsilon-squared or Kendall's W),
results_tidy (one row per kept element, schema as in
av_compare_target), params, dropped_groups.
See also
av_compare_target for a single element.
Examples
example_file <- system.file("extdata", "bcl11a_rs1427407_example.csv", package = "AlphaVaR")
obj <- av_import_csv(example_file)
#> Reading bcl11a_rs1427407_example.csv ...
obj <- av_calc_shift(obj, group_by = "output_type")
#> Warning: 1 element(s) dropped (n < 5 in one arm): SPLICE_SITES
res <- obj$stats$comparisons[["output_type_vs_rest"]]
res$method
#> [1] "Kruskal-Wallis Omnibus + Wilcoxon One-vs-Rest (BH)"
res$global_p
#> [1] 4.609982e-130
res$results_tidy
#> # A tibble: 10 × 12
#> group n_group n_rest median_group median_rest median_diff p_val p_adj
#> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 ATAC 320 2356 0.814 0.454 0.360 0 0
#> 2 CAGE 320 2356 0.696 0.465 0.231 0 0
#> 3 SPLICE… 320 2356 0.277 0.532 -0.254 7.46e-43 2.49e-42
#> 4 RNA_SEQ 320 2356 0.270 0.517 -0.247 2.29e-26 5.72e-26
#> 5 SPLICE… 320 2356 0.336 0.517 -0.181 2.71e-14 5.43e-14
#> 6 CHIP_H… 320 2356 0.640 0.465 0.175 5.02e-10 8.37e-10
#> 7 DNASE 320 2356 0.529 0.483 0.0461 3.05e- 1 4.36e- 1
#> 8 CONTAC… 56 2620 0.535 0.488 0.0461 5.30e- 1 6.63e- 1
#> 9 CHIP_TF 328 2348 0.491 0.494 -0.00288 7.21e- 1 8.01e- 1
#> 10 PROCAP 48 2628 0.565 0.494 0.0705 8.50e- 1 8.50e- 1
#> # ℹ 4 more variables: effect_size <dbl>, effect_type <chr>, method <chr>,
#> # significance_stars <chr>
# Tracks without a biosample label cannot be blocked by biosample.
obj_bs <- av_subset(obj, !is.na(biosample_name))
obj_bs <- av_calc_shift(obj_bs, group_by = "output_type",
balance_by = "variant_id", per = "biosample_name")
#> Warning: Friedman omnibus: dropped 620 of 620 'biosample_name' units with incomplete element coverage.
obj_bs$stats$comparisons[["output_type_vs_rest"]]$results_tidy
#> # A tibble: 10 × 12
#> group n_group n_rest median_group median_rest median_diff p_val p_adj
#> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 SPLICE… 177 177 0.294 0.467 -0.178 1.07e-17 1.07e-16
#> 2 CHIP_H… 126 126 0.623 0.415 0.192 1.67e-13 8.37e-13
#> 3 ATAC 69 69 0.686 0.416 0.289 2.05e-11 6.83e-11
#> 4 RNA_SEQ 172 172 0.298 0.469 -0.146 8.05e- 6 2.01e- 5
#> 5 CAGE 78 78 0.607 0.425 0.187 2.57e- 5 5.14e- 5
#> 6 SPLICE… 202 202 0.365 0.441 -0.0787 1.67e- 4 2.79e- 4
#> 7 CONTAC… 8 8 0.491 0.412 0.129 5.87e- 2 8.39e- 2
#> 8 DNASE 149 149 0.535 0.475 0.0296 2.35e- 1 2.93e- 1
#> 9 PROCAP 6 6 0.511 0.482 -0.00843 8.34e- 1 9.27e- 1
#> 10 CHIP_TF 45 45 0.506 0.497 0.0223 9.64e- 1 9.64e- 1
#> # ℹ 4 more variables: effect_size <dbl>, effect_type <chr>, method <chr>,
#> # significance_stars <chr>