# for the current variable of interest, xname,
# create the list of variables to condition on:
create_cond_list <- function(cond, threshold, xname, input) {
stopifnot(is.logical(cond))
if (!cond) return(NULL)
if (threshold > 0 & threshold < 1) {
ctrl <- ctree_control(teststat = "quad", testtype = "Univariate", stump = TRUE)
xnames <- names(input)
xnames <- xnames[xnames != xname]
ct <- ctree(as.formula(paste(xname, "~", paste(xnames, collapse = "+"), collapse = "")),
data = input, controls = ctrl)
crit <- ct@tree$criterion[[2]]
crit[which(is.na(crit))] <- 0
return(xnames[crit > threshold])
}
stop()
}
## mincriterion = 0 so that complete tree is evaluated;
## regulate size of considered tree here via, e.g., mincriterion = 0.95
## or when building the forest in the first place via cforest_control(mincriterion = 0.95)
varimp <- function (object, mincriterion = 0, conditional = FALSE,
threshold = 0.2, nperm = 1, OOB = TRUE, pre1.0_0 = conditional)
{
response <- object@responses
if (length(response@variables) == 1 &&
inherits(response@variables[[1]], "Surv"))
return(varimpsurv(object, mincriterion, conditional, threshold, nperm, OOB, pre1.0_0))
input <- object@data@get("input")
xnames <- colnames(input)
inp <- initVariableFrame(input, trafo = NULL)
y <- object@responses@variables[[1]]
if(length(response@variables) != 1)
stop("cannot compute variable importance measure for multivariate response")
if (conditional || pre1.0_0) {
if(!all(complete.cases(inp@variables)))
stop("cannot compute variable importance measure with missing values")
}
CLASS <- all(response@is_nominal)
ORDERED <- all(response@is_ordinal)
if (CLASS) {
error <- function(x, oob) mean((levels(y)[sapply(x, which.max)] !=
y)[oob])
}
else {
if (ORDERED) {
error <- function(x, oob) mean((sapply(x, which.max) !=
y)[oob])
}
else {
error <- function(x, oob) mean((unlist(x) - y)[oob]^2)
}
}
w <- object@initweights
if (max(abs(w - 1)) > sqrt(.Machine$double.eps))
warning(sQuote("varimp"), " with non-unity weights might give misleading results")
## list for several permutations
perror <- matrix(0, nrow = nperm*length(object@ensemble), ncol = length(xnames))
## this matrix is initialized with values 0 so that a tree that does not
## contain the current variable adds importance 0 to its average importance
colnames(perror) <- xnames
for (b in 1:length(object@ensemble)){
tree <- object@expand(object@ensemble[[b]])
## if OOB == TRUE use only oob observations, otherwise use all observations in learning sample
if(OOB){oob <- object@weights[[b]] == 0} else{ oob <- rep(TRUE, length(y))}
p <- .Call("R_predict", tree, inp, mincriterion, -1L, PACKAGE = "atlantisPartyMod")
eoob <- error(p, oob)
## for all variables (j = 1 ... number of variables)
for(j in unique(varIDs(tree))){
for (per in 1:nperm){
if (conditional || pre1.0_0) {
tmp <- inp
ccl <- create_cond_list(conditional, threshold, xnames[j], input)
if (is.null(ccl)) {
perm <- sample(which(oob))
} else {
perm <- conditional_perm(ccl, xnames, input, tree, oob)
}
tmp@variables[[j]][which(oob)] <- tmp@variables[[j]][perm]
p <- .Call("R_predict", tree, tmp, mincriterion, -1L,
PACKAGE = "atlantisPartyMod")
} else {
p <- .Call("R_predict", tree, inp, mincriterion, as.integer(j),
PACKAGE = "atlantisPartyMod")
}
## run through all rows of perror
perror[(per+(b-1)*nperm), j] <- (error(p, oob) - eoob)
} ## end of for (per in 1:nperm)
} ## end of for(j in unique(varIDs(tree)))
} ## end of for (b in 1:length(object@ensemble))
perror <- as.data.frame(perror)
#return(MeanDecreaseAccuracy = perror) ## return the whole matrix (= nperm*ntree values per variable)
return(MeanDecreaseAccuracy = colMeans(perror)) ## return only averages over permutations and trees
}
varimpsurv <- function (object, mincriterion = 0, conditional = FALSE,
threshold = 0.2, nperm = 1, OOB = TRUE, pre1.0_0 = conditional)
{
cat("\n")
cat("Variable importance for survival forests; this feature is _experimental_\n\n")
response <- object@responses
input <- object@data@get("input")
xnames <- colnames(input)
inp <- initVariableFrame(input, trafo = NULL)
y <- object@responses@variables[[1]]
weights <- object@initweights
stopifnot(inherits(y, "Surv"))
if (conditional || pre1.0_0) {
if(!all(complete.cases(inp@variables)))
stop("cannot compute variable importance measure with missing values")
}
stopifnot(require("ipred", quietly = TRUE))
error <- function(x, oob) sbrier(y[oob,,drop = FALSE], x[oob])
pred <- function(tree, newinp, j = -1L) {
where <- R_get_nodeID(tree, inp, mincriterion)
wh <- .Call("R_get_nodeID", tree, newinp, mincriterion, as.integer(j), PACKAGE = "atlantisPartyMod")
swh <- sort(unique(wh))
RET <- vector(mode = "list", length = length(wh))
for (i in 1:length(swh)) {
w <- weights * (where == swh[i])
RET[wh == swh[i]] <- list(mysurvfit(y, weights = w))
}
return(RET)
}
w <- object@initweights
if (max(abs(w - 1)) > sqrt(.Machine$double.eps))
warning(sQuote("varimp"), " with non-unity weights might give misleading results")
## list for several permutations
perror <- matrix(0, nrow = nperm*length(object@ensemble), ncol = length(xnames))
## this matrix is initialized with values 0 so that a tree that does not
## contain the current variable adds importance 0 to its average importance
colnames(perror) <- xnames
for (b in 1:length(object@ensemble)){
tree <- object@ensemble[[b]]
## if OOB == TRUE use only oob observations, otherwise use all observations in learning sample
if(OOB){oob <- object@weights[[b]] == 0} else{ oob <- rep(TRUE, length(y))}
p <- pred(tree, inp)
eoob <- error(p, oob)
## for all variables (j = 1 ... number of variables)
for(j in unique(varIDs(tree))){
for (per in 1:nperm){
if (conditional || pre1.0_0) {
tmp <- inp
ccl <- create_cond_list(conditional, threshold, xnames[j], input)
if (is.null(ccl)) {
perm <- sample(which(oob))
} else {
perm <- conditional_perm(ccl, xnames, input, tree, oob)
}
tmp@variables[[j]][which(oob)] <- tmp@variables[[j]][perm]
p <- pred(tree, tmp, -1L)
} else {
p <- pred(tree, inp, as.integer(j))
}
## run through all rows of perror
perror[(per+(b-1)*nperm), j] <- (error(p, oob) - eoob)
} ## end of for (per in 1:nperm)
} ## end of for(j in unique(varIDs(tree)))
} ## end of for (b in 1:length(object@ensemble))
perror <- as.data.frame(perror)
#return(MeanDecreaseAccuracy = perror) ## return the whole matrix (= nperm*ntree values per variable)
return(MeanDecreaseAccuracy = colMeans(perror)) ## return only averages over permutations and trees
}
# cutpoints_list() returns:
# - vector of cutpoints (length=number of cutpoints)
# if variable is continuous
# - vector of indicators (length=number of categories x number of cutpoints)
# if variable is categorical (nominal or ordered)
cutpoints_list <- function(tree, variableID) {
cutp <- function(node) {
if (node[[4]]) return(NULL)
cp <- NULL
if (node[[5]][[1]] == variableID)
cp <- node[[5]][[3]]
nl <- cutp(node[[8]])
nr <- cutp(node[[9]])
return(c(cp, nl, nr))
}
return(cutp(tree))
}
conditional_perm <- function(cond, xnames, input, tree, oob){
## get cutpoints of all conditioning variables of the current variable of interest
## and generate design matrix for permutation from factors in help
blocks <- vector(mode = "list", length = length(cond))
for (i in 1:length(cond)) {
## varID is variable index or column number of input (predictor matrix)
## not variable name!
varID <- which(xnames == cond[i])
## if conditioning variable is not used for splitting in current tree
## proceed with next conditioning variable
cl <- cutpoints_list(tree, varID)
if (is.null(cl)) next
## proceed cutpoints for different types of variables
x <- input[, varID]
xclass <- class(x)[1]
if (xclass == "integer") xclass <- "numeric"
block <- switch(xclass, "numeric" = cut(x, breaks = c(-Inf, sort(unique(cl)), Inf)),
"ordered" = cut(as.numeric(x), breaks = c(-Inf, sort(unique(cl)), Inf)),
"factor" = {
CL <- matrix(as.logical(cl), nrow = nlevels(x))
rs <- rowSums(CL)
dlev <- (1:nrow(CL))[rs %in% rs[duplicated(rs)]]
fuse <- c()
for (ii in dlev) {
for (j in dlev[dlev > ii]) {
if (all(CL[ii,] == CL[j,])) fuse <- rbind(fuse, c(ii, j))
}
}
xlev <- 1:nlevels(x)
newl <- nlevels(x) + 1
block <- as.integer(x)
for (l in xlev) {
if (NROW(fuse) == 0) break
if (any(fuse[, 1] == l)) {
f <- c(l, fuse[fuse[, 1] == l, 2])
fuse <- fuse[!fuse[,1] %in% f, , drop = FALSE]
block[block %in% f] <- newl
newl <- newl + 1
}
}
as.factor(block)
})
blocks[[i]] <- block
}
## remove non-splitting variables
names(blocks) <- cond
blocks <- blocks[!sapply(blocks, is.null)]
## if none of the conditioning variables are used in the tree
if (!length(blocks)>0){
perm <- sample(which(oob))
return(perm)
} else {
blocks <- as.data.frame(blocks)
## from factors blocks create design matrix
f <- paste("~ - 1 + ", paste(colnames(blocks), collapse = ":", sep = ""))
des <- model.matrix(as.formula(f), data = blocks)
## one conditional permutation
perm <- 1:nrow(input)
for (l in 1:ncol(des)) {
index <- which(des[,l] > 0 & oob)
if (length(index) > 1)
perm[index] <- sample(index)
}
return(perm[oob])
}
}
varimpAUC <- function(object, mincriterion = 0, conditional = FALSE,
threshold = 0.2, nperm = 1, OOB = TRUE, pre1.0_0 = conditional)
{
response <- object@responses
input <- object@data@get("input")
xnames <- colnames(input)
inp <- initVariableFrame(input, trafo = NULL)
y <- object@responses@variables[[1]]
if(length(response@variables) != 1)
stop("cannot compute variable importance measure for multivariate response")
if (conditional || pre1.0_0) {
if(!all(complete.cases(inp@variables)))
stop("cannot compute variable importance measure with missing values")
}
CLASS <- all(response@is_nominal)
ORDERED <- all(response@is_ordinal)
if (CLASS) {
if (nlevels(y)>2) {
warning("AUC=TRUE works only for binary y\n error rate is used instead of AUC")
error <- function(x, oob) mean((levels(y)[sapply(x, which.max)] != y)[oob])
}
else {
error <- function(x, oob) {
xoob <- sapply(x, function(x) x[1])[oob]
yoob <- y[oob]
which1 <- which(yoob==levels(y)[1])
noob1 <- length(which1)
noob <- length(yoob)
if (noob1==0|noob1==noob) { return(NA) } # AUC cannot be computed if all OOB-observations are from one class
return(1-sum(kronecker(xoob[which1] , xoob[-which1],">"))/(noob1*(length(yoob)-noob1))) # calculate AUC
}
}
### stop
}
else {
if (ORDERED) {
error <- function(x, oob) mean((sapply(x, which.max) !=
y)[oob])
}
else {
error <- function(x, oob) mean((unlist(x) - y)[oob]^2)
}
}
w <- object@initweights
if (max(abs(w - 1)) > sqrt(.Machine$double.eps))
warning(sQuote("varimp"), " with non-unity weights might give misleading results")
perror <- matrix(0, nrow = nperm*length(object@ensemble), ncol = length(xnames))
colnames(perror) <- xnames
for (b in 1:length(object@ensemble)){
tree <- object@ensemble[[b]]
if(OOB){oob <- object@weights[[b]] == 0} else{ oob <- rep(TRUE, length(xnames))}
p <- .Call("R_predict", tree, inp, mincriterion, -1L, PACKAGE = "atlantisPartyMod")
eoob <- error(p, oob)
for(j in unique(varIDs(tree))){
for (per in 1:nperm){
if (conditional || pre1.0_0) {
tmp <- inp
ccl <- create_cond_list(conditional, threshold, xnames[j], input)
if (is.null(ccl)) {
perm <- sample(which(oob))
} else {
perm <- conditional_perm(ccl, xnames, input, tree, oob)
}
tmp@variables[[j]][which(oob)] <- tmp@variables[[j]][perm]
p <- .Call("R_predict", tree, tmp, mincriterion, -1L,
PACKAGE = "atlantisPartyMod")
} else {
p <- .Call("R_predict", tree, inp, mincriterion, as.integer(j),
PACKAGE = "atlantisPartyMod")
}
perror[(per+(b-1)*nperm), j] <- (error(p, oob) - eoob)
}
}
}
perror <- as.data.frame(perror)
return(MeanDecreaseAccuracy = colMeans(perror, na.rm = TRUE)) ## na.rm = TRUE because with AUC-perm. VIM NA values occur whenever a tree's OOB-observations are all from the same class
}
Datasets
Models
