SHAP (SHapley Additive exPlanations)及DALEX预测单样本变量情况和shapviz可视化学习

rm(list = ls())
load("TCGA_HNSC_practice.Rdata")

library(DALEX)
library(shapviz)

# prepare model
colnames(meta)
head(meta)[1:3,1:8]
#                                ID OS  OS.time gender  race site neoadjuvant age
# TCGA-CR-7374-01A TCGA-CR-7374-01A  0 1.000000 FEMALE WHITE   CR          No  67
# TCGA-CV-A45V-01A TCGA-CV-A45V-01A  1 1.066667 FEMALE WHITE   CV          No  87
# TCGA-CV-7102-01A TCGA-CV-7102-01A  1 1.866667 FEMALE WHITE   CV          No  76
meta$OS <- factor(meta$OS)

# 提取数据，去除NA
# 要注意在建立模型的时候变量可以使非数字
# shap解释变量的时候一定要把变量内容设定为数字
dat <- na.omit(meta[,c(2,8,11,12,17,41:43)])
head(dat)
#                  OS age day_completion month_completion lymph_count T N M
# TCGA-CV-7102-01A  1  76              6               10          24 3 1 0
# TCGA-MT-A67D-01A  0  55             21                5          23 2 0 0
# TCGA-P3-A6T4-01A  1  54             21               10          14 4 1 0
# TCGA-CV-6934-01A  1  66             13                9          19 3 2 0
# TCGA-D6-6824-01A  0  61              7                9          25 4 1 0
# TCGA-CN-A642-01A  1  57             17                6         110 3 2 0
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) # 设置输出概率
model
# Ranger result
# 
# Call:
#  ranger(OS ~ age + day_completion + month_completion + lymph_count +      T + N + M, data = dat, classification = TRUE, probability = TRUE) 
# 
# Type:                             Probability estimation 
# Number of trees:                  500 
# Sample size:                      369 
# Number of independent variables:  7 
# Mtry:                             2 
# Target node size:                 10 
# Variable importance mode:         none 
# Splitrule:                        gini 
# OOB prediction error (Brier s.):  0.2189129 

explainer <- explain(model, 
                     data = dat[,-1], 
                     y = dat$OS==1,
                     label = "Random forest",
                     colorize = FALSE)
# Preparation of a new explainer is initiated
#   -> model label       :  Random forest 
#   -> data              :  369  rows  7  cols 
#   -> target variable   :  369  values 
#   -> predict function  :  yhat.ranger  will be used (  default  )
#   -> predicted values  :  No value for predict function target column. (  default  )
#   -> model_info        :  package ranger , ver. 0.16.0 , task classification (  default  ) 
#   -> model_info        :  Model info detected classification task but 'y' is a logical . Converted to numeric.  (  NOTE  )
#   -> predicted values  :  numerical, min =  0.0617535 , mean =  0.4353077 , max =  0.8829342  
#   -> residual function :  difference between y and yhat (  default  )
#   -> residuals         :  numerical, min =  -0.5614393 , mean =  -0.001703359 , max =  0.6008169  
#   A new explainer has been created!   

# 创建一个空列表来存储每个样本的 SHAP 值
shap_values_list <- list()

# 循环遍历每一行数据并计算 SHAP 值
for (i in 1:nrow(dat)) {
  # 提取当前样本的特征值
  new_observation <- dat[i, -1] # 去掉目标列
  names(new_observation) <- colnames(dat)[-1] # 确保列名一致
  
  # 计算 SHAP 值并存储
  shap_values_list[[i]] <- predict_parts(
    explainer = explainer,
    new_observation = new_observation,
    type = "shap",
    B = 25 # 选择多少个排列组合
  )
  
  # 打印当前进度
  cat("Processed observation:", i, "\n")
}

# 查看 SHAP 值列表的结构
str(shap_values_list)

shap_total <- model_parts(
    explainer = explainer,
    new_observation = dat[,-1],
    type = "variable_importance", # difference,ratio, raw
    B = 25 # 选择多少个排列组合
  )
shap_total

plot(shap_values_list[1])

# 初始化一个空列表，用于存储 shapviz 结果
shapviz_res_list <- list()
# 遍历 shap_values_list
for (i in seq_along(shap_values_list)) {
  # 对第 i 个 SHAP 值计算 shapviz
  shapviz_res_list[[i]] <- shapviz(shap_values_list[[i]])
  # 打印进度
  cat("Processed shapviz for sample:", i, "\n")
}
# 检查生成的 shapviz 对象列表
print(shapviz_res_list)

sv_waterfall(shapviz_res_list[[2]])
sv_force(shapviz_res_list[[2]])
sv_importance(shapviz_res_list[[2]])