# \donttest{
## ------------------------------------------------------------
## RFQ, BRF, and standard classification
## ------------------------------------------------------------
data(breast, package = "randomForestSRC")
dta <- na.omit(breast)
f <- status ~ .
## Use 300 trees for illustration; imbalanced() defaults to 3000.
set.seed(19)
rfq <- imbalanced(f, dta, ntree = 300)
print(rfq)
print(get.imbalanced.performance(rfq, confusion = TRUE))
## Equivalent direct RFQ call using rfsrc().
## RFQ uses equal case weights and sampling without replacement.
set.seed(19)
rfq.direct <- rfsrc(f, dta, ntree = 300,
rfq = TRUE, splitrule = "auc", perf.type = "brier",
block.size = NULL, bootstrap = "by.root",
samptype = "swor", sampsize = function(n) n * .632,
case.wt = NULL)
print(all.equal(rfq$predicted.oob, rfq.direct$predicted.oob))
## Extract the class-frequency threshold from the observed responses.
print(get.rfq.threshold(rfq$yvar))
set.seed(20)
brf <- imbalanced(f, dta, ntree = 300, method = "brf")
## Equivalent direct BRF call using rfsrc().
## Each class has equal total sampling weight; sampling is with replacement.
class.frequency <- table(dta$status)
case.wt <- min(class.frequency) /
as.numeric(class.frequency[as.character(dta$status)])
set.seed(20)
brf.direct <- rfsrc(f, dta, ntree = 300,
rfq = FALSE, splitrule = "auc", perf.type = "brier",
block.size = NULL, bootstrap = "by.root",
samptype = "swr", sampsize = 2L * min(class.frequency),
case.wt = case.wt)
print(all.equal(brf$predicted.oob, brf.direct$predicted.oob))
standard <- imbalanced(f, dta, ntree = 300, method = "standard")
## Compare all three using the helper's default evaluation rule.
metrics <- c("sens", "spec", "gmean", "misclass", "brier.norm", "auc")
comparison <- rbind(
RFQ = as.numeric(get.imbalanced.performance(rfq)[metrics]),
BRF = as.numeric(get.imbalanced.performance(brf)[metrics]),
Standard = as.numeric(get.imbalanced.performance(standard)[metrics])
)
colnames(comparison) <- metrics
print(comparison)
## Supply a cutoff explicitly to inspect another probability rule.
print(get.imbalanced.performance(standard, threshold = .5))
## Let the saved classifier setting select prevalence or 0.5.
print(get.imbalanced.performance(rfq, threshold = "auto"))
print(get.imbalanced.performance(standard, threshold = "auto"))
print(get.imbalanced.performance(brf, threshold = "auto"))
## Request the 0.5 Bayes cutoff explicitly, including for RFQ.
print(get.imbalanced.performance(rfq, threshold = "bayes"))
## Vectors, one-column objects, and named two-column matrices.
class.frequency <- table(rfq$yvar)
minority <- names(class.frequency)[which.min(class.frequency)]
p <- rfq$predicted.oob[, minority]
print(get.imbalanced.performance(rfq$yvar, p))
print(get.imbalanced.performance(rfq$yvar, data.frame(p = p)))
print(get.imbalanced.performance(
rfq$yvar, rfq$predicted.oob[, rev(levels(rfq$yvar)), drop = FALSE]
))
## ------------------------------------------------------------
## Train/test assessment with a training-selected threshold
## ------------------------------------------------------------
set.seed(23)
trn.id <- sort(unlist(lapply(
split(seq_len(nrow(dta)), dta$status),
function(idx) idx[sample.int(length(idx),
size = floor(2 * length(idx) / 3))]
)))
trn <- dta[trn.id, , drop = FALSE]
tst <- dta[-trn.id, , drop = FALSE]
print(rbind(training = table(trn$status), test = table(tst$status)))
fit <- imbalanced(f, trn, ntree = 300)
selected <- get.imbalanced.optimize(fit, measure = "gmean",
ngrid = 501, plot.it = FALSE)
print(selected)
threshold.selected <- selected$threshold
pred <- predict(fit, newdata = tst)
print(get.imbalanced.performance(pred))
print(get.imbalanced.performance(pred, threshold = threshold.selected))
## Predicting labels for new rows does not require their outcomes.
class.frequency <- table(trn$status)
minority <- names(class.frequency)[which.min(class.frequency)]
majority <- setdiff(levels(trn$status), minority)
pred.x <- predict(fit, newdata = tst[, fit$xvar.names, drop = FALSE])
selected.class <- factor(
ifelse(pred.x$predicted[, minority] >= threshold.selected,
minority, majority),
levels = levels(trn$status)
)
print(head(selected.class))
## ------------------------------------------------------------
## Explicit permutation VIMP and blocking
## ------------------------------------------------------------
## RFQ uses normalized Brier performance by default.
rfq.vimp <- imbalanced(f, dta, ntree = 300,
importance = "permute", block.size = 1)
rfq.block <- imbalanced(f, dta, ntree = 300,
importance = "permute", block.size = 10)
print(cbind(unblocked = rfq.vimp$importance[, 1],
blocked = rfq.block$importance[, 1]))
## Small B is for illustration; increase it for final inference.
smp <- subsample(rfq.block, B = 25, verbose = FALSE)
print(extract.subsample(smp)$var.sel.Z)
plot.subsample(smp, cex.axis = .7, alpha = .5)
## ------------------------------------------------------------
## Strong imbalance: probability signals and permutation VIMP
## ------------------------------------------------------------
## x1 identifies a small subgroup; x2 changes risk within it.
## Both affect the outcome probability. The other variables are noise.
## Expected minority prevalence is .04 * (.15 + .5 / 2) = .016.
## Within the subgroup, all true probabilities exceed this prevalence,
## so the population RFQ class rule uses x1 alone. Brier performance
## can also detect the probability information supplied by x2.
local({
## One thread makes the blocked VIMP comparison reproducible.
old.options <- options(rf.cores = 1)
on.exit(options(old.options))
set.seed(47)
n <- 4000
p <- 20
x <- matrix(runif(n * p), ncol = p)
x[, 1] <- as.numeric(x[, 1] < .04)
colnames(x) <- paste0("x", seq_len(p))
prob <- x[, 1] * (.15 + .5 * x[, 2])
y <- factor(rbinom(n, 1, prob), levels = c(0, 1),
labels = c("majority", "minority"))
sim.data <- data.frame(y = y, x)
## Grow one RFQ forest, then use the same trees for both metrics.
## Brier is explicit here; it is also the RFQ default.
sim.fit <- imbalanced(y ~ ., sim.data, ntree = 500, mtry = p,
splitrule = "auc", perf.type = "brier",
block.size = 50, seed = -103)
brier <- vimp(sim.fit, importance = "permute", perf.type = "brier",
block.size = 50, seed = -107)
gmean <- vimp(sim.fit, importance = "permute", perf.type = "gmean",
block.size = 50, seed = -107)
counts <- table(y)
cat("\nClass counts and majority/minority ratio\n")
print(c(n.majority = unname(counts["majority"]),
n.minority = unname(counts["minority"]),
iratio = unname(counts["majority"] / counts["minority"])))
## Compare signals with noise within each metric's own scale.
vimp.values <- cbind(Brier = brier$importance[, 1],
Gmean = gmean$importance[, 1])
cat("\nPermutation VIMP: true signals and largest noise value\n")
print(rbind(vimp.values[1:2, , drop = FALSE],
largest.noise = apply(vimp.values[-(1:2), , drop = FALSE],
2, max)), digits = 5)
## A smaller rank is better; the ideal ranks are 1 and 2.
vimp.ranks <- apply(vimp.values, 2, function(z) {
rank(-z, ties.method = "average")
})
cat("\nRanks of the two true signals among all 20 variables\n")
print(vimp.ranks[1:2, , drop = FALSE])
})
# }
Run the code above in your browser using DataLab