Package {hetNSC}

hetNSC: Heteroscedastic Nearest Shrunken Centroid Classification

Description

Nearest shrunken centroid classifiers for high-dimensional data, with either the pooled-variance (homoscedastic) test statistic and discriminant function of PAM or the heteroscedastic (class-specific variance) ones of Sahtout and Wang, each combined with soft, hard or order thresholding.

The six classifiers

* **STh** (soft, homoscedastic): the original PAM classifier of Tibshirani et al. (2002). * **HTh**, **OTh** (hard / order, homoscedastic): the homoscedastic classifiers of Sahtout, Wang and Ghimire (2024), which replace PAM's soft thresholding with hard or order thresholding. * **STh3**, **HTh3**, **OTh3**: the heteroscedastic classifiers of Sahtout and Wang.

So 'variance = "homo"' is the original PAM only when 'method = "soft"'.

Main functions

* [hnsc()]: fit a classifier, choosing the threshold by refined cross-validation; then [predict.hnsc()]. * [hnsc_train()], [hnsc_cv()], [hnsc_select()]: lower-level steps. * [compare_nsc()]: the repeated-run comparison protocol of the paper; [summary.hnsc_comparison()] and [recall_table()] summarise it. * [simulate_hetero_data()]: the correlated-gene simulation designs. * [compute_rho2()]: the bias-variance crossover diagnostic of Section 2.4. * [make_cm()], [per_class_recall()], [balanced_acc()]: class-level metrics.

Data layout

As in pamr, expression matrices have **genes in rows and samples in columns**.

Author(s)

Maintainer: Haiyan Wang hwang@ksu.edu

Authors:

• Mohammad Omar Sahtout

Other contributors:

• Trevor Hastie (Author of the 'pamr' code this package is derived from) [contributor]

• Robert Tibshirani (Author of the 'pamr' code this package is derived from) [contributor]

• Balasubramanian Narasimhan (Author of the 'pamr' code this package is derived from) [contributor]

• Gilbert Chu (Author of the 'pamr' code this package is derived from) [contributor]

Class-balanced cross-validation folds

Description

Splits sample indices into folds so that every class is spread as evenly as possible across folds. This is the fold generator of the pamr package, which is also used inside [hnsc_cv()].

Usage

balanced_folds(y, nfolds = min(min(table(y)), 10))

Arguments

Value

A list of integer vectors of sample indices, one per fold.

Examples

set.seed(1)
balanced_folds(rep(c("A", "B"), c(6, 9)))

Compare homoscedastic and heteroscedastic NSC classifiers on one dataset

Description

Reproduces the evaluation protocol of the paper: for each thresholding rule, the homoscedastic classifier (STh/HTh/OTh) and its heteroscedastic counterpart (STh3/HTh3/OTh3) are fitted 'n.runs' times, each time with the threshold chosen by [hnsc_select()] after a random reordering of the training samples (which changes the cross-validation folds). Each fit is evaluated on the test set.

Usage

compare_nsc(
  x,
  y,
  xtest,
  ytest,
  n.runs = 100,
  thresholding = c("hard", "soft", "order"),
  seed_fun = function(i) 1000 + 10 * i,
  file = NULL,
  verbose = TRUE
)

Arguments

Details

The training samples are re-shuffled cumulatively across runs and thresholding rules, exactly as in the script used for the paper, so that 'compare_nsc()' reproduces those numbers.

Value

A list of class '"hnsc_comparison"' with * 'error': mean test error rate of each method, * 'ave_num_genes': mean number of genes used, * 'errors', 'n_genes': the per-run values (runs x methods matrices), * 'cm_list': per method, a list of per-run confusion matrices, * 'classes', 'n' (training sample size) and 'K'.

Use [summary()] on the result for a table that includes balanced accuracy, and [recall_table()] for per-class recall.

Examples

set.seed(8900)
dat <- simulate_hetero_data(ntr1 = 15, ntr2 = 20, nte1 = 50, nte2 = 50,
                            p = 300, a = 15)
res <- compare_nsc(dat$xtr, dat$ytr, dat$xte, dat$yte, n.runs = 2,
                   verbose = FALSE)
summary(res)

Empirical relative heterogeneity rho_k^2 and the crossover condition

Description

Computes, for each class k, the empirical relative heterogeneity rho_k^2 = (sigma_k^2 / sigma_bar^2 - 1)^2 of Section 2.4 of the paper, where sigma_bar^2 = sum_k n_k sigma_k^2 / n, summarised across genes, and compares it with the simple crossover bound 2 / (n_k - 1) (obtained with w_k = 1 in eq. 2.11). When 'crossover_met' is 'TRUE' the class-specific variance estimate is expected to have lower mean squared error than the pooled one, i.e. the heteroscedastic classifiers are expected to help.

Usage

compute_rho2(
  X_train,
  y_train,
  X_test = X_train,
  y_test = y_train,
  summary_fn = stats::median
)

Arguments

Value

A data frame with one row per class: 'class', 'n_k', 'sigma2_hat', 'sigma_bar2', 'rho2_hat', 'threshold' (= 2/(n_k - 1)) and 'crossover_met'.

Examples

set.seed(1)
dat <- simulate_hetero_data(ntr1 = 15, ntr2 = 20, nte1 = 5, nte2 = 5,
                            p = 500, a = 12)
compute_rho2(dat$xtr, dat$ytr)

Fit a nearest shrunken centroid classifier with the threshold chosen by cross-validation

Description

The main user-facing function. It selects the threshold with the refined cross-validation search of [hnsc_select()] (unless 'threshold' is given) and then fits the classifier on all training data at that threshold. Use [predict.hnsc()] to classify new samples.

Usage

hnsc(
  x,
  y,
  method = c("hard", "soft", "order"),
  variance = c("hetero", "homo"),
  threshold = NULL,
  offset.percent = 50,
  prior = NULL,
  max.iter = 100
)

Arguments

Value

An object of class '"hnsc"' (see [hnsc_train()]) fitted at the selected threshold, with an extra component 'selection' holding the [hnsc_select()] result (or 'NULL' if 'threshold' was supplied) and 'selected.genes', the row indices of the genes used by the classifier.

References

Sahtout, M. O. and Wang, H. Heteroscedastic nearest shrunken centroid classification: methods and applications to high-dimensional data. *Statistical Applications in Genetics and Molecular Biology* (accepted).

Sahtout, M. O., Wang, H. and Ghimire, S. (2024). Different thresholding methods on nearest shrunken centroid algorithm. *Communications in Statistics - Simulation and Computation*, 53(3), 1444-1460. doi:10.1080/03610918.2022.2047201

Tibshirani, R., Hastie, T., Narasimhan, B. and Chu, G. (2002). Diagnosis of multiple cancer types by shrunken centroids of gene expression. *PNAS*, 99(10), 6567-6572.

Examples

set.seed(1)
dat <- simulate_hetero_data(ntr1 = 15, ntr2 = 20, nte1 = 50, nte2 = 50,
                            p = 500, a = 15)
fit <- hnsc(dat$xtr, dat$ytr, method = "HTh3")
fit
pred <- predict(fit, dat$xte)
mean(pred != dat$yte)            # test error
head(predict(fit, dat$xte, type = "posterior"))

Cross-validation for a fitted NSC classifier

Description

Estimates the misclassification error at each threshold of a fitted [hnsc_train()] object by class-balanced K-fold cross-validation (the procedure of pamr's 'pamr.cv'). The number of folds is the size of the smallest class, capped at 10.

Usage

hnsc_cv(fit, x, y, folds = NULL)

Arguments

Value

A list of class '"hnsc_cv"' with components 'threshold', 'size' (number of genes kept by the full-data fit), 'error' (cross-validated error rate), 'yhat' (matrix of cross-validated predictions) and 'folds'.

Examples

set.seed(1)
dat <- simulate_hetero_data(ntr1 = 10, ntr2 = 10, nte1 = 5, nte2 = 5,
                            p = 200, a = 15)
fit <- hnsc_train(dat$xtr, dat$ytr, method = "HTh3")
cv <- hnsc_cv(fit, dat$xtr, dat$ytr)
cv

The six nearest shrunken centroid classifiers

Description

Lists the six classifiers implemented in the package, with the labels used in the paper and the corresponding 'method' / 'variance' arguments.

Usage

hnsc_methods()

Details

Only STh (soft thresholding with pooled variance) is the original PAM classifier; HTh and OTh are the homoscedastic classifiers of Sahtout, Wang and Ghimire (2024).

Value

A data frame with columns 'label', 'method', 'variance', 'description' and 'reference'.

Examples

hnsc_methods()

Choose the threshold by refined ("deep") cross-validation

Description

Implements the threshold search of Sahtout, Wang and Ghimire (2024) used in the paper: a cross-validation over the default threshold grid, followed by repeated refinement of the grid between the neighbours of the current best threshold until the number of selected genes can no longer change. Among thresholds with equal cross-validated error, the one examined last (for soft and hard thresholding typically the larger, sparser one) is chosen.

Usage

hnsc_select(
  x,
  y,
  method = c("hard", "soft", "order"),
  variance = c("hetero", "homo"),
  offset.percent = 50,
  max.iter = 100
)

Arguments

Value

A list of class '"hnsc_select"' with 'threshold' (selected value), 'size' (genes kept), 'error' (its CV error), 'path' (matrix of every threshold / size / CV error evaluated), 'method', 'variance' and 'label'.

Examples

set.seed(1)
dat <- simulate_hetero_data(ntr1 = 10, ntr2 = 10, nte1 = 5, nte2 = 5,
                            p = 200, a = 15)
sel <- hnsc_select(dat$xtr, dat$ytr, method = "OTh3")
sel$threshold

Train a nearest shrunken centroid classifier over a grid of thresholds

Description

Fits one of the six nearest shrunken centroid (NSC) classifiers for a sequence of threshold values. 'variance = "homo"' uses the PAM test statistic and discriminant function with a pooled within-class standard deviation; combined with 'method = "soft"' this is the original PAM classifier (STh, Tibshirani et al. 2002), while 'method = "hard"' or '"order"' gives the homoscedastic classifiers of Sahtout, Wang and Ghimire (2024) (HTh, OTh). 'variance = "hetero"' gives the heteroscedastic statistic of Sahtout and Wang with class-specific variances in both the test statistic (eq. 2.7) and the discriminant function (eq. 2.9) (STh3, HTh3, OTh3).

Usage

hnsc_train(
  x,
  y,
  method = c("hard", "soft", "order"),
  variance = c("hetero", "homo"),
  threshold = NULL,
  n.threshold = 30,
  offset.percent = 50,
  prior = NULL,
  sign.contrast = c("both", "positive", "negative"),
  remove.zeros = TRUE
)

Arguments

Details

Most users will want [hnsc()], which also chooses the threshold by cross-validation.

Value

An object of class '"hnsc"': a list containing, among others, 'threshold', 'nonzero' (number of genes kept at each threshold), 'errors' (training resubstitution errors), 'centroids', 'centroid.overall', 'sd' (denominator of d_ik), 'sik2' (s_ik^2 + s_0k^2, heteroscedastic only), 'delta' (the unshrunken test statistics d_ik), 'prior', 'method', 'variance' and 'label'.

See Also

[hnsc()], [hnsc_cv()], [predict.hnsc()]

Examples

set.seed(1)
dat <- simulate_hetero_data(ntr1 = 10, ntr2 = 10, nte1 = 5, nte2 = 5,
                            p = 200, a = 15)
fit <- hnsc_train(dat$xtr, dat$ytr, method = "hard", variance = "hetero")
fit

Confusion matrix, per-class recall and balanced accuracy

Description

'make_cm()' builds a K x K confusion matrix (rows = true class, columns = predicted class). 'per_class_recall()' returns the diagonal divided by the row totals (sensitivity of each class) and 'balanced_acc()' their mean (macro-averaged recall).

Usage

make_cm(
  true_labels,
  pred_labels,
  classes = sort(unique(as.character(true_labels)))
)

per_class_recall(cm)

balanced_acc(cm)

Arguments

Value

'make_cm()': an integer matrix; 'per_class_recall()': a named numeric vector; 'balanced_acc()': a number.

Examples

cm <- make_cm(c("A", "A", "B", "B"), c("A", "B", "B", "B"), c("A", "B"))
cm
per_class_recall(cm)
balanced_acc(cm)

Predict with a nearest shrunken centroid classifier

Description

Predict with a nearest shrunken centroid classifier

Usage

## S3 method for class 'hnsc'
predict(
  object,
  newx = NULL,
  threshold = NULL,
  type = c("class", "posterior", "centroid", "nonzero"),
  ...
)

Arguments

Value

A factor, matrix or integer vector depending on 'type'.

Examples

set.seed(1)
dat <- simulate_hetero_data(ntr1 = 10, ntr2 = 10, nte1 = 5, nte2 = 5,
                            p = 200, a = 15)
fit <- hnsc_train(dat$xtr, dat$ytr, method = "STh3")
predict(fit, dat$xte, threshold = fit$threshold[5])

Simulate two-class data with correlated informative genes

Description

Generates the two-class designs with correlated informative genes used in the paper's simulation study. Class A has mean zero for every gene; class B has mean 'a/10 * i' for the first 'ngene' (informative) genes and zero elsewhere. All non-informative genes are independent with standard deviation 's[i]' in both classes.

Usage

simulate_hetero_data(
  ntr1,
  ntr2,
  nte1,
  nte2,
  ngene = 20,
  p = 10000,
  s = rep_len(1:20, p),
  a = 10,
  rho = 0.95,
  correlated = c("B", "both")
)

Arguments

Details

The informative genes are drawn jointly from a multivariate normal distribution with covariance

\Sigma = v v^\top, \qquad v_i = \sqrt{\rho^{\,i-1}}\, s_i,

i.e. \Sigma_{ij} = \rho^{(i+j-2)/2} s_i s_j. This is the matrix used to produce the results in the paper. It has rank one, so the informative genes are perfectly correlated with one another, and gene 'i' has variance \rho^{i-1} s_i^2.

Random numbers are drawn in the same order as in the original scripts, so with the same seed the data are identical (after transposition) to those used for the paper.

Value

A list with 'xtr', 'xte' (genes x samples matrices), 'ytr', 'yte' (character label vectors with values '"A"' and '"B"').

Examples

set.seed(8900)
dat <- simulate_hetero_data(ntr1 = 15, ntr2 = 20, nte1 = 50, nte2 = 50,
                            p = 1000, a = 12)
dim(dat$xtr)
table(dat$ytr)

Summaries of a classifier comparison

Description

'summary()' gives, for each method, the mean and standard deviation of the test error over runs, the balanced accuracy and the mean number of genes used. 'recall_table()' gives the recall (sensitivity) of every class. Both are computed from the confusion matrix averaged over runs, as in the paper.

Usage

## S3 method for class 'hnsc_comparison'
summary(object, ...)

recall_table(object)

Arguments

Value

A data frame. 'summary()': one row per method. 'recall_table()': one row per class and one 'recall_<method>' column per method.

Examples

set.seed(8900)
dat <- simulate_hetero_data(ntr1 = 15, ntr2 = 20, nte1 = 50, nte2 = 50,
                            p = 300, a = 15)
res <- compare_nsc(dat$xtr, dat$ytr, dat$xte, dat$yte, n.runs = 2,
                   thresholding = "hard", verbose = FALSE)
summary(res)
recall_table(res)