I am doing a project research predicting top_20 label for cyclists and races dataset, that is 1 when a cyclist has arrived in the top 20 position, and 0 otherwise. The dataset is heavily unbalanced (92000 instances of class 1 and about 400000 of class 0). I run a basic Neural Network with a simple architecture. I am also using class weights. My NN is not performing well but this is not the point, since I know, from project purpose, that the data is not so good. The f1 score is good for class 0 (majority) and is about 0.87. For class 1 (minority) is about 0.55. This is fine, but the problem arises when trying to explain the results with SHAP.

Using SHAP I have a base value, for all of the instances, of 1 (or 0.98 somethimes). The summary plot seems to be reasonable, but the force and waterfall plot are not. How the base values can lead torwards 1 if the model is overconfident in predicting class 0 instead ? I also thought that the label is inverted in shap values and actually the expected value of 1 reflect class 0. But I am confused on how the structure of my shap value object is, since I cannot select one class to analyze, but it is divided by instances. So for example shap_value[1] refers to instance number 1, and not to class 1. And so on. The code is this:

def build_model(optimizer='adam', dropout_rate=0.5, num_units_1=128, num_units_2=64):
model = Sequential([
    Dense(num_units_1, activation='relu', input_shape=(X_train.shape[1],)),  # Input layer
    BatchNormalization(),
    Dropout(dropout_rate),  # Dropout layer
    Dense(num_units_2, activation='relu', kernel_regularizer=l2(0.01)),  # Hidden layer
    #Dropout(dropout_rate),  # Dropout layer
    Dense(1, activation='sigmoid')  # Output layer for binary classification
])

# For convergence
lr_schedule = ExponentialDecay(
    initial_learning_rate=0.001,
    decay_steps=10000,
    decay_rate=0.9
)

# Compile the model
optimizer_instance = {
    'adam': Adam(learning_rate=lr_schedule),
    'rmsprop': RMSprop(learning_rate=lr_schedule),
    'sgd': SGD(learning_rate=lr_schedule)
}[optimizer]

model.compile(optimizer=optimizer_instance,
              loss='binary_crossentropy',
              metrics=['accuracy'])
return model

For SHAP:

background = X_train_df.sample(200) 
test_sample = X_test_df.sample(100)   

# Explainer
explainer = shap.Explainer(best_model.predict, background)
shap_values = explainer(test_sample)

# Visualize SHAP values for each prediction
shap.summary_plot(shap_values, test_sample)

The plot highliths this that seems reasonable to me: lower delta values can reflect the fact that the cyclist is near the first position (delta is the difference in time w.r.t the first position) and so it pushes to 1

But for example the waterfall/force plots seems totally inconsistent to me:

And finally the structure of my shap values is this:

shap_values[2]
.values =
array([-5.55111512e-17,  1.00000000e-02,  1.73472348e-17,  3.12250226e-17,
    2.42861287e-17, -3.46944695e-18,  1.04083409e-17,  0.00000000e+00])

.base_values =
 0.99

.data =
array([1.20000000e+02, 1.57000000e+05, 1.55100000e+03, 3.28000000e+02,
   1.71790754e-02, 3.83373426e-06, 2.40317094e+01, 1.54471545e-01])