.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_meta_modeling/general_purpose_metamodels/plot_stepwise.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_meta_modeling_general_purpose_metamodels_plot_stepwise.py: Perform stepwise regression =========================== In this example we perform the selection of the most suitable function basis for a linear regression model with the help of the stepwise algorithm. We consider te so-called Linthurst data set, which contains measures of aerial biomass (BIO) as well as 5 five physicochemical properties of the soil: salinity (SAL), pH, K, Na, and Zn. The data set is taken from the book [rawlings2001]_ and is provided below: .. GENERATED FROM PYTHON SOURCE LINES 13-19 .. code-block:: Python import openturns as ot from openturns.viewer import View import numpy as np import matplotlib.pyplot as plt from openturns.usecases import linthurst .. GENERATED FROM PYTHON SOURCE LINES 20-22 We define the data. .. GENERATED FROM PYTHON SOURCE LINES 24-43 .. code-block:: Python ds = linthurst.Linthurst() dimension = ds.data.getDimension() - 1 input_sample = ds.data[:, 1: dimension + 1] print("Input :") print(input_sample[:5]) output_sample = ds.data[:, 0] print("Output :") print(output_sample[:5]) input_description = input_sample.getDescription() output_description = output_sample.getDescription() n = input_sample.getSize() print("n = ", n) dimension = ds.data.getDimension() - 1 print("dimension = ", dimension) .. rst-class:: sphx-glr-script-out .. code-block:: none Input : [ SAL pH K Na Zn ] 0 : [ 33 5 1441.67 35185.5 16.4524 ] 1 : [ 35 4.75 1299.19 28170.4 13.9852 ] 2 : [ 32 4.2 1154.27 26455 15.3276 ] 3 : [ 30 4.4 1045.15 25072.9 17.3128 ] 4 : [ 33 5.55 521.62 31664.2 22.3312 ] Output : [ BIO ] 0 : [ 676 ] 1 : [ 516 ] 2 : [ 1052 ] 3 : [ 868 ] 4 : [ 1008 ] n = 45 dimension = 5 .. GENERATED FROM PYTHON SOURCE LINES 44-54 Complete linear model --------------------- We consider a linear model with the purpose of predicting the aerial biomass as a function of the soil physicochemical properties, and we wish to identify the predictive variables which result in the most simple and precise linear regression model. We start by creating a linear model which takes into account all of the physicochemical variables present within the Linthrust data set. Let us consider the following linear model :math:`\tilde{Y} = a_0 + \sum_{i = 1}^{d} a_i X_i + \epsilon`. If all of the predictive variables are considered, the regression can be performed with the help of the :class:`~openturns.LinearModelAlgorithm` class. .. GENERATED FROM PYTHON SOURCE LINES 56-62 .. code-block:: Python algo_full = ot.LinearModelAlgorithm(input_sample, output_sample) algo_full.run() result_full = ot.LinearModelResult(algo_full.getResult()) print("R-squared = ", result_full.getRSquared()) print("Adjusted R-squared = ", result_full.getAdjustedRSquared()) .. rst-class:: sphx-glr-script-out .. code-block:: none R-squared = 0.677310820565376 Adjusted R-squared = 0.6359404129455524 .. GENERATED FROM PYTHON SOURCE LINES 63-70 Forward stepwise regression --------------------------- We now wish to perform the selection of the most important predictive variables through a stepwise algorithm. It is first necessary to define a suitable function basis for the regression. Each variable is associated to a univariate basis and an additional basis is used in order to represent the constant term :math:`a_0`. .. GENERATED FROM PYTHON SOURCE LINES 72-77 .. code-block:: Python functions = [] functions.append(ot.SymbolicFunction(input_description, ["1.0"])) for i in range(dimension): functions.append(ot.SymbolicFunction(input_description, [input_description[i]])) basis = ot.Basis(functions) .. GENERATED FROM PYTHON SOURCE LINES 78-81 Plese note that this example uses a linear basis with respect to the various predictors for the sake of clarity. However, this is not a necessity, and more complex and non linear relations between predictors may be considered (e.g., polynomial bases). .. GENERATED FROM PYTHON SOURCE LINES 84-90 We now perform a forward stepwise regression. We suppose having no information regarding the given data set, and therefore the set of minimal indices only contains the constant term (indexed by 0). The first regression is performed by relying on the Akaike Information Criterion (AIC), which translates into a penalty term equal to 2. In practice, the algorithm selects the functional basis subset that minimizes the AIC by iteratively adding the single function which provides the largest improvement until convergence is reached. .. GENERATED FROM PYTHON SOURCE LINES 92-105 .. code-block:: Python minimalIndices = [0] direction = ot.LinearModelStepwiseAlgorithm.FORWARD penalty = 2.0 algo_forward = ot.LinearModelStepwiseAlgorithm( input_sample, output_sample, basis, minimalIndices, direction ) algo_forward.setPenalty(penalty) algo_forward.run() result_forward = algo_forward.getResult() print("Selected basis: ", result_forward.getCoefficientsNames()) print("R-squared = ", result_forward.getRSquared()) print("Adjusted R-squared = ", result_forward.getAdjustedRSquared()) .. rst-class:: sphx-glr-script-out .. code-block:: none Selected basis: [[SAL,pH,K,Na,Zn]->[1.0],[SAL,pH,K,Na,Zn]->[pH],[SAL,pH,K,Na,Zn]->[Na]] R-squared = 0.658432822226285 Adjusted R-squared = 0.6421677185227748 .. GENERATED FROM PYTHON SOURCE LINES 106-113 With this first forward stepwise regression, the results show that the selected optimal basis contains a constant term, plus two linear terms depending respectively on the pH value (pH) and on the sodium concentration (Na). As can be expected, the R-squared value diminishes if compared to the regression on the entire basis, as the stepwise regression results in a lower number of predictive variables. However, it can also be seen that the adjusted R-squared, which is a metric that also takes into account the ratio between the amount of training data and the number of explanatory variables, is improved if compared to the complete model. .. GENERATED FROM PYTHON SOURCE LINES 115-121 Backward stepwise regression ---------------------------- We now perform a backward stepwise regression, meaning that rather than iteratively adding predictive variables, we will be removing them, starting from the complete model. This regression is performed by relying on the Bayesian Information Criterion (BIC), which translates into a penalty term equal to :math:`log(n)`. .. GENERATED FROM PYTHON SOURCE LINES 123-136 .. code-block:: Python minimalIndices = [0] direction = ot.LinearModelStepwiseAlgorithm.BACKWARD penalty = np.log(n) algo_backward = ot.LinearModelStepwiseAlgorithm( input_sample, output_sample, basis, minimalIndices, direction ) algo_backward.setPenalty(penalty) algo_backward.run() result_backward = algo_backward.getResult() print("Selected basis: ", result_backward.getCoefficientsNames()) print("R-squared = ", result_backward.getRSquared()) print("Adjusted R-squared = ", result_backward.getAdjustedRSquared()) .. rst-class:: sphx-glr-script-out .. code-block:: none Selected basis: [[SAL,pH,K,Na,Zn]->[1.0],[SAL,pH,K,Na,Zn]->[pH],[SAL,pH,K,Na,Zn]->[K]] R-squared = 0.6475759074104157 Adjusted R-squared = 0.6307938077632926 .. GENERATED FROM PYTHON SOURCE LINES 137-139 It is interesting to point out that although both approaches converge towards a model characterized by 2 predictive variables, the selected variables do not coincide. .. GENERATED FROM PYTHON SOURCE LINES 141-146 Both directions stepwise regression ----------------------------------- A third available option consists in performing a stepwise regression in both directions, meaning that at each iteration the predictive variables can be either added or removed. In this case, a set of starting indices must be provided to the algorithm. .. GENERATED FROM PYTHON SOURCE LINES 148-162 .. code-block:: Python minimalIndices = [0] startIndices = [0, 2, 3] penalty = np.log(n) direction = ot.LinearModelStepwiseAlgorithm.BOTH algo_both = ot.LinearModelStepwiseAlgorithm( input_sample, output_sample, basis, minimalIndices, direction, startIndices ) algo_both.setPenalty(penalty) algo_both.run() result_both = algo_both.getResult() print("Selected basis: ", result_both.getCoefficientsNames()) print("R-squared = ", result_both.getRSquared()) print("Adjusted R-squared = ", result_both.getAdjustedRSquared()) .. rst-class:: sphx-glr-script-out .. code-block:: none Selected basis: [[SAL,pH,K,Na,Zn]->[1.0],[SAL,pH,K,Na,Zn]->[pH],[SAL,pH,K,Na,Zn]->[K]] R-squared = 0.6475759074104157 Adjusted R-squared = 0.6307938077632926 .. GENERATED FROM PYTHON SOURCE LINES 163-165 It is interesting to note that the basis varies depending on the selected set of starting indices, as is shown below. An informed initialization might therefore improve the model selection and the resulting regression .. GENERATED FROM PYTHON SOURCE LINES 167-181 .. code-block:: Python minimalIndices = [0] startIndices = [0, 1] penalty = np.log(n) direction = ot.LinearModelStepwiseAlgorithm.BOTH algo_both = ot.LinearModelStepwiseAlgorithm( input_sample, output_sample, basis, minimalIndices, direction, startIndices ) algo_both.setPenalty(penalty) algo_both.run() result_both = algo_both.getResult() print("Selected basis: ", result_both.getCoefficientsNames()) print("R-squared = ", result_both.getRSquared()) print("Adjusted R-squared = ", result_both.getAdjustedRSquared()) .. rst-class:: sphx-glr-script-out .. code-block:: none Selected basis: [[SAL,pH,K,Na,Zn]->[1.0],[SAL,pH,K,Na,Zn]->[pH],[SAL,pH,K,Na,Zn]->[Na]] R-squared = 0.658432822226285 Adjusted R-squared = 0.6421677185227748 .. GENERATED FROM PYTHON SOURCE LINES 182-187 Graphical analyses ------------------ Finally, we can rely on the LinearModelAnalysis class in order to analyse the predictive differences between the obtained models. .. GENERATED FROM PYTHON SOURCE LINES 189-208 .. code-block:: Python analysis_full = ot.LinearModelAnalysis(result_full) analysis_full.setName("Full model") analysis_forward = ot.LinearModelAnalysis(result_forward) analysis_forward.setName("Forward selection") analysis_backward = ot.LinearModelAnalysis(result_backward) analysis_backward.setName("Backward selection") fig = plt.figure(figsize=(12, 8)) for k, analysis in enumerate([analysis_full, analysis_forward, analysis_backward]): graph = analysis.drawModelVsFitted() ax = fig.add_subplot(3, 1, k + 1) ax.set_title(analysis.getName(), fontdict={"fontsize": 16}) graph.setXTitle("Exact values") ax.xaxis.label.set_size(12) ax.yaxis.label.set_size(14) graph.setTitle("") v = View(graph, figure=fig, axes=[ax]) plt.tight_layout() .. image-sg:: /auto_meta_modeling/general_purpose_metamodels/images/sphx_glr_plot_stepwise_001.png :alt: , Full model, Forward selection, Backward selection :srcset: /auto_meta_modeling/general_purpose_metamodels/images/sphx_glr_plot_stepwise_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 209-213 For illustrative purposes, we show the Bayesian Information Criterion (BIC) and Akaike Information Criterion (AIC) values which are computed during the iterations of the forward step-wise regression. Please note that in order to do so, we set the penalty parameter to a negligible value (meaning that the basis selection only takes into account the model likelihood, and not the number of parameters characterizing the linear model). .. GENERATED FROM PYTHON SOURCE LINES 215-248 .. code-block:: Python minimalIndices = [0] penalty = 1e-10 direction = ot.LinearModelStepwiseAlgorithm.FORWARD BIC = [] AIC = [] for iterations in range(1, 6): algo_forward = ot.LinearModelStepwiseAlgorithm( input_sample, output_sample, basis, minimalIndices, direction ) algo_forward.setPenalty(penalty) algo_forward.setMaximumIterationNumber(iterations) algo_forward.run() result_forward = algo_forward.getResult() RSS = np.sum( np.array(result_forward.getSampleResiduals()) ** 2 ) # Residual sum of squares LL = n * np.log(RSS / n) # Log-likelihood BIC.append(LL + iterations * np.log(n)) # Bayesian Information Criterion AIC.append(LL + iterations * 2) # Akaike Information Criterion print("Selected basis: ", result_forward.getCoefficientsNames()) plt.figure() plt.plot(np.arange(1, 6), BIC, label="BIC") plt.plot(np.arange(1, 6), AIC, label="AIC") plt.xticks(np.arange(1, 6)) plt.xlabel("Basis size", fontsize=14) plt.ylabel("Information criterion", fontsize=14) plt.legend(fontsize=14) plt.tight_layout() .. image-sg:: /auto_meta_modeling/general_purpose_metamodels/images/sphx_glr_plot_stepwise_002.png :alt: plot stepwise :srcset: /auto_meta_modeling/general_purpose_metamodels/images/sphx_glr_plot_stepwise_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Selected basis: [[SAL,pH,K,Na,Zn]->[1.0],[SAL,pH,K,Na,Zn]->[pH]] Selected basis: [[SAL,pH,K,Na,Zn]->[1.0],[SAL,pH,K,Na,Zn]->[pH],[SAL,pH,K,Na,Zn]->[Na]] Selected basis: [[SAL,pH,K,Na,Zn]->[1.0],[SAL,pH,K,Na,Zn]->[pH],[SAL,pH,K,Na,Zn]->[Na],[SAL,pH,K,Na,Zn]->[Zn]] Selected basis: [[SAL,pH,K,Na,Zn]->[1.0],[SAL,pH,K,Na,Zn]->[pH],[SAL,pH,K,Na,Zn]->[Na],[SAL,pH,K,Na,Zn]->[Zn],[SAL,pH,K,Na,Zn]->[SAL]] Selected basis: [[SAL,pH,K,Na,Zn]->[1.0],[SAL,pH,K,Na,Zn]->[pH],[SAL,pH,K,Na,Zn]->[Na],[SAL,pH,K,Na,Zn]->[Zn],[SAL,pH,K,Na,Zn]->[SAL],[SAL,pH,K,Na,Zn]->[K]] .. GENERATED FROM PYTHON SOURCE LINES 249-252 The graphic above shows that the optimal linear model in terms of compromise between prediction likelihood and model complexity should take into account the influence of 2 regession variables as well as the constant term. This is coherent with the results previously obtained .. _sphx_glr_download_auto_meta_modeling_general_purpose_metamodels_plot_stepwise.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_stepwise.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_stepwise.py `