Principes fondamentaux du bootstrapping

Rappelez-vous que pour un ensemble de données bivariées, c'est-à-dire un ensemble de points $X=\left\{(x_i,y_i)\right\}$ dans le plan, la régression linéaire ou droite du meilleur ajustement est une droite d'équation $y=a+bx$, dans laquelle les valeurs de $a$ et $b$ sont choisies de manière à minimiser la distance verticale au carré entre chaque $y_i$ et la valeur prédite, ${\hat{y}}_i=a+b{x}_i$. Les valeurs de $a$ et $b$ dépendent manifestement des points du jeu de données $X$, de sorte que ces paramètres ne sont que nos meilleures estimations des « vraies » valeurs de $a$ et $b$. Un échantillonnage différent de $X$ produirait peut-être des valeurs légèrement différentes de $a$ et $b$. Nous n'avons peut-être pas le luxe de rééchantillonner à partir de la population, et nous faisons donc ce qui s'en rapproche le plus : rééchantillonner à partir de l'ensemble $X$ lui-même ! Ce processus est appelé bootstrapping parce que l'on tente de « se hisser par ses propres bottes », comme le dit le vieux dicton. Voici un aperçu du processus.

1. Choisir un grand nombre $N$ d'échantillons prélevés dans le jeu de données $X$, avec remise (les répétitions de données sont donc autorisées).
2. Déterminer la régression linéaire pour cet ensemble de $N$ échantillons. Enregistrer les valeurs de $a$ (ordonnée à l'origine) et de $b$ (pente) dans des listes $A$ et $B$, respectivement.
3. Répéter les étapes 1 et 2 de nombreuses fois.
4. Les ensembles $A$ et $B$ représentent les distributions des valeurs possibles des paramètres $a$ et $b$ pour le véritable ajustement linéaire, $y=a+bx$.
5. Poser $\bar{a}=MEAN\left(A\right)$ et $\bar{b}=MEAN\left(B\right)$. Ce sont les paramètres que nous choisirons pour notre régression linéaire (obtenue par bootstrapping), $y=\bar{a}+\bar{b}x$.
6. Calculer également l'écart type de $A$ et de $B$ afin de produire un intervalle de confiance pour chaque paramètre.

Ces étapes vont ensuite être illustrées à l'aide de Python.

Régression linéaire avec bootstrapping en Python

Dans le basketball universitaire de division I de la NCAA, une statistique appelée BARTHAG fournit une mesure de nombreux facteurs différents qui influent sur les chances d'une équipe de l'emporter contre une équipe moyenne de la même division. Cette statistique (ainsi que beaucoup d'autres) se trouve dans le jeu de données NCAA-2021-stats.csv. Soit $X$ le jeu de données constitué des paires (BARTHAG, W), où W représente le nombre de victoires en 2021. Cet ensemble comporte 346 points de données. Voici le code qui génère le nuage de points et la droite de régression.

Code Python

      # Import libraries
      import pandas as pd  ## for dataset management
      import matplotlib.pyplot as plt ## for data visualization
      from sklearn.linear_model import LinearRegression
      
      # Read data
      data = pd.read_csv('NCAA-2021-stats.csv')
      data = data[['BARTHAG', 'W']].dropna()
      
      # Find the linear regression model
      model = LinearRegression()
      x = data['BARTHAG'].values.reshape(-1,1)
      y = data['W'].values.reshape(-1,1)
      model.fit(x,y)
      y_pred = model.predict(x)
      
      # Visualizing the linear regression on the scatterplot
      plt.scatter(x,y)
      plt.plot(x, y_pred, color='red')
      plt.show()
      
      a = model.intercept_
      b = model.coef_
      print("Linear model: y = %1.3f + %2.3f x" % (a, b))
    

Le résultat obtenu ressemble à ceci :

    Linear model: y = 4.820 + 14.682 x
    

Ensuite, l'ensemble est rééchantillonné 50 fois, en utilisant 100 points de données à chaque fois. Une régression linéaire est exécutée à chaque fois, produisant un ensemble de 50 valeurs de $a$ et de 50 valeurs de $b$. L'ensemble des 50 droites de régression est représenté sur le même graphique.

Code Python

      from sklearn.utils import resample

      data.plot.scatter(x='BARTHAG', y='W')
      
      # Build arrays for slopes and intercepts
      slopes = []
      intercepts = []
      
      # Resampling 100 points 50 times
      for i in range(50):
        data1 = resample(data, n_samples=100, replace=True)
      
        x = data1['BARTHAG'].values.reshape(-1,1)
        y = data1['W'].values.reshape(-1,1)
        model.fit(x,y)
        y_pred = model.predict(x)
      
        slopes.append(model.coef_[0][0])
        intercepts.append(model.intercept_[0])
      
        plt.plot(x, y_pred, color='red')
      
      plt.show()
    

Le résultat obtenu ressemble à ceci (les résultats pour les modèles de bootstrapping peuvent légèrement varier) :

Remarquez à quel point toutes les droites sont très proches les unes des autres, bien qu'il existe une certaine variation. La moyenne et l'écart type des ordonnées à l'origine et des pentes peuvent être obtenus à l'aide de la bibliothèque statistics de Python.

      import statistics as stats # for mean and standard deviation

      a_bar = stats.mean(intercepts)
      a_std = stats.stdev(intercepts)
      b_bar = stats.mean(slopes)
      b_std = stats.stdev(slopes)
      
      print("Intercept: Mean = %1.3f, St.Dev. = %1.3f " % (a_bar, a_std) )
      print("Slope: Mean = %1.3f, St.Dev. = %1.3f " % (b_bar, b_std) )
    

Intercept: Mean = 4.947, St.Dev. = 0.645
Slope: Mean = 14.489, St.Dev. = 1.267

Rappelez-vous qu'avant le bootstrapping, nous avions obtenu un modèle linéaire avec une ordonnée à l'origine de 4,82 et une pente de 14,682. Ces valeurs sont très proches de l'ordonnée à l'origine et de la pente obtenues à partir de la méthode du bootstrapping. L'avantage du bootstrapping réside toutefois dans le fait que nous disposons désormais de mesures de dispersion (l'écart type), ce qui nous permet de produire un intervalle de confiance pour les deux paramètres. Rappelez-vous que l'intervalle de confiance à 95 % est déterminé à plus ou moins $\left(1.96\right)\sigma /\sqrt{N}$, où $\sigma$ est l'écart type (voir Méthodes de bootstrapping). Cela peut être réalisé en Python à l'aide du paquet scipy.stats.

      import scipy.stats as st
      import numpy as np
      
      a_min, a_max = st.t.interval(0.95, df=len(intercepts)-1, loc=np.mean(intercepts), scale=st.sem(intercepts))
      b_min, b_max = st.t.interval(0.95, df=len(slopes)-1, loc=np.mean(slopes), scale=st.sem(slopes))
      
      print("95%% confidence interval for intercept: %1.2f < a < %1.2f " % (a_min, a_max) )
      print("95%% confidence interval for slope: %1.2f < b < %1.2f " % (b_min, b_max) )
    

95% confidence interval for intercept: 4.76 < a < 5.13
95% confidence interval for slope: 14.13 < b < 14.85

Ainsi, nous nous attendons à ce que la vraie valeur de l'ordonnée à l'origine soit comprise entre 4,76 et 5,13 avec une probabilité de 95 %, et que la vraie valeur de la pente soit comprise entre 14,13 et 14,85 avec une probabilité de 95 %.

Régression linéaire multiple

Un modèle de régression linéaire multiple implique un certain nombre de variables prédictives, que nous noterons $X_1,X_2,\dots X_r$. Nous utilisons ici des lettres majuscules pour distinguer ces variables de caractéristiques des valeurs de données telles que $x_i$. Le modèle de régression multiple prend la forme suivante :

Dans cette équation, $a$ est l'ordonnée à l'origine (valeur de $y$ lorsque tous les $X_k=0$), et les valeurs de $b_k$ sont les coefficients des prédicteurs, qui fonctionnent de manière très similaire à la pente $b$ dans le modèle de régression linéaire simple. Un $b_k$ positif indique que la caractéristique mesurée par $X_k$ a une corrélation positive avec les valeurs de $y$, tandis qu'un $b_k$ négatif indique une corrélation négative. Des modèles plus sophistiqués introduisent également un terme d'erreur, $E$, pour tenir compte de la variabilité inexpliquée de $y$.

L'objectif de la régression linéaire multiple est de déterminer (ou d'estimer) les valeurs des paramètres $a$ et $b_k$ (pour $k=1,2,\dots ,r$) qui minimisent la distance verticale au carré entre les sorties réelles et prédites. Il existe des formules bien connues pour cette tâche, et la régression est généralement traitée par un ordinateur.

Certaines hypothèses doivent être vérifiées pour pouvoir utiliser la régression multiple. Les deux premières concernent la nature même des données d'entrée et de sortie. En règle générale, ces données doivent être continues.

1. La sortie $y$ doit être continue et non discrète. Toutefois, la régression peut toujours être utilisée si la sortie est une variable discrète avec une plage de valeurs très étendue et régulière (comme le nombre d'habitants d'un pays donné, qui peut être n'importe quel entier positif).
2. Les variables d'entrée, $X_k$, doivent être continues et non discrètes. Toutefois, la même remarque qu'au point (1) s'applique ici : si la variable possède une plage de valeurs possibles très étendue et régulière, alors la régression fonctionnera sans problème.

Les deux hypothèses suivantes garantissent que la régression multiple fournira le meilleur modèle. Si ces hypothèses ne sont pas vérifiées, un modèle plus flexible sera nécessaire.

3. Il doit exister une relation linéaire entre $X_k$ et y pour chaque $k$, en considérant les autres variables comme constantes. Cela peut être vérifié en exécutant des régressions linéaires simples sur la paire $(X_k,y)$ pour chaque $k$ et en vérifiant la linéarité.
4. Les données ne doivent pas présenter de multicollinéarité. La multicollinéarité signifie qu'au moins deux des variables d'entrée $X_j$, $X_k$ sont fortement corrélées entre elles, auquel cas une (ou plusieurs) des variables peut être supprimée.

Enfin, un certain nombre d'hypothèses techniques doivent être vérifiées, notamment la recherche de valeurs aberrantes, la vérification que la variance est plus ou moins uniforme sur l'ensemble du domaine (ce que l'on appelle l'homoscédasticité), et l'analyse des résidus pour déterminer s'ils sont relativement indépendants et distribués normalement. Une analyse plus poussée pourrait conduire à déterminer que certaines variables ne sont pas significatives dans le modèle. Ces considérations supplémentaires sortent toutefois du cadre de cet ouvrage.

Régression linéaire multiple en Python

Dans cette section, nous étudions le code Python permettant de réaliser une régression linéaire multiple à l'aide du jeu de données NCAA-2021-stats.csv. Comme dans L'apprentissage automatique dans l'analyse de régression, nous souhaitons prédire le nombre de victoires (W), mais nous inclurons désormais plusieurs caractéristiques prédictives. Dans le premier bloc de code, nous utiliserons BARTHAG et WAB (connu sous le nom de wins above bubble, une mesure du nombre de victoires supplémentaires obtenues par cette équipe par rapport à une équipe typique qui se qualifie de justesse pour le tournoi avec un calendrier similaire) pour prédire W.

      # Import libraries
      import pandas as pd  ## for dataset management
      import matplotlib.pyplot as plt ## for data visualization
      from sklearn.linear_model import LinearRegression
      
      # Read data
      data = pd.read_csv('NCAA-2021-stats.csv')
      data = data.drop(['SEED'], axis=1).dropna() ## Remove the SEED column
      
      X = data[['BARTHAG', 'WAB']]  # Feature list
      y = data[['W']]  # Response variable
      
      # Create the multiple linear regression model
      model = LinearRegression()
      model.fit(X, y)
      
      b = model.intercept_[0]
      a1, a2 = model.coef_[0]
      
      print("Intercept: ", round(b,3))
      print("Coefficients: ", [round(a1,3), round(a2,3)])
      print("Score: ", round(model.score(X,y),2))
    

Intercept:  12.868
Coefficients:  [4.794, 0.545]
Score:  0.58

D'après le résultat, l'équation de régression multiple pour ce modèle est :

Nous constatons que BARTHAG et WAB ont tous deux des coefficients positifs, de sorte que lorsque l'une ou l'autre variable augmente, il y a une augmentation correspondante de W. Le coefficient de détermination R-carré indique qu'environ 58 % de la variance totale de W est attribuable conjointement à BARTHAG et WAB. Examinons les données et observons dans quelle mesure les points se conforment au plan de régression. Tout d'abord, voici le code permettant de générer le nuage de points (BARTHAG, WAB, W) en trois dimensions. (Assurez-vous d'avoir préalablement exécuté le code précédent.)

Code Python

      # import libraries for graphing
      import matplotlib.pyplot as plt
      from mpl_toolkits.mplot3d import Axes3D
      
      X1 = X[['BARTHAG']]
      X2 = X[['WAB']]
      
      # Create a 3D figure
      fig = plt.figure(figsize=(10, 6))
      ax = fig.add_subplot(111, projection='3d')
      
      # Scatterplot of actual data points
      ax.set_xlabel('BARTHAG')
      ax.set_ylabel('WAB', rotation=45)
      ax.set_zlabel('W')
      ax.set_title('Actual Values')
      ax.set_box_aspect(aspect=None, zoom=0.85)
      ax.scatter(X1, X2, y, c='b', marker='o')
    

Le résultat obtenu ressemble à ceci :

Examinons ensuite les valeurs prédites obtenues en injectant BARTHAG et WAB dans la formule de régression multiple présentée précédemment dans ce chapitre. Nous incluons également le plan de régression. Le code Python nécessaire pour créer le plan de régression est relativement complexe, nécessitant un autre paquet, numpy, pour générer une grille de maillage de points $(x,y)$ ; ne vous inquiétez donc pas si vous ne comprenez pas entièrement cette partie du code.

Code Python

      import numpy as np

      # Create a 3D figure
      fig = plt.figure(figsize=(10, 6))
      ax = fig.add_subplot(111, projection='3d')
      
      # Generate prediction values using the model
      y_pred = model.predict(X)
      
      # Create a grid of points for the regression plane
      X1min = min(X1.values)[0]
      X1max = max(X1.values)[0]
      X2min = min(X2.values)[0]
      X2max = max(X2.values)[0]
      X1_space = np.linspace(X1min, X1max, 100)
      X2_space = np.linspace(X2min, X2max, 100)
      x1, x2 = np.meshgrid(X1_space, X2_space)
      
      # Use the regression equation to find z-values
      z = b + a1 * x1 + a2 * x2
      
      # Plot the regression plane
      ax.set_xlabel('BARTHAG')
      ax.set_ylabel('WAB', rotation=45)
      ax.set_zlabel('W')
      ax.set_title('Multiple Linear Regression Plane and Predicted Values')
      ax.set_box_aspect(aspect=None, zoom=0.85)
      ax.scatter(X1, X2, y_pred, c='r', marker='x')
      ax.plot_surface(x1, x2, z, alpha=0.3, cmap='plasma')
    

Le résultat obtenu ressemble à ceci :

Ce n'est peut-être pas tout à fait évident sur le graphique, mais chaque croix rouge « x » est un point situé sur le plan de régression. Pour voir avec quelle précision la formule de régression modélise les données réelles, nous allons retirer les valeurs prédites et insérer les valeurs réelles accompagnées de segments de droite montrant la proximité de ces valeurs par rapport au plan.

Code Python

      # Create a 3D figure
      fig = plt.figure(figsize=(10, 6))
      ax = fig.add_subplot(111, projection='3d')
      
      # Add line segments from actual data points to predicted points
      for i in range(len(X1)):
          ax.plot([X1.values[i,0], X1.values[i,0]], 
           [X2.values[i,0], X2.values[i,0]], 
            [y.values[i,0], y_pred[i,0]], c='g', linewidth=1)
      
      # Plot the regression plane
      ax.set_xlabel('BARTHAG')
      ax.set_ylabel('WAB', rotation=45)
      ax.set_zlabel('W')
      ax.set_title('Multiple Linear Regression')
      ax.set_box_aspect(aspect=None, zoom=0.85)
      ax.scatter(X1, X2, y, c='b', marker='o')
      ax.plot_surface(x1, x2, z, alpha=0.3, cmap='plasma')
    

Le résultat obtenu ressemble à ceci :

Voyons maintenant si nous pouvons améliorer le modèle en ajoutant d'autres caractéristiques. Ajoutons TOR (pourcentage de pertes de balle), 2P_O (pourcentage de réussite aux tirs à deux points) et 3P_O (pourcentage de réussite aux tirs à trois points).

      X = data[['BARTHAG', 'WAB', 'TOR']]  # Feature list
      y = data[['W']]  # Response variable
      
      # Create the multiple linear regression model
      model = LinearRegression()
      model.fit(X, y)
      
      b = model.intercept_[0]
      a1, a2, a3 = model.coef_[0]
      
      print("Intercept: ", round(b,3))
      print("Coefficients: ", [round(a1,3), round(a2,3), round(a3,3)])
      print("Score: ", round(model.score(X,y),2))
    

Intercept:  16.116
Coefficients:  [4.502, 0.526, -0.169]
Score:  0.59

La formule s'étend désormais à :

La mesure du R-carré a légèrement augmenté pour atteindre 0,59, ce qui signifie que le nouveau modèle explique un peu plus de variance. Toutefois, cette faible amélioration de l'exactitude des prédictions se fait au détriment d'une complexité accrue du modèle. Dans la pratique, les modèles plus simples parviennent généralement mieux à éviter le surapprentissage ; nous nous en tiendrons donc aux deux variables de caractéristiques initiales pour cet exemple.

Régression logistique multiple

Comme son nom l'indique, la régression logistique multiple est une extension de la méthode de régression logistique à une seule variable à deux variables d'entrée ou plus. Le modèle estime la probabilité que la variable dépendante $\left(y\right)$ prenne une catégorie particulière telle que oui ou non. Ainsi, le modèle est un classificateur binaire doté de plusieurs entrées indépendantes. Le modèle de régression logistique multiple standard prend la forme suivante :

Ici, la fonction logit est la même que celle décrite dans la section 6.2.1. La régression logistique multiple fonctionne exactement comme la régression logistique simple, en ce sens que la fonction logit est utilisée pour transformer les données, une régression linéaire (multiple) est effectuée, puis les résultats sont retransformés à l'aide d'une fonction logistique. Notez que les caractéristiques peuvent être des données numériques continues ou des données discrètes, y compris des données booléennes (1 ou 0 pour Vrai/Faux ou Oui/Non).

Régression logistique multiple en Python

Revenons au jeu de données CollegeCompletionData.csv de la section La matrice de confusion. Il comporte deux colonnes de caractéristiques, « GPA » et « In-State ». La dernière fois, nous n'avions utilisé qu'une seule caractéristique pour prédire l'obtention du diplôme. Il est tout aussi simple d'utiliser plusieurs caractéristiques !

Code Python

      # Import libraries
      import pandas as pd  ## for dataset management
      import matplotlib.pyplot as plt ## for data visualization
      import seaborn as sns ## for heatmap visualization
      from sklearn.linear_model import LogisticRegression
      from sklearn.metrics import confusion_matrix
      
      # Read data
      data = pd.read_csv('CollegeCompletionData.csv').dropna()
      X = data[['GPA','In_State']]
      y = data['Completion']
      
      # Build the logistic model
      model = LogisticRegression()
      model.fit(X,y)
      
      # Generate model predictions
      y_pred = model.predict(X)
      
      # Generate the confusion matrix
      cf = confusion_matrix(y, y_pred)
      
      # Plot the heatmap using seaborn and matplotlib
      sns.heatmap(cf, annot=True, fmt='d', cmap='Blues', cbar=True)
      plt.xlabel('Predicted')
      plt.ylabel('True')
      plt.yticks(rotation=0)
      plt.title('Confusion Matrix')
      plt.show()
    

Le résultat obtenu ressemble à ceci :

Rappelez-vous que les lignes de la matrice de confusion correspondent aux valeurs réelles, tandis que les colonnes correspondent aux valeurs prédites. Ainsi, il y a 25 vrais positifs (TP) et 23 vrais négatifs (TN), ce qui porte l'exactitude à $48/62=77\%$. Le nombre de faux négatifs (FN) est passé de 10 à 5, mais le nombre de faux positifs (FP) a augmenté de 1, ce qui constitue un compromis tout à fait acceptable.