---
title: "Getting started with uncertR"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with uncertR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(uncertR)
Overview

uncertR provides bootstrap-based estimates of prediction uncertainty for
gene-expression-based drug-response models. The package generates
out-of-bag (OOB) predictions from bootstrap models and summarizes the
variability of those predictions.

This vignette demonstrates the basic workflow using a small simulated
gene-expression dataset. The example is intentionally small so that it runs
quickly and does not require a large external dataset.

Simulating gene-expression and drug-response data

We first create a gene-expression matrix with genes as rows and samples as
columns. We then simulate a continuous drug-response phenotype that depends
on several genes.

set.seed(42)

n_genes <- 100
n_samples <- 80

expr <- matrix(
  rnorm(n_genes * n_samples),
  nrow = n_genes,
  ncol = n_samples,
  dimnames = list(
    paste0("Gene", seq_len(n_genes)),
    paste0("Sample", seq_len(n_samples))
  )
)

y <- 2 * expr[1, ] -
  1.5 * expr[2, ] +
  0.8 * expr[3, ] +
  rnorm(n_samples, sd = 1)

dim(expr)
length(y)
Bootstrap prediction uncertainty

boot_uncertainty() resamples the samples with replacement, fits an
elastic-net model to each bootstrap sample, and predicts only samples that
were not selected in that bootstrap replicate. These are the out-of-bag
predictions.

For this example, we use 40 bootstrap replicates to keep computation
relatively short.

fit <- boot_uncertainty(
  expr = expr,
  y = y,
  B = 40,
  seed = 42
)

fit

The OOB prediction matrix contains one column for each bootstrap replicate.
A sample has an NA for a replicate when it was included in that replicate's
bootstrap training set.

dim(fit$oob)
head(fit$oob[, 1:5])
Calculating uncertainty scores

uncertainty_score() summarizes the variability of the OOB predictions for
each sample. Three metrics are available: standard deviation ("sd"),
interquartile range ("iqr"), and percentile confidence-interval width
("ci_width").

unc_sd <- uncertainty_score(
  fit,
  metric = "sd",
  min_oob = 5
)

head(unc_sd$scores)

The same OOB predictions can be summarized using the other metrics:

unc_iqr <- uncertainty_score(
  fit,
  metric = "iqr",
  min_oob = 5
)

unc_ci <- uncertainty_score(
  fit,
  metric = "ci_width",
  min_oob = 5
)

head(unc_iqr$scores)
head(unc_ci$scores)

The resulting unc_scores object contains the mean OOB prediction,
uncertainty score, observed response, absolute prediction error, and number
of OOB predictions for each sample.

Validating uncertainty

Prediction uncertainty and prediction error are distinct quantities.
uncertainty_score() measures the dispersion of the OOB predictions, while
validate_uncertainty() evaluates whether estimated uncertainty is
associated with absolute OOB prediction error.

validation <- validate_uncertainty(unc_sd)

validation

The result includes a Spearman correlation between uncertainty and absolute
prediction error, its p-value, and mean prediction error across uncertainty
quartiles.

Investigating drivers of uncertainty

uncertainty_drivers() can be used to examine whether sample-level
covariates explain variation in prediction uncertainty.

The function fits a joint linear model and reports the partial R-squared of
each covariate, permutation p-values, and Benjamini-Hochberg false discovery
rates (FDR). By default, uncertainty is adjusted for the predicted value
using a spline because uncertainty can be mechanically related to prediction
extremeness.

Here, we create two example covariates.

set.seed(42)

covariates <- data.frame(
  purity = runif(n_samples),
  technical_factor = rnorm(n_samples),
  row.names = rownames(unc_sd$scores)
)

drivers <- uncertainty_drivers(
  unc_sd,
  covariates,
  adjust_for = "prediction",
  n_perm = 100,
  seed = 42
)

drivers

The resulting table can be accessed directly:

drivers$table

For numeric covariates, spearman_rho reports the Spearman correlation
between the covariate and the adjusted uncertainty. The partial_r2
measures the contribution of each covariate in the joint model, while p_perm
and fdr provide permutation-based significance measures.

Preparing additional covariates

tumor_covariates() can combine technical expression-derived covariates with
user-supplied sample-level covariates.

For example:

extra <- data.frame(
  tissue = sample(
    c("Tissue_A", "Tissue_B", "Tissue_C"),
    n_samples,
    replace = TRUE
  ),
  growth_rate = rnorm(n_samples),
  row.names = colnames(expr)
)

covariate_data <- tumor_covariates(
  expr,
  extra = extra,
  technical = TRUE
)

head(covariate_data)

The technical covariates include log library size and the number of detected
genes. These measures are most appropriate for raw or TPM-like RNA-seq data.
For RMA-normalized microarray data such as GDSC, technical = FALSE should
generally be used.

Optional immune-cell deconvolution methods can also be supplied through the
methods argument when immunedeconv is installed. These methods are
intended for patient tumor expression data and are not generally meaningful
for cell-line data.

Transforming compositional fractions

Cell-type fractions are compositional because their components sum to a
constant. clr_transform() applies a centered log-ratio transformation.

fractions <- matrix(
  c(
    0.2, 0.3, 0.5,
    0.1, 0.1, 0.8
  ),
  nrow = 2,
  byrow = TRUE,
  dimnames = list(
    paste0("Sample", 1:2),
    c("CellTypeA", "CellTypeB", "CellTypeC")
  )
)

clr_transform(fractions)

A small pseudocount is added before taking logarithms so that zero fractions
can be handled.

Collapsing rare factor levels

Categorical covariates with many small groups can make downstream models
unstable. collapse_rare_levels() combines levels with fewer than the
specified number of observations into "Other".

tissue <- c(
  rep("Tissue_A", 30),
  rep("Tissue_B", 30),
  rep("Rare", 5),
  rep("Rare_2", 5)
)

collapse_rare_levels(tissue, min_n = 10)
Summary

A typical uncertR workflow is:

Provide a gene-expression matrix and continuous response variable.
Generate bootstrap models and OOB predictions with
boot_uncertainty().
Summarize prediction uncertainty with uncertainty_score().
Assess whether uncertainty tracks prediction error with
validate_uncertainty().
Investigate biological or technical drivers of uncertainty with
uncertainty_drivers().
Use tumor_covariates(), clr_transform(), and
collapse_rare_levels() to prepare covariates for downstream analysis.

The vignette uses simulated data so that it remains lightweight,
reproducible, and suitable for automated package checks without requiring
large external datasets.