diff --git a/src/bayesvalidrox/post_processing/post_processing.py b/src/bayesvalidrox/post_processing/post_processing.py
index 2d3b16c46939be0e276361adefbd58b0e610edde..4ee885b3b1eed290e30ace9563bbdb73be0b356e 100644
--- a/src/bayesvalidrox/post_processing/post_processing.py
+++ b/src/bayesvalidrox/post_processing/post_processing.py
@@ -871,32 +871,34 @@ class PostProcessing:
 
             y_pce_val_ = y_pce_val[key]
             y_val_ = y_val[key]
+            residuals = y_val_ - y_pce_val_
 
             # ------ Residuals vs. predicting variables ------
             # Check the assumptions of linearity and independence
             fig1 = plt.figure()
-            plt.title(key+": Residuals vs. Predicting variables")
-            residuals = y_val_ - y_pce_val_
-            plt.scatter(x=y_val_, y=residuals, color='blue', edgecolor='k')
-            plt.grid(True)
-            xmin, xmax = min(y_val_), max(y_val_)
-            plt.hlines(y=0, xmin=xmin*0.9, xmax=xmax*1.1, color='red', lw=3,
-                       linestyle='--')
-            plt.xlabel(key)
-            plt.ylabel('Residuals')
-            plt.show()
+            for i, par in enumerate(MetaModel.ExpDesign.par_names):
+                plt.title(f"{key}: Residuals vs. {par}")
+                plt.scatter(
+                    x=samples[:, i], y=residuals, color='blue', edgecolor='k')
+                plt.grid(True)
+                xmin, xmax = min(samples[:, i]), max(samples[:, i])
+                plt.hlines(y=0, xmin=xmin*0.9, xmax=xmax*1.1, color='red',
+                           lw=3, linestyle='--')
+                plt.xlabel(par)
+                plt.ylabel('Residuals')
+                plt.show()
 
-            if save_fig:
-                # save the current figure
-                fig1.savefig(f'./{newpath}/Residuals_vs_PredVariables.pdf',
-                             bbox_inches='tight')
-                # Destroy the current plot
-                plt.clf()
+                if save_fig:
+                    # save the current figure
+                    fig1.savefig(f'./{newpath}/Residuals_vs_Par_{i+1}.pdf',
+                                 bbox_inches='tight')
+                    # Destroy the current plot
+                    plt.clf()
 
             # ------ Fitted vs. residuals ------
             # Check the assumptions of linearity and independence
             fig2 = plt.figure()
-            plt.title(key+": Residuals vs. predicting variables")
+            plt.title(f"{key}: Residuals vs. fitted values")
             plt.scatter(x=y_pce_val_, y=residuals, color='blue', edgecolor='k')
             plt.grid(True)
             xmin, xmax = min(y_val_), max(y_val_)
@@ -949,7 +951,7 @@ class PostProcessing:
             plt.yticks()
             plt.xlabel("Theoretical quantiles")
             plt.ylabel("Sample quantiles")
-            plt.title(key+": Q-Q plot of normalized residuals")
+            plt.title(f"{key}: Q-Q plot of normalized residuals")
             plt.grid(True)
             plt.show()