Changes in Scree_plot and scatterplot

fsoares232 · fsoares232 · commit 5be5852eb5ae · 2023-11-23T13:26:17.000Z
Changes in scree_plot: divided into two functions and added more metrics to the plot
and changes in scatterplot to alter the marker size
diff --git a/.ipynb_checkpoints/Multivariate Analysis - Supervised Analysis with PLS-DA-checkpoint.ipynb b/.ipynb_checkpoints/Multivariate Analysis - Supervised Analysis with PLS-DA-checkpoint.ipynb
@@ -223,7 +223,7 @@
    "outputs": [],
    "source": [
     "# Plot the spectra in the dataset\n",
-    "_scatterplots(np.log(np.mean(x_data, axis=0) + 1), xaxis=retention_times, yaxis=mz_values)"
+    "_scatterplots(np.log(np.mean(x_data, axis=0) + 1), xaxis=retention_times, yaxis=mz_values, marker_size=3)"
    ]
   },
   {
@@ -420,16 +420,16 @@
    "source": [
     "# The remaining plots with the new features I added into the pyChemometrics - flsoares232\n",
     "# Plot the weights\n",
-    "pls_da.plot_model_parameters(parameter='w', component=1, instrument = 'lcms', xaxis=retention_times, yaxis=mz_values)\n",
+    "pls_da.plot_model_parameters(parameter='w', component=1, instrument = 'lcms', xaxis=retention_times, yaxis=mz_values, marker_size=5)\n",
     "\n",
     "# Plot the loadings\n",
-    "pls_da.plot_model_parameters(parameter='p', component=1, instrument = 'lcms', xaxis=retention_times, yaxis=mz_values)\n",
+    "pls_da.plot_model_parameters(parameter='p', component=1, instrument = 'lcms', xaxis=retention_times, yaxis=mz_values, marker_size=5)\n",
     "\n",
     "# Plot beta\n",
-    "pls_da.plot_model_parameters(parameter='beta', instrument = 'lcms', xaxis=retention_times, yaxis=mz_values)\n",
+    "pls_da.plot_model_parameters(parameter='beta', instrument = 'lcms', xaxis=retention_times, yaxis=mz_values, marker_size=5)\n",
     "\n",
     "# Plot VIP\n",
-    "pls_da.plot_model_parameters(parameter='VIP', instrument = 'lcms', xaxis=retention_times, yaxis=mz_values)"
+    "pls_da.plot_model_parameters(parameter='VIP', instrument = 'lcms', xaxis=retention_times, yaxis=mz_values, marker_size=5)"
    ]
   },
   {
@@ -488,6 +488,15 @@
     "We apply a similar criterion as the one used with PCA: choosing as the number of components after which the $Q^{2}Y$ value reaches a plateau (less than 5% increase compared to previous number of components). "
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pls_da.scree_cv(x_train_log, y_train, total_comps=10)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -499,7 +508,16 @@
    },
    "outputs": [],
    "source": [
-    "pls_da.scree_plot(x_train_log, y_train, total_comps=10)"
+    "pls_da.scree_plot(metric = ['Q2Y','R2Y','Accuracy']) \n",
+    "# Possible metrics (1 to 3 options):\n",
+    "# 'Q2Y'\n",
+    "# 'R2Y'\n",
+    "# 'AUC'\n",
+    "# 'MCC'\n",
+    "# 'Recall'\n",
+    "# 'Precision'\n",
+    "# 'F1'\n",
+    "# 'Accuracy'"
    ]
   },
   {
diff --git a/Multivariate Analysis - Supervised Analysis with PLS-DA.ipynb b/Multivariate Analysis - Supervised Analysis with PLS-DA.ipynb
@@ -223,7 +223,7 @@
    "outputs": [],
    "source": [
     "# Plot the spectra in the dataset\n",
-    "_scatterplots(np.log(np.mean(x_data, axis=0) + 1), xaxis=retention_times, yaxis=mz_values)"
+    "_scatterplots(np.log(np.mean(x_data, axis=0) + 1), xaxis=retention_times, yaxis=mz_values, marker_size=3)"
    ]
   },
   {
@@ -420,16 +420,16 @@
    "source": [
     "# The remaining plots with the new features I added into the pyChemometrics - flsoares232\n",
     "# Plot the weights\n",
-    "pls_da.plot_model_parameters(parameter='w', component=1, instrument = 'lcms', xaxis=retention_times, yaxis=mz_values)\n",
+    "pls_da.plot_model_parameters(parameter='w', component=1, instrument = 'lcms', xaxis=retention_times, yaxis=mz_values, marker_size=5)\n",
     "\n",
     "# Plot the loadings\n",
-    "pls_da.plot_model_parameters(parameter='p', component=1, instrument = 'lcms', xaxis=retention_times, yaxis=mz_values)\n",
+    "pls_da.plot_model_parameters(parameter='p', component=1, instrument = 'lcms', xaxis=retention_times, yaxis=mz_values, marker_size=5)\n",
     "\n",
     "# Plot beta\n",
-    "pls_da.plot_model_parameters(parameter='beta', instrument = 'lcms', xaxis=retention_times, yaxis=mz_values)\n",
+    "pls_da.plot_model_parameters(parameter='beta', instrument = 'lcms', xaxis=retention_times, yaxis=mz_values, marker_size=5)\n",
     "\n",
     "# Plot VIP\n",
-    "pls_da.plot_model_parameters(parameter='VIP', instrument = 'lcms', xaxis=retention_times, yaxis=mz_values)"
+    "pls_da.plot_model_parameters(parameter='VIP', instrument = 'lcms', xaxis=retention_times, yaxis=mz_values, marker_size=5)"
    ]
   },
   {
@@ -488,6 +488,15 @@
     "We apply a similar criterion as the one used with PCA: choosing as the number of components after which the $Q^{2}Y$ value reaches a plateau (less than 5% increase compared to previous number of components). "
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pls_da.scree_cv(x_train_log, y_train, total_comps=10)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -499,7 +508,16 @@
    },
    "outputs": [],
    "source": [
-    "pls_da.scree_plot(x_train_log, y_train, total_comps=10)"
+    "pls_da.scree_plot(metric = ['Q2Y','R2Y','Accuracy']) \n",
+    "# Possible metrics (1 to 3 options):\n",
+    "# 'Q2Y'\n",
+    "# 'R2Y'\n",
+    "# 'AUC'\n",
+    "# 'MCC'\n",
+    "# 'Recall'\n",
+    "# 'Precision'\n",
+    "# 'F1'\n",
+    "# 'Accuracy'"
    ]
   },
   {
@@ -1118,8 +1136,8 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
-   "language": "python",
+   "display_name": "Python 3 (Spyder)",
+   "language": "python3",
    "name": "python3"
   },
   "language_info": {
diff --git a/pyChemometrics/ChemometricsOrthogonalPLS.py b/pyChemometrics/ChemometricsOrthogonalPLS.py
@@ -670,7 +670,7 @@ def plot_scores(self, orthogonal_component=1, color=None, discrete=False, label_
         return None
 
     ####### flsoares232 version - Updated 20-10-2023
-    def plot_model_parameters(self, parameter='w_pred', orthogonal_component=1, cross_val=False, sigma=2, bar=False, xaxis=None, yaxis=None, instrument=None):
+    def plot_model_parameters(self, parameter='w_pred', orthogonal_component=1, cross_val=False, sigma=2, bar=False, xaxis=None, yaxis=None, instrument=None, marker_size=3):
         
         """
         Plot different model parameters related with the variables
@@ -685,23 +685,23 @@ def plot_model_parameters(self, parameter='w_pred', orthogonal_component=1, cros
         :return:
         """
     
-        choices = {'w_pred': self.w_pred, 'p_pred': self.p_pred, 'w_ortho': self.w_ortho, 'p_ortho': self.p_ortho}
+        choices = {'w_pred': self.w_pred, 'p_pred': self.p_pred, 'w_ortho': self.w_ortho, 'p_ortho': self.p_ortho, 'beta': self.beta_coeffs, 'VIP': self.VIP()}
         choices_cv = {'wpred': 'Wpred_w_pred', }
     
         # decrement component to adjust for python indexing
         orthogonal_component -= 1
         
         # Beta and VIP don't depend on components so have an exception status here
         if cross_val is True:
-            if parameter in ['w_pred', 'p_pred']:
+            if parameter in ['w_pred', 'p_pred', 'beta', 'VIP']:
                 mean = self.cvParameters['Mean_' + choices_cv[parameter]].squeeze()
                 error = sigma * self.cvParameters['Stdev_' + choices_cv[parameter]].squeeze()
             else:
                 mean = self.cvParameters['Mean_' + choices_cv[parameter]][:, orthogonal_component]
                 error = sigma * self.cvParameters['Stdev_' + choices_cv[parameter]][:, orthogonal_component]
         else:
             error = None
-            if parameter in ['w_pred', 'p_pred']:
+            if parameter in ['w_pred', 'p_pred', 'beta', 'VIP']:
                 mean = choices[parameter].squeeze()
             else:
                 mean = choices[parameter][:, orthogonal_component]
@@ -725,7 +725,7 @@ def plot_model_parameters(self, parameter='w_pred', orthogonal_component=1, cros
                 plt.xlabel("Variable No")    
         elif instrument == 'lcms':
             if xaxis is not None and yaxis is not None:
-                _scatterplots(mean, xaxis=xaxis, yaxis=yaxis)
+                _scatterplots(mean, xaxis=xaxis, yaxis=yaxis, marker_size=marker_size)
         else:
             if bar is False:
                 _lineplots(mean, error=error, xaxis=xaxis)
diff --git a/pyChemometrics/ChemometricsPLS.py b/pyChemometrics/ChemometricsPLS.py
@@ -1160,46 +1160,65 @@ def plot_scores(self, comps=[0, 1], color=None, discrete=False, label_outliers=F
             ax.set_ylabel("T[{0}]".format((comps[1] + 1)))
         plt.show()
         return None
-
-    def scree_plot(self, x, y, total_comps=5):
+    
+    def scree_cv(self, x, y, total_comps=5):
         """
 
-        :param x:
+        :param x: Data to use in the scree plot
         :param y:
         :param total_comps:
         :return:
         """
-        fig, ax = plt.subplots()
-
         models = list()
         for n_components in range(1, total_comps + 1):
             currmodel = deepcopy(self)
             currmodel.n_components = n_components
             currmodel.fit(x, y)
             currmodel.cross_validation(x, y)
             models.append(currmodel)
-            q2 = np.array([x.cvParameters['PLS']['Q2Y'] for x in models])
-            r2 = np.array([x.modelParameters['PLS']['R2Y'] for x in models])
 
-        ax.bar([x - 0.1 for x in range(1, total_comps + 1)], height=r2, width=0.2)
-        ax.bar([x + 0.1 for x in range(1, total_comps + 1)], height=q2, width=0.2)
+        q2 = np.array([x.cvParameters['PLS']['Q2Y'] for x in models])
+        r2 = np.array([x.modelParameters['PLS']['R2Y'] for x in models])
+        
+        # Store everything...
+        self.screeCV = {'Q2Y': q2, 'R2Y': r2}
+        
+        return None
+
+    def scree_plot(self, metric = ['Q2Y','R2Y']):
+        """
+
+        :param x:
+        :param y:
+        :param total_comps:
+        :return:
+        """
+        fig, ax = plt.subplots()
+            
+        b1 = self.screeCV[metric[0]]
+        b2 = self.screeCV[metric[1]]
+        
+        total_comps = np.shape(self.screeCV[metric[2]])[0]      
+
+        ax.bar([x - 0.1 for x in range(1, total_comps + 1)], height=b2, width=0.2)
+        ax.bar([x + 0.1 for x in range(1, total_comps + 1)], height=b1, width=0.2)
         ax.legend(['R2', 'Q2'])
         ax.set_xlabel("Number of components")
         ax.set_ylabel("R2/Q2Y")
 
         # Specific case where n comps = 2
-        if q2.size == 2:
-            plateau_index = np.where(np.diff(q2) / q2[0] < 0.05)[0]
+        if b1.size == 2:
+            plateau_index = np.where(np.diff(b1) / b1[0] < 0.05)[0]
             if plateau_index.size == 0:
                 print("Consider exploring a higher level of components")
             else:
-                plateau = np.min(np.where(np.diff(q2)/q2[0] < 0.05)[0])
+                plateau = np.min(np.where(np.diff(b1)/b1[0] < 0.05)[0])
                 ax.vlines(x=(plateau + 1), ymin=0, ymax=1, colors='red', linestyles='dashed')
                 print("Q2Y measure stabilizes (increase of less than 5% of previous value or decrease) "
                       "at component {0}".format(plateau + 1))
 
         else:
-            plateau_index = np.where((np.diff(q2) / q2[0:-1]) < 0.05)[0]
+            plateau_index = np.where((np.diff(b1) / b1[0:-1]) < 0.05)[0]
             if plateau_index.size == 0:
                 print("Consider exploring a higher level of components")
             else:
@@ -1258,7 +1277,7 @@ def plot_permutation_test(self, permt_res, metric='Q2Y'):
             raise exp
     
     ####### flsoares232 version - Updated 20-10-2023
-    def plot_model_parameters(self, parameter='w', component=1, cross_val=False, sigma=2, bar=False, xaxis=None, yaxis=None, instrument=None):
+    def plot_model_parameters(self, parameter='w', component=1, cross_val=False, sigma=2, bar=False, xaxis=None, yaxis=None, instrument=None, marker_size=3):
         
         """
         Plot different model parameters related with the variables
@@ -1317,7 +1336,7 @@ def plot_model_parameters(self, parameter='w', component=1, cross_val=False, sig
                 plt.xlabel("Variable No")    
         elif instrument == 'lcms':
             if xaxis is not None and yaxis is not None:
-                _scatterplots(mean, xaxis=xaxis, yaxis=yaxis)
+                _scatterplots(mean, xaxis=xaxis, yaxis=yaxis, marker_size=marker_size)
         else:
             if bar is False:
                 _lineplots(mean, error=error, xaxis=xaxis)
diff --git a/pyChemometrics/ChemometricsPLSDA.py b/pyChemometrics/ChemometricsPLSDA.py
@@ -1251,16 +1251,15 @@ def _residual_ssx(self, x):
         xscaled = self.scaler.transform(x)
         residuals = np.sum(np.square(xscaled - x_reconstructed), axis=1)
         return residuals
-
-    def scree_plot(self, x, y, total_comps=5):
+    
+    def scree_cv(self, x, y, total_comps=5):
         """
 
         :param x: Data to use in the scree plot
         :param y:
         :param total_comps:
         :return:
         """
-        fig, ax = plt.subplots()
         models = list()
         for n_components in range(1, total_comps + 1):
             currmodel = deepcopy(self)
@@ -1272,34 +1271,71 @@ def scree_plot(self, x, y, total_comps=5):
         q2 = np.array([x.cvParameters['Q2Y'] for x in models])
         r2 = np.array([x.modelParameters['R2Y'] for x in models])
         auc = np.array([x.cvParameters['DA']['Mean_AUC'][0] for x in models])
+        mcc = np.array([x.cvParameters['DA']['Mean_MCC'] for x in models])
+        recall = np.array([x.cvParameters['DA']['Mean_Recall'] for x in models])
+        precision = np.array([x.cvParameters['DA']['Mean_Precision']for x in models])
+        f1 = np.array([x.cvParameters['DA']['Mean_f1'] for x in models])
+        accuracy = np.array([x.cvParameters['DA']['Mean_Accuracy'] for x in models])
+        
+        # Store everything...
+        self.screeCV = {'Q2Y': q2, 'R2Y': r2, 'AUC': auc,
+                        'MCC': mcc, 'Recall': recall,
+                        'Precision': precision, 'F1': f1,
+                        'Accuracy': accuracy}
+        
+        return None
+        
+
+    def scree_plot(self,metric = ['Q2Y','R2Y','AUC']):
+        """
+
+        :param x: Data to use in the scree plot
+        :param y:
+        :param total_comps:
+        :return:
+        """
+        fig, ax = plt.subplots()
+
+        if np.shape(metric)[0] == 1:
+            sup_list = ['Q2Y','R2Y']
+            metric = sup_list + metric
+        elif np.shape(metric)[0] == 2:
+            sup_list = ['Q2Y']
+            metric = sup_list + metric
+            
+        b1 = self.screeCV[metric[0]]
+        b2 = self.screeCV[metric[1]]
+        b3 = self.screeCV[metric[2]]
+        
+        total_comps = np.shape(self.screeCV[metric[2]])[0]
 
-        ax.bar([x - 0.2 for x in range(1, total_comps + 1)], height=r2, width=0.2)
-        ax.bar([x for x in range(1, total_comps + 1)], height=q2, width=0.2)
-        ax.bar([x + 0.2 for x in range(1, total_comps + 1)], height=auc, width=0.2)
+        ax.bar([x - 0.2 for x in range(1, total_comps + 1)], height=b1, width=0.2)
+        ax.bar([x for x in range(1, total_comps + 1)], height=b2, width=0.2)
+        ax.bar([x + 0.2 for x in range(1, total_comps + 1)], height=b3, width=0.2)
 
-        ax.legend(['R2', 'Q2', 'Mean_AUC'])
+        ax.legend(metric)
         ax.set_xlabel("Number of components")
-        ax.set_ylabel("R2/Q2Y/AUC")
+        ax.set_ylabel((metric[0] + "/" + metric[1] + "/" + metric[2]))
 
         # Specific case where n comps = 2 #
-        if q2.size == 2:
-            plateau_index = np.where(np.diff(q2) / q2[0] < 0.05)[0]
+        if b1.size == 2:
+            plateau_index = np.where(np.diff(b1) / b1[0] < 0.05)[0]
             if plateau_index.size == 0:
                 print("Consider exploring a higher level of components")
             else:
-                plateau = np.min(np.where(np.diff(q2)/q2[0] < 0.05)[0])
+                plateau = np.min(np.where(np.diff(b1)/b1[0] < 0.05)[0])
                 ax.vlines(x=(plateau + 1), ymin=0, ymax=1, colors='red', linestyles='dashed')
-                print("Q2Y measure stabilizes (increase of less than 5% of previous value or decrease) "
+                print(metric[0] + " measure stabilizes (increase of less than 5% of previous value or decrease) "
                       "at component {0}".format(plateau + 1))
 
         else:
-            plateau_index = np.where((np.diff(q2) / q2[0:-1]) < 0.05)[0]
+            plateau_index = np.where((np.diff(b1) / b1[0:-1]) < 0.05)[0]
             if plateau_index.size == 0:
                 print("Consider exploring a higher level of components")
             else:
                 plateau = np.min(plateau_index)
                 ax.vlines(x=(plateau + 1), ymin=0, ymax=1, colors='red', linestyles='dashed')
-                print("Q2Y measure stabilizes (increase of less than 5% of previous value or decrease) "
+                print(metric[0] + "measure stabilizes (increase of less than 5% of previous value or decrease) "
                       "at component {0}".format(plateau + 1))
 
         plt.show()
diff --git a/pyChemometrics/plotting_utils.py b/pyChemometrics/plotting_utils.py
