diff --git a/src/bayesvalidrox/post_processing/post_processing.py b/src/bayesvalidrox/post_processing/post_processing.py
index d67e3529fcc6eec98b63e390d905722faeb719e3..08bb5eb25cff244788359a3aeb6091617692e1d5 100644
--- a/src/bayesvalidrox/post_processing/post_processing.py
+++ b/src/bayesvalidrox/post_processing/post_processing.py
@@ -844,11 +844,11 @@ class PostProcessing:
         y_val = outputs
         if y_val is None:
             y_val, _ = self.engine.Model.run_model_parallel(samples, key_str="valid")
-        y_pce_val, _ = self.engine.eval_metamodel(samples=samples)
+        y_metamod_val, _ = self.engine.eval_metamodel(samples=samples)
 
         # Fit the data(train the model)
-        for key in y_pce_val.keys():
-            residuals = y_val[key] - y_pce_val[key]
+        for key in y_metamod_val.keys():
+            residuals = y_val[key] - y_metamod_val[key]
 
             # ------ Residuals vs. predicting variables ------
             # Check the assumptions of linearity and independence
@@ -874,7 +874,8 @@ class PostProcessing:
 
             # ------ Fitted vs. residuals ------
             # Check the assumptions of linearity and independence
-            plt.scatter(x=y_pce_val[key], y=residuals, color="blue", edgecolor="k")
+            for i in range(y_metamod_val[key].shape[0]):
+                plt.scatter(x=y_metamod_val[key][i,:], y=residuals[i,:], color="blue", edgecolor="k")
             plt.title(f"{key}: Residuals vs. fitted values")
             plt.grid(True)
             plt.hlines(