---
title: "RFlocalfdr A Simulated Data Example \n"
author: "Rob Dunne"
date: "Wednesday, January 29, 2025"
output:
rmarkdown::html_vignette:
toc: true
toc_depth: 2
vignette: >
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{RFlocalfdr Simulated Data Example}
%\VignetteEncoding{UTF-8}
bibliography: paper.bib
header-includes:
\usepackage{dcolumn}
\usepackage{placeins}
\usepackage{float}
\usepackage{biblatex}
\usepackage{caption}
\usepackage{subcaption}
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(fig.width=6, fig.height=6, dpi=300,echo = FALSE)
knitr::opts_chunk$set(fig.pos = "H", out.extra = "")
```
\clearpage
# Introduction
Variable importance measures are often used to perform variable selection. It is possible to identify two final aims (@genuerVariableSelectionUsing2010):
1. find a small number of variables with a maximized accuracy, or
2. detect and rank all influential variables to focus on for interpretation and further exploration with domain experts.
It is the second objective that is more often our aim in biomolecular work. After the important variables are selected, they
can be joined with other data from the literature or databases such as gene annotations. They can also inform further rounds of experiments.
In additions, the reduced data set can be used in other analysis methods such as extraction of structural causal models.
The problem is characterized by:
1. a relatively small $n$ and a $p$ that is typically two to three orders of magnitude larger than $n$;
2. a lot of correlation between the variables;
# Discussion and Conclusion
We simulate a classification problem with 6350 variables of which 127 are related to the classification and a strong correlation structure.
We consider a number of variable ranking and selection methods. We give the results for variable selection, not for the performance of the classifier.
We consider two problems:
- ranking the variables by variable importance. We have taken the point on the ROC curve that
gives the minimum cost (with equal cost for positives and negatives). So these are the best results that could be
obtained for these methods.
- where there is a method of thresholding the importances to do variable selection, we have used it.
- where there is a method of threshold the variable importance to do variable selection, we consider the accuracy of the selected set of variables.
The results are given in the table below.
Considering the second part of the table where we select variables based on a RF impurity measure.
RFlocalfdr has done reasonably well on the FDR. Recursive Feature
Elimination has performed better but only selected a small number of
variables.
I could not get many of the SHAP implementations to converge in a 2 1/2 hour run time
limit. Shaff converged and returned one variable.
This is not surprising. The Shapley values are an "all subsets" calculation and this is NP hard. All of the methods
available use some heuristic to get around the hardness of the problem.
As this is quite interesting, I tried a number of SHAP implementations on the smaller that set that we used for visualizing the problem.
This goes beyond the original intent of a vignette for the **RFlocalfdr** package but seems interesting to pursue.
\marginpar{ select variables based on a RF impurity measure}
```{r table3, echo=FALSE, message=FALSE, warnings=FALSE, eval=FALSE}
tabl <- "
Method | true positives | false positives | Sensitivity | Specificity | FDR | AUC |
rank variables
RF impurity | 55 | 6 | 0.99 | 0.43 | 0.1 | 0.99 |
RF perm. | 61 | 3 | 0.99 | 0.48 | 0.05 | 0.74 |
RF corrected | 87 | 37 | 0.99 | 0.68 | 0.3 | 0.99 |
cforest | 61 | 3 | 0.99 | 0.84 | 0.12 | 0.99 |
SobolMDA | 61 | 3 | 0.99 | 0.48 | 0.05 | 0.99 |
Shaff | 1 | 0 | | | | |
iml Shapley | - | - | - | - | - | - |
fastshap | - | - | - | - | - | - |
select variables based on a RF impurity measure
RFlocalfdr | 59 | 36 | 0.99 | 0.46 | 0.34 | 0.73 |
Boruta | 2 | 2 | 0.99 | 0.016 | 0.5 | 0.51 |
AIR | 127 | 273 | 0.96 | 1 | 0.68 | 0.98 |
AIR+locfdr | 123 | 186 | 0.97 | 0.97 | 0.61 | 0.97 |
PIMP | 39 | 556 | 0.91 | 0.31 | 0.58 | 0.98 |
RFE | 22 | 2 | 0.99 | 0.17 | 0.08 | 0.58 |
"
kable(tabl, format = "markdown") #disasterous
```
\clearpage
```{R , echo=FALSE, eval=FALSE, result='asis',fig.cap='...'}
cat(tabl) # output the table in a format good for HTML/PDF/docx conversion
```
\clearpage
\begin{verbatim}
Method | true positives | false positives | Sensitivity | Specificity | FDR | AUC |
rank variables
RF impurity| 55 | 6 | 0.99 | 0.43 | 0.1 | 0.99 |
RF perm. | 61 | 3 | 0.99 | 0.48 | 0.045 | 0.74 |
RF corrected | 87 | 37 | 0.99 | 0.68 | 0.3 | 0.99 |
cforest | 61 | 3 | 0.99 | 0.84 | 0.12 | 0.99 |
SobolMDA | 61 | 3 | 0.99 | 0.48 | 0.05 | 0.99 |
Shaff | 1 | 0 | | | | |
fastshap | - | - | - | - | - | - |
iml Shapley | - | - | - | - | - | - |
select variables based on RF
RFlocalfdr | 59 | 36 | 0.99 | 0.46 | 0.34 | 0.73 |
Boruta | 2 | 2 | 0.99 | 0.016 | 0.5 | 0.51 |
AIR | 127 | 273 | 0.96 | 1 | 0.68 | 0.98 |
AIR+locfdr | 123 | 186 | 0.97 | 0.97 | 0.61 | 0.97 |
PIMP | 39 | 556 | 0.91 | 0.31 | 0.58 | 0.98 |
RFE | 22 | 2 | 0.99 | 0.17 | 0.08 | 0.58 |
The small data set (508 variables)
fastshap | 109 | 18 | 0.97 | 0.858 | 0.1 | - |
iml Shapley | - | - | - | - | - | - |
\end{verbatim}
\clearpage
# The Simulation
We simulate a small data set with a covariance structure to explore the actions of **RFlocalfdr**.
This package implements an empirical Bayes method of determining a significance level for the Random Forest
MDI (Mean Decrease in Impurity) variable importance measure.
See @dunneThresholdingGiniVariable2022 for details.
We
- simulate the data set
- fit a Random Forest model
- use RFlocalfdr to estimate the significant variables
- for comparison, we also estimate the significant variables using
- Boruta (@Kursa.and.Rudnicki.2010)
- Recursive Fearure Elimination
- AIR (@Nembrini.et.al.2018)
- PIMP (@Altmann.et.al.2010).
## data setup
This dataset consist of _bands_, with _blocks_ of $\{1, 2, 4, 8, 16, 32, 64\}$ of
identical variables. The variables are $\in \{0,1,2\}$, a common encoding for
genomic data where the numbers represent the number of copies of the minor allele. Only band 1 is used to calculate the
$y$ vector and $y$ is 1 if any of $X[, c(1, 2, 4, 8, 16, 32, 64)]$ is non-zero, so $y$ is unbalanced, containing more
1's than 0's.
We plot a small data set of $300\times 508$ to explore the data generation method, see figure \@ref(fig:simulation2).
For the problem we generate a data set with 50 bands and 300 observations so $X$ is $300 \times 6350$ with 127 non-null variables. We fit a
standard RF to this dataset and record the resulting MDI importance score.
```{r eval=FALSE,echo=TRUE}
library(ranger)
library(RFlocalfdr)
library(caret)
library(pROC)
```
```{r , eval=FALSE,echo=FALSE}
packageDescription("RFlocalfdr")$GithubSHA1
#source("/home/dun280/Dropbox/R_libraries/RF_localfdr/RFlocalfdr/R/my.pimp.s")
```
```{r eval=FALSE,echo=FALSE}
#just plot a small data set to show the structure
set.seed(13)
num_samples <- 300
num_bands <- 4
band_rank <- 6
num_vars <- num_bands * (2 ** (band_rank+1) -1)
print(num_vars)
X <- matrix(NA, num_samples , num_vars)
set.seed(123)
temp<-matrix(0,508,3)
var_index <- 1
for(band in 1:num_bands) {
for (rank in 0:band_rank) {
for (i in 1:2**rank) {
temp[var_index,]<-c(band , rank, var_index)
var_index <- var_index +1
}
}
}
#png("./supp_figures/small_simulated_data_set.png")
plot(temp[,1],ylim=c(0,10),type="p")
lines(temp[,2],type="b",col="red")
legend(0,10,c("band","rank"),pch=1,col=c(1,2))
#dev.off()
abline(v=c( 1, 2 , 4 , 8 ,16 ,32, 64 ))
table(temp[temp[,1]==1,2])
## # 0 1 2 3 4 5 6
## # 1 2 4 8 16 32 64
X <- matrix(NA, num_samples , num_vars)
set.seed(123)
var_index <- 1
for(band in 1:num_bands) {
for (rank in 0:band_rank) {
# cat("rank=",rank,"\n")
var <- sample(0:2, num_samples, replace=TRUE)
for (i in 1:2**rank) {
X[,var_index] <- var
var_index <- var_index +1
}
# print(X[1:2,1:140])
# system("sleep 3")
}
}
y <- as.numeric(X[,1] > 1 | X[,2] > 1 | X[,4] > 1 | X[,8] > 1 | X[,16] > 1 |
X[,32] > 1 | X[,64] > 1 )
data <- cbind(y,X)
colnames(data) <- c("y",make.names(1:num_vars))
rfModel <- ranger(data=data,dependent.variable.name="y", importance="impurity",
classification=TRUE, mtry=100,num.trees = 10000, replace=TRUE,
seed = 17 )
zz2 <-log(importance(rfModel))
plot(zz2,type="b",col=temp[,1]+1)
roccurve <- roc(c(rep(1,127),rep(0,508-127)),zz2)
plot(roccurve)
auc(roccurve) # 0.993
```
```{r simulation2, echo=FALSE, fig.cap="A small simulated data set. Each band contains blocks of size {1, 2, 4, 8, 16, 32, 64}, and each block consists of correlated (identical variables).", fig.align="center", out.width = '50%'}
knitr::include_graphics("./supp_figures/small_simulated_data_set.png")
```
```{r eval=FALSE,echo=TRUE}
# generate the data
set.seed(13)
num_samples <- 300
num_bands <- 50
band_rank <- 6
num_vars <- num_bands * (2 ** (band_rank+1) -1)
print(num_vars)
X <- matrix(NA, num_samples , num_vars)
set.seed(123)
var_index <- 1
for(band in 1:num_bands) {
for (rank in 0:band_rank) {
# cat("rank=",rank,"\n")
var <- sample(0:2, num_samples, replace=TRUE)
for (i in 1:2**rank) {
X[,var_index] <- var
var_index <- var_index +1
}
# print(X[1:2,1:140])
# system("sleep 3")
}
}
y <- as.numeric(X[,1] > 1 | X[,2] > 1 | X[,4] > 1 | X[,8] > 1 | X[,16] > 1 |
X[,32] > 1 | X[,64] > 1 )
```
## Random forest
### parameter setting for a RF
- **ranger** sets **mtry** to \(\lfloor \sqrt{p} \rfloor \)
- **randomForest** sets **mtry** \(\lfloor \sqrt{p} \rfloor\) for classification and \(\ \lfloor p/3 \rfloor \) for regression
- if there are only a few relevant variables out of many, which is the case in many genetic datasets, **mtry** should be set
high, so that the algorithm can find the relevant variables (@goldsteinRandomForestsGenetic2011). A large **mtry**
ensures that there is (with high probability) at least one strong variable in the set of **mtry** candidate variables.
- For high dimensional data they observe lower error rates for higher **mtry** values for both classification and
regression, corroborating (@goldsteinRandomForestsGenetic2011)
- **min.node.size** specifies the minimum number of observations in a terminal node. Setting it lower
leads to trees with a larger depth which means that more splits are performed until the terminal nodes. Both **ranger**
and **randomForest** set the default value to 1 for classification and 5 for regression.
- **num.trees** we have set this to a large value.
- **replace=FALSE**
We fit a Random Forest (@Breiman.2001) model and recover the Mean Decrease in Impurity variable importance.
```{r eval=FALSE,echo=TRUE}
data <- cbind(y,X)
colnames(data) <- c("y",make.names(1:num_vars))
rfModel <- ranger(data=data,dependent.variable.name="y", importance="impurity",
classification=TRUE,num.trees = 10000, replace=FALSE, mtry=100, seed = 17 )
imp <-log(ranger::importance(rfModel))
t2 <-count_variables(rfModel)
plot(density(imp))
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),imp)
plot(roccurve)
auc(roccurve) # 0.993
palette("default")
col<-c(1, rep(2,2), rep(3,4), rep(4, 8), rep(5,16 ), rep(6,32 ), rep(7,64) )
plot(1:1016,imp[1:1016],type="p",col=rep(col,8),pch=16,cex=0.8,
xlab="variable number",ylab="log importances")
plot(imp,type="p",col=rep(col,8),pch=16,cex=0.8,xlab="variable number",
ylab="log importances")
predictions <- imp
labels <- c(rep(1,127),rep(0,6350-127))
pred <- prediction(predictions, labels)
perf <- performance(pred, "tpr", "fpr")
plot(perf)
cost_perf = performance(pred, "cost")
pred@cutoffs[[1]][which.min(cost_perf@y.values[[1]])]
# X22
#-3.485894
predicted_values <-rep(0,6350);predicted_values[ which(imp> -3.485894) ]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,6350-127)))
conf_matrix
#predicted_values 0 1
# 0 6217 72
# 1 6 55
conf_matrix[2,1]/sum(conf_matrix[2,]) #[1]] 0.09836066 FDR
sensitivity(conf_matrix) #] 0.9990358 TP/(TP+FN)
specificity(conf_matrix) # 0.4330709 TN/(FP+TN)
```
We have set **num.trees** to a large value and **mtry** to a larger value than the default. This is because of the large number of variables.
\marginpar{refs for setting parameters for large p}
We plot the log importances for the first 8 bands (figure \@ref(fig:simulation2zz2)). The plot for all 50 bands it too compressed.
The blocks are shown in different colors. The bias described by @Strobl.et.al.2007 is clearly visible.
Blocks 1 to 7 of band 1 should be of equal expected influence on $y$, but the importance is decreasing as the number of
variables in the block is increased.
```{r simulation2zz2, echo=FALSE, fig.cap="The log importances for the first 8 bands.", fig.align="center", out.width = '50%'}
knitr::include_graphics("./supp_figures/simulation2_zz2.png")
```
## RFlocalfdr
### determine cutoff
```{r eval=FALSE,echo=TRUE}
cutoffs <- c(0,1,4,10)
#png("./supp_figures/simulated_data_determine_cutoff.png")
res.con<- determine_cutoff(imp,t2,cutoff=cutoffs,plot=c(0,1,4,10))
#dev.off()
```
We plot the kernel density estimate of the histogram of the data $y$, and the skew normal fit, $t_1$, using the data up to the
quantile $Q$, shown in red.
```{r , echo=FALSE, fig.cap="For this data set, the selected cutoff value is 0.", fig.align="center", out.width = '50%'}
knitr::include_graphics("./supp_figures/simulated_data_determine_cutoff.png")
```
```{r eval=FALSE,echo=TRUE}
plot(cutoffs,res.con[,3])
cutoffs[which.min(res.con[,3])]
```
By plotting $\max(|y - t_1|)$ versus the cutoff values, we determine the appropriate cutoff. In this case it is
just $t2>0$
### recover variables
```{r eval=FALSE,echo=TRUE}
temp<-imp[t2 > 0]
palette("R3")
qq <- plotQ(temp,debug.flag = 1)
ppp<-run.it.importances(qq,temp,debug=0)
aa<-significant.genes(ppp,temp,cutoff=0.05,do.plot=2)
length(aa$probabilities) # 95
tt1 <-as.numeric(gsub("X([0-9]*)","\\1",names(aa$probabilities)))
tt2 <- table(ifelse((tt1 < 127),1,2))
# 1 2
# 59 36
# The false discovery rate is 0.3789474
tt2[2]/(tt2[1]+tt2[2])
#59 36 36/(36+59) = 0.3789474
predicted_values<-rep(0, 6350);predicted_values[tt1]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,6350-127)))
conf_matrix
conf_matrix[2,1]/sum(conf_matrix[2,]) # 0.3789474 FDR
sensitivity(conf_matrix) #0.994215 TP/(TP+FN)
specificity(conf_matrix) #0.4645669 TN/(FP+TN)
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),predicted_values)
plot(roccurve)
auc(roccurve) #0.7294
accuracy<-(conf_matrix[1,1]+conf_matrix[2,2])/(sum(conf_matrix))
accuracy # 0.983622
```
The FDR is 0.379. We can also calculate the sensitivity, sensitivity and accruacy.
```{r eval=FALSE,echo=FALSE}
temp <- temp - min(temp) + .Machine$double.eps
palette(topo.colors(n = 7))
col<-c(1, rep(2,2), rep(3,4), rep(4, 8), rep(5,16 ), rep(6,32 ), rep(7,64) )
plot(1:127,temp[1:127],type="p",col=rep(col,2),pch=16,cex=0.8,
xlab="variable number",ylab="log importances")
plot(1:1016,temp[1:1016],type="p",col=rep(col,8),pch=16,cex=0.8,
xlab="variable number",ylab="log importances")
lines(1:1016,temp[1:1016],col = "gray62",lwd=0.5)
abline(h=3.699622,col="red")
abline(v=127,col="green")
legend("topright",legend=c("RFlocalfdr cutoff","non-null variables"),lty=1,col=c("red","green" ))
```
In order to make the comparisons with other methods, the following packages may need to be installed.
```{r eval=FALSE,echo=TRUE}
install.packages("Boruta")
install.packages("locfdr")
install.packages("vita")
install.packages("locfdr")
devtools::install_github("silkeszy/Pomona")
```
```{r eval=FALSE,echo=FALSE}
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("twilight")
```
## Boruta
We try the Boruta algorithm (@Kursa.and.Rudnicki.2010).
```{r eval=FALSE,echo=TRUE}
library(Boruta)
set.seed(120)
system.time(boruta1 <- Boruta(X,as.factor(y), num.threads = 7,getImp=getImpRfGini,
classification=TRUE,num.trees = 10000, replace=FALSE, mtry=100, seed = 17))
print(boruta1)
#Boruta performed 99 iterations in 19.54727 secs.
#4 attributes confirmed important: X4859, X58, X6132, X7;
# 6346 attributes confirmed unimportant: X1, X10, X100, X1000, X1001 and 6341 more;
plotImpHistory(boruta1)
aa <- which(boruta1$finalDecision=="Confirmed")
bb <- which(boruta1$finalDecision=="Tentative")
predicted_values <-rep(0,6350);predicted_values[c(aa,bb)]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,6350-127)))
conf_matrix
conf_matrix[2,1]/sum(conf_matrix[2,]) # 0.3789474 FDR
sensitivity(conf_matrix) #0.9996786 TP/(TP+FN)
specificity(conf_matrix) #0.01574803 TN/(FP+TN)
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),predicted_values)
plot(roccurve)
auc(roccurve) #0.5077
accuracy<-(conf_matrix[1,1]+conf_matrix[2,2])/(sum(conf_matrix))
accuracy #0.98
```
```{verbatim, eval=FALSE,echo=FALSE}
D !D
+ TP FP
- FN TN
- Sensitivity=[TP/(TP+FN)]×100
- Specificity=[TN/(FP+TN)]×100
- Positive predictive value(PPV)= [TP/(TP+FP)]×100
- Negative predictive value(NPV)=[TN/(FN+TN)]×100.
```
## Recursive Feature Elimination
- This is provided by the package Pomona (@Pomona.2022)
- recalculates importance scores in each step.
- does not recognize the **ranger** option **replace=FALSE**
\marginpar{why is ntree = 500}
```{r eval=FALSE,echo=TRUE}
library(Pomona)
colnames(X) <- make.names(1:dim(X)[2])
set.seed(111)
res <- var.sel.rfe(X, y, prop.rm = 0.2, recalculate = TRUE, tol = 10,
ntree = 500, mtry.prop = 0.2, nodesize.prop = 0.1, no.threads = 7,
method = "ranger", type = "classification", importance = "impurity",
case.weights = NULL)
res$var
#[1] "X1" "X106" "X11" "X12" "X127" "X13" "X15" "X16" "X2"
#[10] "X23" "X24" "X3" "X4" "X44" "X46" "X4885" "X5" "X54"
#[19] "X5474" "X6" "X7" "X72" "X9" "X91"
tt <-as.numeric(gsub("X([0-9]*)","\\1", res$var))
table(ifelse((tt < 127),1,2))
# 1 2
#21 3 0.0833
res<-c(1,106, 11, 12, 127, 13, 15, 16, 2, 23, 24, 3, 4, 44, 46, 4885, 5, 54, 5474, 6,
7, 72, 9, 91)
predicted_values <-rep(0,6350);predicted_values[c(res)]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,6350-127)))
conf_matrix
#predicted_values 0 1
# 0 6221 105
# 1 2 22
conf_matrix[2,1]/sum(conf_matrix[2,]) # 0.083 FDR
sensitivity(conf_matrix) # 0.9996786 TP/(TP+FN)
specificity(conf_matrix) #0 0.1732283 TN/(FP+TN)
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),predicted_values)
plot(roccurve)
auc(roccurve) #0.5865
accuracy<-(conf_matrix[1,1]+conf_matrix[2,2])/(sum(conf_matrix))
accuracy # 0.9831496
```
## AIR
See @Nembrini.et.al.2018.
- AIR is provided in the package _ranger_, using the option _impurity_corrected_.
- the suggested approach is to apply the **janitza* method and select the variables with $p < 0.05$
```{r eval=FALSE,echo=TRUE}
rfModel2 <- ranger(data=data,dependent.variable.name="y", importance="impurity_corrected",
classification=TRUE, mtry=100,num.trees = 10000, replace=TRUE, seed = 17 )
ww<- importance_pvalues( rfModel2, method = "janitza")
p <- ww[,"pvalue"]
cc <- which(p< 0.05)
predicted_values <-rep(0,6350);predicted_values[cc]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,6350-127)))
conf_matrix
#predicted_values 0 1
# 0 5950 0
# 1 273 127
#FDR is 273/(127+273) = 0.6825
sensitivity(conf_matrix) #0.9561305 TP/(TP+FN)
specificity(conf_matrix) #1 TN/(FP+TN)
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),predicted_values)
plot(roccurve)
auc(roccurve) # 0.9781
accuracy<-(conf_matrix[1,1]+conf_matrix[2,2])/(sum(conf_matrix))
accuracy # 0.9570079
```
As the null values, as modelled bt the AIR procedure, are symmetric around 0, the
question arises as to whether we can apply the locfdr (@locfdr.2015) procedure.
In this case, it reduces the FDR from 0.682 to 0.601.
```{r eval=FALSE,echo=TRUE}
plot(density(ww[,"importance"]))
imp <- ww[,"importance"]
#imp <-imp/sqrt(var(imp))
#plot(density(imp))
library(locfdr)
aa<-locfdr(imp,df=13)
which( (aa$fdr< 0.05) & (imp>0))
cc2 <- which( (aa$fdr< 0.05) & (imp>0))
length(cc2) # [1] 309
length(intersect(cc,cc2)) #[1] 309
(length(cc2) - length(which(cc2 <= 127)))/length(cc2) #[1] 0.6019417 fdr
predicted_values <-rep(0,6350);predicted_values[cc2]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,6350-127)))
conf_matrix
conf_matrix[2,1]/sum(conf_matrix[2,]) # 0.6019417 FDR
sensitivity(conf_matrix) # 0.9701109 TP/(TP+FN)
specificity(conf_matrix) #0.9685039 TN/(FP+TN)
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),predicted_values)
plot(roccurve)
auc(roccurve) # 0.9693
accuracy<-(conf_matrix[1,1]+conf_matrix[2,2])/(sum(conf_matrix))
accuracy # 0.9700787
```
## PIMP
note that PIMP
- permutes the response vector $y$ $S$ times to give $y^{*s}$
- For each permutation, a new RF is grown and the permutation
variable importance measure (VarImp$^{*s}$) for all predictor variables X is computed.
- The vector ‘perVarImp’ of $S$ VarImp measures for every predictor variable is used to approximate the null importance distributions
for each variable (see **PimpTest**)
- we are doing $p$ tests, hence a multiple testing correction is in order
- we base our work on the PIMP implementation provided by @vita.2015 _Variable Importance Testing Approaches_
- pump as described by Altman is only applicable to permutation importance.
- however we see no impediment in extending the method to Gini importance
- we have done that in our code and also
- fixed a small error
- extended the code to use **ranger** as well as **randomForest** implementations of RF
```{r eval=FALSE,echo=TRUE}
#vita: Variable Importance Testing Approaches
library(vita)
y<-factor(y)
X<-data.frame(X)
set.seed(117)
cl.ranger <- ranger(y~. , X,mtry = 100,num.trees = 10000, classification=TRUE, replace=FALSE,
importance = 'impurity')
system.time(pimp.varImp.cl<-my_ranger_PIMP(X,y,cl.ranger,S=10, parallel=TRUE, ncores=10))
pimp.t.cl <- vita::PimpTest(pimp.varImp.cl,para = FALSE)
aa <- summary(pimp.t.cl,pless = 0.1)
tt<-as.numeric(gsub("X([0-9]*)","\\1", names(which(aa$cmat2[,"p-value"]< 0.05))))
table(ifelse((tt < 127),1,2))
# 1 2
# 126 176 176 /( 126+ 176 ) = 0.582
predicted_values <-rep(0,6350);predicted_values[which(aa$cmat2[,"p-value"]< 0.05)]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,6350-127)))
conf_matrix
conf_matrix[2,1]/sum(conf_matrix[2,]) #[1] 0.5794 FDR
sensitivity(conf_matrix) #0.971 TP/(TP+FN)
specificity(conf_matrix) #1 TN/(FP+TN)
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),predicted_values)
plot(roccurve)
auc(roccurve) # 0.9859
accuracy<-(conf_matrix[1,1]+conf_matrix[2,2])/(sum(conf_matrix))
accuracy # 0.9724409
```
```{r eval=FALSE,echo=FALSE}
# twilight uses a stochastic downhill search algorithm for estimating the local false discovery rate
library(twilight)
p.values <- aa$cmat2[,"p-value"]
ans <- twilight(p.values) #Twilight cannot run properly.
fdr <- ans$result$fdr
sum(fdr < 0.05) #[1] 0
```
```{r eval=FALSE,echo=FALSE}
#how to make a neat html table in Rmarkdown
library(tables)
X <- rnorm(125, sd=100)
Group <- factor(sample(letters[1:5], 125, rep=TRUE))
tab <- tabular( Group ~ (N=1)+Format(digits=2)*X*((Mean=mean) + Heading("Std Dev")*sd) )
table_options(knit_print = FALSE)
tab # This chunk uses the default results = 'markup'
toHTML(tab) # This chunk uses results = 'asis'
table_options(CSS =
"<style>
#ID .center {
text-align:center;
background-color: aliceblue;
}
</style>", doCSS = TRUE)
tab
table_options(doCSS = FALSE)
```
\clearpage
# Other Methods
## permutation importance
```{R, echo=TRUE, eval=FALSE}
rfModel2 <- ranger(data=data,dependent.variable.name="y", importance="permutation",
classification=TRUE, mtry=100,num.trees = 10000, replace=FALSE, seed = 17 )
imp <-ranger::importance(rfModel2)
plot(density(imp))
palette(topo.colors(n = 7))
plot(1:1016,imp[1:1016],type="p",col=rep(col,2),pch=16,cex=0.8,
xlab="variable number",ylab="log importances")
lines(1:1016,imp[1:1016],col = "gray62",lwd=0.5)
abline(v=127,col="green")
# could apply Briemans and Cutlers argument that the permutation importance is Gaussian --
# or could use empirical Bayes
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),imp)
plot(roccurve)
auc(roccurve) # 0.9924
library(ROCR)
predictions <- imp
labels <- c(rep(1,127),rep(0,6350-127))
pred <- prediction(predictions, labels)
cost_perf = performance(pred, "cost")
pred@cutoffs[[1]][which.min(cost_perf@y.values[[1]])]
table(imp> 0.0001063063 ,c(rep(1,127),rep(0,6350-127)))
# 0 1
# FALSE 6220 66
# TRUE 3 61
predicted_values <-rep(0,6350);predicted_values[which(imp> 0.0001063063 )]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,6350-127)))
conf_matrix
#predicted_values 0 1
# 0 6220 66
# 1 3 61
conf_matrix[2,1]/sum(conf_matrix[2,]) #[1] 0.046875 FDR
sensitivity(conf_matrix) # 0.9995179 TP/(TP+FN)
specificity(conf_matrix) # 0.480315 TN/(FP+TN)
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),predicted_values)
plot(roccurve)
auc(roccurve) # 0.7399
accuracy<-(conf_matrix[1,1]+conf_matrix[2,2])/(sum(conf_matrix))
accuracy # 0.9724409
```
## impurity corrected
```{R, echo=TRUE, eval=FALSE}
rfModel3 <- ranger(data=data,dependent.variable.name="y", importance="impurity_corrected",
classification=TRUE, mtry=100,num.trees = 10000, replace=FALSE, seed = 17 )
imp <-ranger::importance(rfModel3)
plot(density(imp))
palette(topo.colors(n = 7))
plot(1:1016,imp[1:1016],type="p",col=col,pch=16,cex=0.8,
xlab="variable number",ylab="log importances")
lines(1:1016,imp[1:1016],col = "gray62",lwd=0.5)
abline(v=127,col="green")
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),imp)
plot(roccurve)
auc(roccurve) # 0.9951
library(ROCR)
predictions <- imp
labels <- c(rep(1,127),rep(0,6350-127))
pred <- prediction(predictions, labels)
cost_perf = performance(pred, "cost")
pred@cutoffs[[1]][which.min(cost_perf@y.values[[1]])]
table(imp> 0.0106588 ,c(rep(1,127),rep(0,6350-127)))
# 0 1
# FALSE 6186 40
# TRUE 37 87
#best FDR is
37/(37 + 87) #[1] ] 0.2983871
#
predicted_values <-rep(0,6350);predicted_values[ which(imp> 0.0106588)]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,6350-127)))
conf_matrix
conf_matrix[2,1]/sum(conf_matrix[2,]) #[1]] 0.2983871 FDR
sensitivity(conf_matrix) #0.9940543 TP/(TP+FN)
specificity(conf_matrix) #0.6850394 TN/(FP+TN)
```
## cforest
we consider the cforest forest fitting procedure (@hothorn_unbiased_2006) and the Conditional Permutation Importance (CPI)
(@Strobl.et.al.2008). It appears that the implementation of CPI provided by **party::varimp** is excessively slow.
We use the implementation provided by the **permimp** package.
```{R, echo=TRUE, eval=FALSE}
library(party)
library(permimp) #
data <- data.frame(y,X)
system.time(mod1.cf <- party::cforest(y ~ ., data = data,
control = party::cforest_unbiased(ntree = 10, mtry = 100)))
# 12.848 0.000 12.885
system.time(aa<-party::varimp(mod1.cf, conditional = TRUE))
# user system elapsed
#3089.301 9.484 3100.554
system.time(aa3<-permimp(mod1.cf, conditional = TRUE))
# user system elapsed
# 4.685 0.000 4.707
system.time(mod1.cf <- party::cforest(y ~ ., data = data,
control = party::cforest_unbiased(ntree = 100, mtry = 100)))
# 27.233 0.464 27.703
system.time(aa<-party::varimp(mod1.cf, conditional = TRUE))
# user system elapsed
#29380.980 81.511 29475.138
save(aa,file="aa.Rdata")
system.time(aa3<-permimp(mod1.cf, conditional = FALSE))
# user system elapsed
# 6.084 0.000 6.089
system.time(mod1.cf <- party::cforest(y ~ ., data = data,
ontrol = party::cforest_unbiased(ntree = 1000, mtry = 100)))
#35.438 0.636 36.095
system.time(aa4<-permimp(mod1.cf, conditional = FALSE))
# user system elapsed
# 47.135 0.429 47.578
palette("default")
col<-c(1, rep(2,2), rep(3,4), rep(4, 8), rep(5,16 ), rep(6,32 ), rep(7,64) )
plot(1:1016,aa4$values[1:1016],type="p",col=col,pch=16,cex=0.8,
xlab="variable number",ylab="log importances")
system.time(mod1.cf <- party::cforest(y ~ ., data = data,
control = party::cforest_unbiased(ntree = 10000, replace=FALSE, mtry = 100)))
#35.438 0.636 36.095
system.time(aa4<-permimp(mod1.cf, conditional = FALSE))
# user system elapsed
# 47.135 0.429 47.578
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),aa4$values)
plot(roccurve)
auc(roccurve) #0.9486
system.time(mod1.cf <- party::cforest(y ~ ., data = data,
control = party::cforest_unbiased(ntree = 10000, mtry = 100)))
#Process R killed at Mon Feb 26 07:18:23 2024 on robubuntu
# user system elapsed
# 98.673 2.832 101.977 on petra, seems to have used 26GB of RAM
system.time(aa4<-permimp(mod1.cf, conditional = TRUE))
# user system elapsed
#352.681 2.262 356.232
head(aa4$values)
# X1 X2 X3 X4 X5 X6
#2.600188e-05 9.454988e-05 1.217855e-04 1.682538e-04 1.572708e-04 1.706311e-04
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),aa4$values)
plot(roccurve)
auc(roccurve) # 0.9978
library(ROCR)
predictions <- aa4$values
labels <- c(rep(1,127),rep(0,6350-127))
pred <- prediction(predictions, labels)
perf <- performance(pred, "tpr", "fpr")
plot(perf)
cost_perf <- performance(pred, "cost")
pred@cutoffs[[1]][which.min(cost_perf@y.values[[1]])]
table(predictions> 3.264138e-05 ,c(rep(1,127),rep(0,6350-127)))
# 0 1
# FALSE 6220 66
# TRUE 3 61
predicted_values <-rep(0,6350);predicted_values[ which(predictions> 3.264138e-05 )]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,6350-127)))
conf_matrix
conf_matrix[2,1]/sum(conf_matrix[2,]) #[1]] 0.1229508 FDR
sensitivity(conf_matrix) # 0.9975896 TP/(TP+FN)
specificity(conf_matrix) # 0.8425197 TN/(FP+TN)
#https://stats.stackexchange.com/questions/61521/cut-off-point-in-a-roc-curve-is-there-a-simple-function
library(pROC)
rocobj <- roc(labels,predictions)
plot(rocobj)
coords(rocobj, "best")
coords(rocobj, x="best", input="threshold", best.method="youden") # Sa
auc(rocobj) # 0.9978
predicted_values <-rep(0,6350);predicted_values[ which(predictions> 8.976665e-06 )]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,6350-127)))
conf_matrix
#see also
#https://stats.stackexchange.com/questions/132547/roc-curves-for-unbalanced-datasets
#https://cran.r-project.org/web/packages/qeML/vignettes/Unbalanced_Classes.html
#https://juandelacalle.medium.com/how-and-why-i-switched-from-the-roc-curve-to-the-precision-recall-curve-to-analyze-my-imbalanced-6171da91c6b8
#https://machinelearningmastery.com/roc-curves-and-precision-recall-curves-for-imbalanced-classification/
```
## sobolMDA
@benardMDARandomForests2022
```{R, echo=FALSE, eval=FALSE}
#library(remotes)
#install_gitlab("drti/sobolmda")
library(sobolMDA)
data <- data.frame(y,X)
system.time(rng.sobolmda <- sobolMDA::ranger(dependent.variable.name = "y", classification = TRUE, data = data, replace=FALSE, mtry = 100,
num.trees = 10000, importance = "sobolMDA", seed = 125))
palette("default")
col<-c(1, rep(2,2), rep(3,4), rep(4, 8), rep(5,16 ), rep(6,32 ), rep(7,64) )
plot(1:1016,rng.sobolmda$variable.importance[1:1016],type="p",col=col,pch=16,cex=0.8,
xlab="variable number",ylab="log importances")
plot(rng.sobolmda$variable.importance,type="p",col=col,pch=16,cex=0.8,
xlab="variable number",ylab="log importances")
imp<- rng.sobolmda$variable.importance
roccurve <- roc(c(rep(1,127),rep(0,6350-127)),imp)
plot(roccurve)
auc(roccurve) #0.9954
library(ROCR)
predictions <- rng.sobolmda$variable.importance
labels <- c(rep(1,127),rep(0,6350-127))
pred <- prediction(predictions, labels)
perf <- performance(pred, "tpr", "fpr")
plot(perf)
#perf <- performance(pred, "acc")
#plot(perf, avg= "vertical", spread.estimate="boxplot", lwd=3,col='blue',
# show.spread.at= seq(0.1, 0.9, by=0.1),
# main= "Accuracy across the range of possible cutoffs")
#plot(0,0,type="n", xlim= c(0,0.001), ylim=c(0,10000),
# xlab="Cutoff", ylab="Density",
# main="How well do the predictions separate the classes?")
#
#lines(density(pred@predictions[[1]][pred@labels[[1]]=="1"]), col= "red")
#lines(density(pred@predictions[[1]][pred@labels[[1]]=="0"]), col="green")
cost_perf = performance(pred, "cost")
pred@cutoffs[[1]][which.min(cost_perf@y.values[[1]])]
table(rng.sobolmda$variable.importance> 0.001048671 ,c(rep(1,127),rep(0,6350-127)))
# 0 1
# FALSE 6220 66
# TRUE 3 61
predicted_values <-rep(0,6350);predicted_values[ which(imp> 0.001048671 )]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,6350-127)))
conf_matrix
#predicted_values 0 1
# 0 6220 66
# 1 3 61
conf_matrix[2,1]/sum(conf_matrix[2,]) #[1]] 0.046875 FDR
sensitivity(conf_matrix) # 0.9995179 TP/(TP+FN)
specificity(conf_matrix) # 0.480315 TN/(FP+TN)
#best FDR is
3/(61+3) #[1] 0.046875
```
```{R, echo=FALSE, eval=FALSE}
# example from Variable Importance in Random Forests: Traditional Methods and New Developments
#https://towardsdatascience.com/variable-importance-in-random-forests-20c6690e44e0
library(kernlab)
library(drf) #An implementation of distributional random forests as introduced in Cevid & Michel & Meinshausen & Buhlmann (2020)
library(Matrix)
library(DescTools)
library(mice)
library(sobolMDA)
source("./simulated2/compute_drf_vimp.R")
source("./simulated2/evaluation.R")
##Simulate Data that experiences both a mean as well as sd shift
n <- 200
# Simulate from X
x1 <- runif(n,-1,1)
x2 <- runif(n,-1,1)
X0 <- matrix(runif(7*n,-1,1), nrow=n, ncol=7)
x3 <- x1+ runif(n,-1,1)
X <- cbind(x1,x2, x3, X0)
# Simulate dependent variable Y
Y <- as.matrix(rnorm(n,mean = 0.8*(x1 > 0), sd = 1 + 1*(x2 > 0)))
colnames(X)<-paste0("X", 1:10)
head(cbind(Y,X))
Xfull <-X
## Variable importance for conditional Expectation Estimation
XY <- as.data.frame(cbind(Xfull, Y))
colnames(XY) <- c(paste('X', 1:(ncol(XY)-1), sep=''), 'Y')
num.trees <- 500
forest <- sobolMDA::ranger(Y ~., data = XY, num.trees = num.trees, importance = 'sobolMDA') #was this a classification model? check
sobolMDA <- forest$variable.importance
names(sobolMDA) <- colnames(X)
sort(sobolMDA, decreasing = T)
forest <- sobolMDA::ranger(Y ~., data = XY, num.trees = num.trees, importance = 'permutation')
MDA <- forest$variable.importance
names(MDA) <- colnames(X)
sort(MDA, decreasing = T)
MMDVimp <- compute_drf_vimp(X=X,Y=Y)
sort(MMDVimp, decreasing = T)
load("~/Downloads/wage_benchmark.Rdata")
##Define the training data
n<-1000
Xtrain<-X[1:n,]
Ytrain<-Y[1:n,]
Xtrain<-cbind(Xtrain,Ytrain[,"male"])
colnames(Xtrain)[ncol(Xtrain)]<-"male"
Ytrain<-Ytrain[,1, drop=F]
##Define the test data
ntest<-2000
Xtest<-X[(n+1):(n+ntest),]
Ytest<-Y[(n+1):(n+ntest),]
Xtest<-cbind(Xtest,Ytest[,"male"])
colnames(Xtest)[ncol(Xtest)]<-"male"
Ytest<-Ytest[,1, drop=F]
# Calculate variable importance for both measures
# 1. Sobol-MDA
XY <- as.data.frame(cbind(Xtrain, Ytrain))
colnames(XY) <- c(paste('X', 1:(ncol(XY)-1), sep=''), 'Y')
num.trees <- 500
forest <- sobolMDA::ranger(Y ~., data = XY, num.trees = num.trees, importance = 'sobolMDA')
SobolMDA <- forest$variable.importance
names(SobolMDA) <- colnames(Xtrain)
# 2. MMD-MDA
MMDVimp <- compute_drf_vimp(X=Xtrain,Y=Ytrain,silent=T)
print("Top 10 most important variables for conditional Expectation estimation")
sort(SobolMDA, decreasing = T)[1:10]
print("Top 5 most important variables for conditional Distribution estimation")
sort(MMDVimp, decreasing = T)[1:10]
# Remove variables one-by-one accoring to the importance values saved in SobolMDA
# and MMDVimp.
evallistSobol<-evalall(SobolMDA, X=Xtrain ,Y=Ytrain ,Xtest, Ytest, metrics=c("MSE"), num.trees )
evallistMMD<-evalall(MMDVimp, X=Xtrain ,Y=Ytrain ,Xtest, Ytest, metrics=c("MMD"), num.trees ) #slow
plot(1:79,evallistSobol$evalMSE, type="l", lwd=2, cex=0.8, col="darkgreen", main="MSE and (MMD+0.8) loss" , xlab="Number of Variables removed", ylab="Values")
lines(1:79,evallistMMD$evalMMD+0.8, type="l", lwd=2, cex=0.8, col="blue")
#Random Forests in 2023: Modern Extensions of a Powerful Method
#https://towardsdatascience.com/random-forests-in-2023-modern-extensions-of-a-powerful-method-b62debaf1d62
```
# Shap values
SHapley Additive exPlanations, (SHAP values), are used as a measure of variable importance. It
is based on Shapley values, which use game theory to assign credit for a model's prediction to each feature or feature
value.
The Shapley value is defined for any value
function and SHAP is a special case of the Shapley value by taking the function to be the conditional
expectation function of our model. A SHAP value is the average marginal contribution of a feature value across all possible sets of features.
- the very fast TreeSHAP is an algorithm to compute SHAP values for tree ensemble models such as decision trees, random forests, and
gradient boosted trees in a polynomial-time proposed by @lundbergConsistentIndividualizedFeature2019
- the R library *treeshap* only accepts regression models (not classification)
- other options
- *fastshap* uses Monte-Carlo sampling to approximate SHAP values,
- *shapr* and *kernelshap* provide implementations of KernelSHAP.
- shapper A port to Python’s “shap†package is provided in shapper
Global explanations
Aside from explaining individual prediction (i.e., local explanation), it can be useful to aggregate the results of
several (i.e., all of the training predictions) into an overall global summary about the model (i.e., global
explanations).
The interpretation of the Shapley value for feature value j is: The value of the j-th feature contributed $\phi_j$ to the prediction of this particular instance compared to the average prediction for the dataset.
## SHAFF: Fast and consistent SHApley eFfect estimates via random Forests
@benardSHAFFFastConsistent2022
1. improve the Monte-Carlo approach by using importance sampling to focus on the most relevant subsets of variables identified by the forest.
2. uses a projected random forest algorithm to compute fast and accurate estimates of the conditional expectations for any variable subset (same algorithm as SobolMDA?)
3. initial forest sets **mtry** to $p/3$, much larger than we are using
Finds one variable.
```{R, echo=FALSE, eval=FALSE}
# devtools::install_gitlab("drti/shaff")
library(shaff)
library("nnls") #needs this but does not load it
data <- cbind(y,X)
names(data)[1] <- "Y"
system.time(rng <- shaff(data = data, mtry = 100, num.trees = 10000, K=500))
which(rng>0)
#[1] 17
```
## iml: Interpretable Machine Learning
- See @Molnar.2021
- did not finish in 3 hours (20 processors but seems to be only using 1, can we make it parallel?)
```{R, echo=FALSE, eval=FALSE}
library(iml)
library(ranger)
set.seed(314)
data <- data.frame(y,X)
rfModel2 <- ranger(data=data,dependent.variable.name="y", importance="permutation",
classification=TRUE,mtry=50,num.trees = 1000, replace=TRUE, seed = 17 )
X <- data.frame(data[,-which(names(data) == "y")])
mod <- Predictor$new(rfModel2, data = X, type = "response")
shapley <- Shapley$new(mod, x.interest = X[1, ])
plot(shapley)
#full data set killed -- using a lot of RAM
# dim(X) [1] 300 508 OK
shapley$plot(sort = TRUE)
shapley$y.hat.interest
phis <- shapley$results$phi
names(phis) <- shapley$results$feature
par(mar=c(5, 10, 4, 2)+0.1) # for wide enough left margin
barplot(sort(phis), horiz=TRUE, las=1)
# ## pick features e.g. by a limit for the absolute value of phis
# values <- shapley$results
N <- 10 # number of features to show
# Capture the ggplot2 object
p <- plot(shapley, sort = TRUE)
# Modify it so it shows only top N features
print(p + scale_x_discrete(limits=rev(p$data$feature.value[order(-p$data$phi)][1:N])))
# dim(X) [1] 300 508 OK
library(iml)
library(ranger)
set.seed(314)
data <- data.frame(y,X)
rfModel2 <- ranger(data=data,dependent.variable.name="y", importance="permutation",
classification=TRUE,mtry=50,num.trees = 1000, replace=TRUE, seed = 17 )
X <- data.frame(data[,-which(names(data) == "y")])
mod <- Predictor$new(rfModel2, data = X, type = "response")
shapley <- Shapley$new(mod, x.interest = X[1, ])
plot(shapley)
```
## mlr3
## fastshap
- https://bgreenwell.github.io/fastshap/articles/fastshap.html
- did not finish in 2 1/2 hours with 20 processors
- tried the samll data set
-
```{R, echo=FALSE, eval=FALSE}
#Source: vignettes/fastshap.Rmd
library(fastshap)
library(ranger)
set.seed(2053) # for reproducibility
pfun <- function(object, newdata) {
predict(object, data = newdata)$predictions
}
library(doParallel)
# With parallelism
#salloc --account=OD-225217 --mem=100GB --nodes=1 --ntasks-per-node=1 --cpus-per-task=20 -J interactive -t 6:00:00 srun --pty /bin/bash -l
registerDoParallel(cores = 20) # use forking
set.seed(5038)
system.time({ # estimate run time
ex.ames.par <- explain(rfModel2, X = X, pred_wrapper = pfun, nsim = 50, adjust = TRUE, parallel = TRUE)
})
######## a test data set from the fastshap documentation
# Load the sample data; see ?datasets::mtcars for details
data(mtcars)
# Fit a projection pursuit regression model
fit <- ppr(mpg ~ ., data = mtcars, nterms = 5)
# Prediction wrapper
pfun <- function(object, newdata) { # needs to return a numeric vector
predict(object, newdata = newdata)
}
# Compute approximate Shapley values using 10 Monte Carlo simulations
set.seed(101) # for reproducibility
shap1 <- explain(fit, X = subset(mtcars, select = -mpg), nsim = 10, pred_wrapper = pfun)
head(shap1)
apply(shap1,2,f<-function(x){mean(abs(x))})
shap <- explain(fit, X = subset(mtcars, select = -mpg), nsim = 10, pred_wrapper = pfun,adjust=TRUE)
x <- iris[which(iris[,5] != "setosa"), c(1,5)]
trials <- 10000
ptime <- system.time({
r <- foreach(icount(trials), .combine=cbind) %dopar% {
ind <- sample(100, 100, replace=TRUE)
result1 <- glm(x[ind,2]~x[ind,1], family=binomial(logit))
coefficients(result1)
}
})[3]
ptime
############################################################################################
######## try the small data set
#salloc --account=OD-225217 --mem=100GB --nodes=1 --ntasks-per-node=1 --cpus-per-task=20 -J interactive -t 6:00:00 srun --pty /bin/bash -l
dim(data) #[1] 300 509
rfModel2 <- ranger(data=data,dependent.variable.name="y", importance="none",
classification=TRUE,mtry=50,num.trees = 1000, replace=TRUE, seed = 17 )
pfun <- function(object, newdata) { # needs to return a numeric vector
predict(object, data = newdata)$predictions
}
system.time(shap <- fastshap::explain(rfModel2, X = subset(data, select = -y), nsim = 10, pred_wrapper = pfun))
# user system elapsed
#201.161 31.122 153.928
head(shap)
plot(apply(shap,2,f<-function(x){mean(abs(x))}))
library(doParallel)
registerDoParallel(cores = 20) # use forking
set.seed(5038)
system.time({ # estimate run time
shap <- fastshap::explain(rfModel2, X = subset(data, select = -y), nsim = 1000, pred_wrapper = pfun, adjust = TRUE, parallel = TRUE)
})
# user system elapsed
#23285.927 5089.726 3168.814
predictions <- SHAP.global <- apply(shap,2,f<-function(x){mean(abs(x))})
palette(topo.colors(n = 7))
col<-c(1, rep(2,2), rep(3,4), rep(4, 8), rep(5,16 ), rep(6,32 ), rep(7,64) )
plot(1:508,predictions,type="p",col=col,pch=16,cex=0.8,
xlab="variable number",ylab="log importances")
lines(1:508,predictions,col = "gray62",lwd=0.5)
abline(v=127,col="green")
roccurve <- roc(c(rep(1,127),rep(0, 508-127)),predictions)
plot(roccurve)
auc(roccurve) # 0.9806
library(ROCR)
labels <- c(rep(1,127),rep(0,508-127))
pred <- prediction(predictions, labels)
cost_perf = performance(pred, "cost")
pred@cutoffs[[1]][which.min(cost_perf@y.values[[1]])]
#0.0002388418
table(SHAP.global> 0.0002388418 ,c(rep(1,127),rep(0,508-127)))
# 0 1
# FALSE 363 11
# TRUE 18 116
predicted_values <-rep(0,508);predicted_values[ which(predictions> 0.0002388418 )]<-1
conf_matrix<-table(predicted_values,c(rep(1,127),rep(0,508-127)))
conf_matrix
conf_matrix[2,1]/sum(conf_matrix[2,]) #[1]] 0.1343284 FDR
sensitivity(conf_matrix) # 0.9527559 TP/(TP+FN)
specificity(conf_matrix) # 0.9133858 TN/(FP+TN)
```