######################################################################################
#                                                                                    #
#                       R-code for reproducing Figure 2                              #
#                                                                                    #
######################################################################################

require(randomForest)

##################################################################
# Function for creating data and computing VIs                   #
##################################################################

get_VIs <- function(nnoise = 50, nobs1 = 10, nobs2 = 10,
                    effect1 = 0.5, effect2 = 0.75, effect3 = 1){

    set.seed(128083)
    
    # Create predictors (5 of each effect strength and nnoise noise pred.)
    
    noise <- matrix(rnorm(nnoise * (nobs1 + nobs2)), ncol = nnoise)
    strong.effects  <- rbind(matrix(rnorm(nobs1*5, mean = effect3), ncol = 5),
                              matrix(rnorm(nobs2*5, mean = 0      ), ncol = 5))
    moderate.effects <- rbind(matrix(rnorm(nobs1*5, mean = effect2), ncol = 5),
                              matrix(rnorm(nobs2*5, mean = 0      ), ncol = 5))
    weak.effects     <- rbind(matrix(rnorm(nobs1*5, mean = effect1), ncol = 5),
                              matrix(rnorm(nobs2*5, mean = 0      ), ncol = 5))

    data <- data.frame(y = as.factor(c(rep(1, nobs1), rep(0, nobs2))),
                           cbind(strong.effects, moderate.effects, weak.effects,
                                 noise))
                                 
    # Grow random forest
    
    forest <- randomForest(y ~ ., data = data, 
                           importance  = TRUE,
                           replace     = FALSE, 
                           ntree       = 1000)
                           
    # Return importance matrix containing class specific VIs
    
    return(imp = importance(forest, scale = FALSE))
    
}
                                                       
##################################################################
# Call function get_VIs()                                        #
##################################################################
                                   
# Get VI matrix for one simulated dataset with a 5% imbalance level

VI_matrix <- get_VIs(nobs1 = 25, nobs2 = 475)  

# Get class specific VIs and original VI

VIs_from_minority_class <- VI_matrix[,"1"]
VIs_from_majority_class <- VI_matrix[,"0"]
VIs_original            <- VI_matrix[,"MeanDecreaseAccuracy"]


##################################################################
# Plot class specific and original VIs (Figure 2)                #
##################################################################

graphics.off()

pdf(file = "AUC_VIM/Class_specific_VIs/Figure2.pdf", width = 25, height = 16)
par(mfrow = c(3,1), cex.main = 4.5, cex.axis = 3, oma = c(0,2,0,0), mar = c(0,0,0,0), mai = c(1.05,1,1,0))

barplot(VIs_from_minority_class, xlab = "", ylab = "", ylim = c(0, 0.06), main = "VIs using only OOB observations from minority class", las = 2)
barplot(VIs_from_majority_class, xlab = "", ylab = "", ylim = c(0, 0.06), main = "VIs using only OOB observations from majority class", las = 2)
barplot(VIs_original,            xlab = "", ylab = "", ylim = c(0, 0.06), main = "VIs using all OOB observations", las = 2)

graphics.off()

