SHAP (SHapley Additive exPlanations)及kernelshap预测单样本/全局情况和shapviz可视化学习

rm(list = ls())
library(kernelshap)
library(ggplot2)
library(ranger)
library(shapviz)

load("TCGA_HNSC_practice.Rdata")

colnames(meta)
head(meta)[1:3,1:8]
meta$OS <- factor(meta$OS)

# 提取数据，去除NA
# 要注意在建立模型的时候变量可以使非数字
# shap解释变量的时候一定要把变量内容设定为数字
dat <- na.omit(meta[,c(2,8,11,12,17,41:43)])
head(dat)
dat[,-1] <- data.frame(lapply(dat[,-1], as.numeric))
dat$age <- ifelse(dat$age > 60, 1,0)
dat$day_completion <- ifelse(dat$day_completion > 15, 1,0)
dat$month_completion <- ifelse(dat$month_completion > 6, 1,0)
dat$lymph_count <- ifelse(dat$lymph_count > 20, 1,0)
dat$`T` <- ifelse(dat$`T` > 2, 1,0)
dat$`N` <- ifelse(dat$`N` > 1, 1,0)
dat$`M` <- ifelse(dat$`M` > 0, 1,0)
str(dat)

# 使用randomForest/ranger去构建随机森林树模型
library(randomForest)
model <- randomForest(OS ~.,data =dat)

library(ranger)
model <- ranger(OS ~  age + day_completion + month_completion + lymph_count + `T` + `N` +`M`, 
                data = dat,classification = TRUE,
                probability = TRUE,
                num.trees = 100,
                seed = 20) # 设置输出概率
model

xvars <- c("age", "day_completion", "month_completion", 
          "lymph_count","T","N","M")

# 1) Sample rows to be explained
# 数据量不大,笔者使用全部数据
set.seed(10)
# X <- dat[sample(nrow(diamonds), 1000), xvars] # 官方代码，建议1000个样本
X <- dat[, xvars]

# 2) Optional: Select background data. If unspecified, 200 rows from X are used
# bg_X <- dat[sample(nrow(dat), 200), ]

# 3) Crunch SHAP values 
# Note: Since the number of features is small, we use permshap()
ps <- permshap(model, X) #, bg_X = bg_X)
ps
# SHAP values of first observations:
# $`0`
#             age day_completion month_completion lymph_count           T
# [1,] 0.01530288    -0.02829701      -0.02270137 0.015690372 -0.01515734
# [2,] 0.04099188     0.07402535       0.03189634 0.009485131  0.11013120
#                N           M
# [1,]  0.03507547 0.001203453
# [2,] -0.00196143 0.001408363
# 
# $`1`
#              age day_completion month_completion  lymph_count           T
# [1,] -0.01530288     0.02829701       0.02270137 -0.015690372  0.01515734
# [2,] -0.04099188    -0.07402535      -0.03189634 -0.009485131 -0.11013120
#                N            M
# [1,] -0.03507547 -0.001203453
# [2,]  0.00196143 -0.001408363


# Kernel SHAP gives almost the same:
ks <- kernelshap(model, X, bg_X = bg_X)
ks

# 4) Analyze with {shapviz}
ps <- shapviz(ps)

# 重要性分析
sv_importance(ps,kind = "bar") # "bar", "beeswarm", "both", "no"
sv_importance(ps,kind = "beeswarm")+  theme_bw() 

# 依赖图
sv_dependence(ps, "day_completion")

# 瀑布图
sv_waterfall(ps, row_id = 1) +
  theme(axis.text = element_text(size = 11))

# 力图
sv_force(ps, row_id = 1)

# 选择特定特征属性
sv_waterfall(ps, ps$`1`$X$`T`=="1") +
  theme(axis.text = element_text(size = 11))