Example usage
## Loading required package: MASS
## Loading required package: foreach
## Loading required package: doSNOW
## Loading required package: iterators
## Loading required package: snow
##
## Attaching package: 'misc'
## The following object is masked from 'package:stats':
##
## integrate
# Create test functions
library(MASS)
# 1. Nonlinear function with interactions
f_polynomial <- function(X) {
X <- as.matrix(X)
X[,1]^2 + 2*X[,2] + 0.5*X[,1]*X[,2] + 0.1*X[,3]^3
}
# 2. Neural network-like function
f_neural <- function(X) {
X <- as.matrix(X)
tanh(X[,1] + 2*X[,2]) + sin(X[,3]) + X[,1]*X[,2]
}
# 3. Exponential function
f_exponential <- function(X) {
X <- as.matrix(X)
exp(0.1*X[,1]) + log(abs(X[,2]) + 1) + sqrt(abs(X[,3]) + 1)
}
# Generate test data
set.seed(42)
n_train <- 100
n_test <- 5
X_train <- matrix(rnorm(n_train * 3), nrow = n_train, ncol = 3)
colnames(X_train) <- c("x1", "x2", "x3")
X_test <- matrix(rnorm(n_test * 3), nrow = n_test, ncol = 3)
colnames(X_test) <- c("x1", "x2", "x3")
cat("=== Testing Gradient-Based SHAP Methods ===\n\n")
## === Testing Gradient-Based SHAP Methods ===
# Test 1: Polynomial function
cat("1. Polynomial Function: x1² + 2x2 + 0.5x1*x2 + 0.1x3³\n")
## 1. Polynomial Function: x1² + 2x2 + 0.5x1*x2 + 0.1x3³
result_grad <- gradient_shap(f_polynomial, X_train, X_test, h = 1e-4)
cat("Gradient SHAP - Max residual:", round(result_grad$max_residual, 4), "\n")
## Gradient SHAP - Max residual: 0.6513
result_ig <- integrated_gradients(f_polynomial, X_train, X_test, n_steps = 20, h = 1e-4)
cat("Integrated Gradients - Max residual:", round(result_ig$max_residual, 4), "\n")
## Integrated Gradients - Max residual: 0.0028
result_lime <- local_linear_shap(f_polynomial, X_train, X_test, n_samples = 500)
cat("Local Linear - Max residual:", round(result_lime$max_residual, 4), "\n")
## Local Linear - Max residual: 1.0293
cat("Predictions:", round(result_grad$predictions, 4), "\n")
## Predictions: -0.1157 1.6824 1.8961 -2.4356 -0.0229
cat("Gradient SHAP attributions:\n")
## Gradient SHAP attributions:
## x1 x2 x3
## [1,] 0.0014 0.0591 -0.0001
## [2,] 1.2820 1.3563 -0.1293
## [3,] 0.0037 2.1819 -0.3287
## [4,] 0.5950 -2.7438 -0.0740
## [5,] 0.0554 0.1040 0.0617
cat("Integrated Gradients attributions:\n")
## Integrated Gradients attributions:
## x1 x2 x3
## [1,] 0.0003 0.0594 0.0000
## [2,] 0.6488 1.2526 -0.0448
## [3,] 0.0019 2.1802 -0.1135
## [4,] 0.3050 -2.5402 -0.0257
## [5,] 0.0258 0.1065 0.0208
# Test 2: Neural network-like function
cat("2. Neural-like Function: tanh(x1 + 2x2) + sin(x3) + x1*x2\n")
## 2. Neural-like Function: tanh(x1 + 2x2) + sin(x3) + x1*x2
result_neural <- gradient_shap(f_neural, X_train, X_test)
cat("Max residual:", round(result_neural$max_residual, 4), "\n")
## Max residual: 0.6285
cat("RMSE residual:", round(result_neural$rmse_residual, 4), "\n")
## RMSE residual: 0.3191
cat("Predictions:", round(result_neural$predictions, 4), "\n")
## Predictions: -0.1905 0.617 0.1449 -2.4482 0.3483
## Attributions:
## x1 x2 x3
## [1,] -0.0345 0.0582 -0.0604
## [2,] 0.4378 0.5692 -0.5431
## [3,] 0.0069 0.1876 -0.5233
## [4,] -0.7969 -1.1119 -0.5010
## [5,] -0.1651 0.0953 0.4973
# Test 3: Exponential function
cat("3. Exponential Function: exp(0.1x1) + log(|x2|+1) + sqrt(|x3|+1)\n")
## 3. Exponential Function: exp(0.1x1) + log(|x2|+1) + sqrt(|x3|+1)
result_exp <- gradient_shap(f_exponential, X_train, X_test, baseline_method = "median")
cat("Max residual:", round(result_exp$max_residual, 4), "\n")
## Max residual: 0.2436
cat("Predictions:", round(result_exp$predictions, 4), "\n")
## Predictions: 2.0907 2.7989 3.1198 3.1626 2.2781
## Attributions:
## x1 x2 x3
## [1,] -0.0094 -0.0108 0.0223
## [2,] 0.0723 0.3721 0.2768
## [3,] -0.0051 0.5330 0.3539
## [4,] 0.0694 0.5240 0.2373
## [5,] -0.0233 -0.0346 0.2427
# Comparison of methods for polynomial function
cat("=== Method Comparison (Polynomial Function) ===\n")
## === Method Comparison (Polynomial Function) ===
cat("True function: x1² + 2x2 + 0.5x1*x2 + 0.1x3³\n")
## True function: x1² + 2x2 + 0.5x1*x2 + 0.1x3³
cat("Test point:", round(X_test[1,], 4), "\n")
## Test point: -0.0046 -0.0579 -0.071
cat("True prediction:", round(f_polynomial(X_test[1,,drop=FALSE]), 4), "\n\n")
## True prediction: -0.1157
cat("Gradient SHAP attributions:", round(result_grad$attributions[1,], 4), "\n")
## Gradient SHAP attributions: 0.0014 0.0591 -1e-04
cat("Integrated Gradients:", round(result_ig$attributions[1,], 4), "\n")
## Integrated Gradients: 3e-04 0.0594 0
cat("Local Linear:", round(result_lime$attributions[1,], 4), "\n")
## Local Linear: 0.0019 0.0589 -3e-04
cat("\nResiduals (smaller = better):\n")
##
## Residuals (smaller = better):
cat("Gradient SHAP:", round(result_grad$residuals[1], 4), "\n")
## Gradient SHAP: -8e-04
cat("Integrated Gradients:", round(result_ig$residuals[1], 4), "\n")
## Integrated Gradients: 0
cat("Local Linear:", round(result_lime$residuals[1], 4), "\n")
## Local Linear: -0.3858
# EXPLICIT DECOMPOSITION VERIFICATION
cat("\n=== EXPLICIT DECOMPOSITION CHECK ===\n")
##
## === EXPLICIT DECOMPOSITION CHECK ===
cat("For first test point:\n")
## For first test point:
# Check each method
methods <- list(
"Gradient SHAP" = result_grad,
"Integrated Gradients" = result_ig,
"Local Linear" = result_lime
)
for (method_name in names(methods)) {
result <- methods[[method_name]]
cat("\n", method_name, ":\n", sep="")
cat(" Prediction:", round(result$predictions[1], 6), "\n")
cat(" Base value:", round(result$baseline_prediction, 6), "\n")
cat(" Sum of attributions:", round(sum(result$attributions[1,]), 6), "\n")
cat(" Base + Sum:", round(result$baseline_prediction + sum(result$attributions[1,]), 6), "\n")
cat(" Residual:", round(result$residuals[1], 6), "\n")
cat(" Decomposition valid?", abs(result$residuals[1]) < 1e-5, "\n") # More reasonable threshold
}
##
## Gradient SHAP:
## Prediction: -0.115655
## Base value: -0.175333
## Sum of attributions: 0.060451
## Base + Sum: -0.114882
## Residual: -0.000774
## Decomposition valid? FALSE
##
## Integrated Gradients:
## Prediction: -0.115655
## Base value: -0.175333
## Sum of attributions: 0.059677
## Base + Sum: -0.115656
## Residual: 1e-06
## Decomposition valid? TRUE
##
## Local Linear:
## Prediction: -0.115655
## Base value: 0.209629
## Sum of attributions: 0.060487
## Base + Sum: 0.270116
## Residual: -0.385771
## Decomposition valid? FALSE
# Show the mathematical relationship more clearly
cat("\n=== MATHEMATICAL VERIFICATION ===\n")
##
## === MATHEMATICAL VERIFICATION ===
pred <- result_ig$predictions[1]
base <- result_ig$baseline_prediction
attrs <- result_ig$attributions[1,]
cat("Mathematical check for Integrated Gradients:\n")
## Mathematical check for Integrated Gradients:
cat("prediction =", round(pred, 6), "\n")
## prediction = -0.115655
cat("base_value + sum(attributions) =", round(base, 6), "+", round(sum(attrs), 6), "=", round(base + sum(attrs), 6), "\n")
## base_value + sum(attributions) = -0.175333 + 0.059677 = -0.115656
## Difference = 5.855e-07
cat("Passes additivity test:", abs(pred - (base + sum(attrs))) < 1e-5, "\n")
## Passes additivity test: TRUE