hetNSC implements six nearest shrunken centroid (NSC)
classifiers for high-dimensional, low-sample-size data:
library(hetNSC)
hnsc_methods()[, c("label", "method", "variance")]
#> label method variance
#> 1 STh soft homo
#> 2 HTh hard homo
#> 3 OTh order homo
#> 4 STh3 soft hetero
#> 5 HTh3 hard hetero
#> 6 OTh3 order heteroAs in pamr, expression matrices have genes
in rows and samples in columns, and y has one
label per column.
set.seed(8900)
dat <- simulate_hetero_data(ntr1 = 15, ntr2 = 20, nte1 = 50, nte2 = 50,
p = 1000, a = 14)
dim(dat$xtr)
#> [1] 1000 35
table(dat$ytr)
#>
#> A B
#> 15 20simulate_hetero_data() generates the correlated-gene
designs of the paper’s simulation study: 20 informative genes with
class-B means a/10 * i, gene standard deviations cycling
through 1, …, 20, and correlated informative genes in class B only
(correlated = "B", the default, heteroscedastic) or in both
classes (correlated = "both", homoscedastic). See
?simulate_hetero_data for the exact covariance matrix.
hnsc() chooses the threshold by refined cross-validation
and fits the classifier at that threshold. Choose the classifier with a
paper label or with method and variance.
set.seed(1)
fit <- hnsc(dat$xtr, dat$ytr, method = "HTh3")
fit
#> Nearest shrunken centroid classifier: HTh3 (hard thresholding, class-specific variance)
#> Classes: A, B | genes: 1000
#> Selected threshold: 5.226 (genes used: 2, CV error: 0.0571)
pred <- predict(fit, dat$xte)
mean(pred != dat$yte) # test error
#> [1] 0.12
head(predict(fit, dat$xte, type = "posterior"))
#> A B
#> te1 0.99961637 0.0003836345
#> te2 0.03044574 0.9695542589
#> te3 0.99403579 0.0059642118
#> te4 0.99479673 0.0052032674
#> te5 0.10279039 0.8972096126
#> te6 0.59259413 0.4074058665
head(fit$selected.genes) # genes used by the classifier
#> [1] 17 19The homoscedastic counterpart is fitted the same way:
hnsc_train() fits a grid of thresholds,
hnsc_cv() cross-validates it, and
hnsc_select() runs the full refined search used by
hnsc().
tr <- hnsc_train(dat$xtr, dat$ytr, method = "OTh3")
set.seed(2)
cv <- hnsc_cv(tr, dat$xtr, dat$ytr)
head(data.frame(threshold = cv$threshold, genes = cv$size, cv.error = cv$error))
#> threshold genes cv.error
#> 1 500 250 0.1428571
#> 2 483 242 0.1428571
#> 3 466 233 0.1714286
#> 4 448 224 0.1714286
#> 5 431 216 0.1714286
#> 6 414 207 0.2000000
set.seed(2)
sel <- hnsc_select(dat$xtr, dat$ytr, method = "OTh3")
sel
#> Deep CV threshold search for OTh3
#> Selected threshold: 13.14 | genes: 7 | CV error: 0.0857 (40 thresholds evaluated)For order thresholding the threshold is the number of largest |d_ik| values kept; for soft and hard thresholding it is a cut-off on |d_ik|.
compare_nsc() runs the evaluation protocol of the paper
(100 runs with different cross-validation folds in the paper; 2 here to
keep the vignette fast) and reports test error, balanced accuracy and
the number of genes.
res <- compare_nsc(dat$xtr, dat$ytr, dat$xte, dat$yte, n.runs = 2,
thresholding = c("hard", "order"), verbose = FALSE)
res
#> NSC comparison: 2 runs, 2 classes, 35 training samples
#> method error error.sd balanced.accuracy genes
#> HTh 0.200 0.042426 0.800 2.0
#> HTh3 0.105 0.035355 0.895 5.5
#> OTh 0.135 0.049497 0.865 6.5
#> OTh3 0.125 0.007071 0.875 2.5
recall_table(res)
#> class recall_HTh recall_HTh3 recall_OTh recall_OTh3
#> 1 A 0.62 0.86 0.82 0.78
#> 2 B 0.98 0.93 0.91 0.97The script
system.file("scripts", "simulation_correlated.R", package = "hetNSC")
contains the full correlated-gene simulation studies.