## ----confounder-setup, include = FALSE----------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")

## ----build-cohort-------------------------------------------------------------
library(rnaSentry)
library(SummarizedExperiment)
library(survival)

set.seed(12)
n <- 180
latent <- rnorm(n)
region <- factor(cut(latent, breaks = 3, labels = c("R1", "R2", "R3")))
batch <- as.character(region)
flip <- runif(n) < 0.2
batch[flip] <- sample(c("R1", "R2", "R3"), sum(flip), replace = TRUE)
batch <- factor(batch)

counts <- matrix(rpois(50 * n, lambda = 200), nrow = 50, ncol = n,
                 dimnames = list(paste0("g", seq_len(50)), paste0("S", seq_len(n))))
counts[1:6, ] <- counts[1:6, ] *
  matrix(1 + 0.6 * as.integer(batch), nrow = 6, ncol = n, byrow = TRUE)

event_time <- rexp(n, rate = 0.05 * exp(0.8 * latent))
censor_time <- rexp(n, rate = 0.03)
cd <- S4Vectors::DataFrame(batch = batch, region = region,
                           time = pmin(event_time, censor_time),
                           event = as.integer(event_time < censor_time),
                           row.names = colnames(counts))
se <- SummarizedExperiment(assays = list(counts = counts), colData = cd)
se

## ----naive-screen, warning = FALSE--------------------------------------------
lg <- log2(counts + 1)
stat <- do.call(rbind, lapply(rownames(lg), function(g) {
    s <- summary(coxph(Surv(cd$time, cd$event) ~ lg[g, ]))
    c(beta = s$coefficients[1, 1], p = s$coefficients[1, 5])
}))
rownames(stat) <- rownames(lg)
naive_top <- head(stat[order(stat[, "p"]), , drop = FALSE], 6)
data.frame(gene = rownames(naive_top),
           planted_batch_effect = rownames(naive_top) %in% paste0("g", 1:6),
           naive_top)

## ----double-count-------------------------------------------------------------
summary(coxph(Surv(time, event) ~ batch + region, data = as.data.frame(cd)))

## ----design-audit-------------------------------------------------------------
da <- design_audit(se, design_vars = c("batch", "region"),
                   redundant_effect_size = 0.5)

## ----pairwise-table-----------------------------------------------------------
da$pairwise_table

## ----confounder-table---------------------------------------------------------
da$confounder_table

## ----pca-batch-tests----------------------------------------------------------
pa <- pca_audit(se, batch_col = "batch")
pa$batch_tests

## ----pca-batch-plot, fig.width = 7, fig.height = 5, fig.align = "center", fig.cap = "PC1-PC2 scatter of the batch-confounded cohort, colored by batch. The three batch clusters separate along PC1 (the surrogate gradient); this is the structure the confounder scan detects."----
plot_pca_audit(pa, color_by = "batch")

## ----pca-region-plot, fig.width = 7, fig.height = 5, fig.align = "center", fig.cap = "The same PC1-PC2 scatter colored by region. The near-identical structure is why the redundancy scan reports Cramer's V = 0.82 and the recommended formula drops `region`."----
plot_pca_audit(pa, color_by = "region")

## ----recovered-design---------------------------------------------------------
da$formula_text

## ----adj-screen, warning = FALSE----------------------------------------------
adj <- do.call(rbind, lapply(rownames(lg), function(g) {
    s <- summary(coxph(Surv(cd$time, cd$event) ~ lg[g, ] + batch,
                       data = as.data.frame(cd)))
    c(beta = s$coefficients[1, 1], p = s$coefficients[1, 5])
}))
rownames(adj) <- rownames(lg)
guarded_top <- head(adj[order(adj[, "p"]), , drop = FALSE], 6)
data.frame(gene = rownames(guarded_top),
           planted_batch_effect = rownames(guarded_top) %in% paste0("g", 1:6),
           guarded_top)

## ----guarded-sig, warning = TRUE----------------------------------------------
sig_guarded <- build_signature(se, "time", "event", top_n = 6,
                               repeats = 3, folds = 3,
                               design_terms = "batch", seed = 12)
sig_guarded$genes
sig_guarded$flags

## ----session-info-------------------------------------------------------------
sessionInfo()