Régressions linéaires en python

Étape 0 - Partons d'un exemple ! Imaginons un échantillon d'une population, dont on a la taille, le poids et l'IMC (Indice de Masse Corporelle) !

Essayons de voir si on peut retrouver la relation entre taille, poids et IMC sachant qu'on connait pour une fois la réponse IMC = poids/(taille/100)²

Simulons ce jeu de données :

import numpy as np # Importe les fonctions de numpy

from numpy.random import normal as norm # Import juste une fonction de numpy pour simplifier son nom par norm

from numpy import sqrt as sqrt # Import juste une fonction de numpy pour simplifier son nom par sqrt

nb = 180 # Nombre d'individus

poids = np.round(norm(size = nb,loc = 65,scale = 11),0) # Simulation de poids moyens autour de 65 (SD : 11)

imc = norm(size = nb,loc = 22,scale = 5) # Simuations des IMC de ces individus

taille = np.round(sqrt(poids/imc)*100,0) # Déduction d'une taille à partir de ces IMC

taille_bis = np.round(norm(size = nb,loc = 170,scale = 11),0) # Simulation de tailles pour ces individus

taille_biss = np.round(norm(size = nb,loc = 170,scale = 11),0) # Simulation d'autres tailles pour ces individus

taille = np.round((taille+taille_bis+taille_biss)/3,0) # Moyenne des 3 jeux de tailles (2 simulés et 1 déduits d'IMC)

imc = np.round(poids/(taille/100)**2,1) # Calculs des IMC à partir de taille et poids

sexe = ["F"]*90+["M"]*90

num = list(range(1,181,1)) # numérotation de 1 à 180

import pandas as df # Importation de pandas pour la gestion des tableaux de données (dataframes)

from pandas import DataFrame as DataFrame # Importation uniquement de la fonction DataFrame pour réduire son écriture

compil = {'num':num,'poids':poids,'taille':taille,'imc':imc,'sexe':sexe}

compil = DataFrame(compil) # Compilation des valeurs en dataframe (il faut que ce soit en dictionnaire ci-avant)

En projet : développer une accélération de la fonction et les intervalles de confiance en polygones grisés.

Développement en cours : suivre ces liens externes :

• https://stackoverrun.com/fr/q/12785742

• source code de l'image 2

• Je peux faire une matrice de corrélation

mycor = compil.corr() ; print(mycor)

• Pour faire un test de corrélation précis avec p-value (l'équivalent de cor.test() de R), il faut faire appel au modul scipy

import scipy

my_cor, pval = scipy.stats.pearsonr(poids,imc)

# même si ma corrélation est faible (valeur <<1), la valeur 2 montre une p-value très significative

print(my_cor)

print(pval)

Dans une relation simple où y dépend d'un seul x (y = f(x)), la corrélation correspond à la racine carrée du R².

A taper dans la console windows pour installer des librairies utiles (Keras est optionnel, à développer sur ce site)

python -m pip install sklearn Keras

python -m pip install --upgrade pip

pip3 install sklearn

Etape 2.1 - Définir les valeurs X et Y (attention, tous les formats ne sont pas acceptés)

Option 1 - Les valeurs X sont des arrays de type numpy. X doit être verticalisé avec la fonction reshape().

X = (np.array(compil.poids)).reshape(-1,1)

Y = compil.imc

Option 2- Même si X est une colonne extraite d'une dataframe, il faut forcer sa reconnaissance dans ce format

X = DataFrame(compil.poids)

Y = compil.imc

Etape 2.2 - Mettre en place le modèle de régression

from sklearn.linear_model import LinearRegression

myreg = LinearRegression()

myreg = myreg.fit(X,Y)

Autre méthode se rapprochant de l'écriture de R et offrant un bilan plus complet

import statsmodels.formula.api as sm

model = sm.ols(formula='imc ~ poids',data=compil)

# model = sm.ols(formula='imc ~ poids +0',data=compil) # Même régression passant par 0

myreg = model.fit()

print(myreg.summary())

Etape 2.4 - Récupérer l'équation de la droite

print(myreg.intercept_) # B0 (b)

print(myreg.coef_) # B1 (coefficient directeur)

print(myreg._residues) # A étudier pour avoir tous les résidus. sinon on peut faire la différence entre Y et myreg.predict(X)

print('coef dir =',myreg.coef_, 'ord orgine =',myreg.intercept_)

La première chose à vérifier est que le modèle à une bonne p-value (significatif) et que ses coefficients ont aussi des bonnes p-values.

Il faut aussi contrôler les p-values des coefficients et du modèle.

• La p-value du modèle se calcule en faisant une simple corrélation entre y et fitted.values.

• Les p-values des coefficients sont à extraire de diverses façon selon statsmodels et sklearn

statsmodels :

myreg.pvalues

# ou encore :

for attributeIndex in range (0, 2):

print(myreg.pvalues[attributeIndex])

sklearn

A voir...

Les résidus doivent avoir à peu près la même valeur moyenne pour les différentes valeurs de y prédites (fitted values).

On peut vérifier que cette condition est respectée par une 1) représentation graphique ou 2) par un test statistique (test de Rainbow).

2) Test statistique de Rainbow

Le test de Rainbow vérifie l'hypothèse nulle suivante : la représentation statistique est bien linéaire. On parle ici de contrôle de l'adéquation. La p-value renvoyée par ce test devrait donc être supérieure à 0,05 pour qu'on considère que le modèle de régression peut être conservé.

from statsmodels.stats.diagnostic import linear_rainbow

Ftest, pval = linear_rainbow(myreg)

print(pval)

Il faut ensuite contrôler l'indépendance des résidus. On peut le faire là-aussi via 1) une représentation graphique ACP (auto-correlation function) ou 2) par le test statistique de Durbin-Watson.

2) Test de Durbin-Watson

Le test de Durbin-Watson (un exemple interne pour en savoir plus) permet de détecter ces autocorrélations par un test objectif. La valeur renvoyée doit s'approcher de 2. Plus elle s'en éloigne, plus on risque d'avoir des autocorrélations.

Seuils admis (et très subjectifs) : pas d'autocorrélation entre 1.6 et 2.4... (quoi que...)

from statsmodels.stats.stattools import durbin_watson

# Durbin-Watson test

print(durbin_watson(myreg.resid))

Les résidus standards doivent être répartis selon la loi normale. On peut vérifier cela 1) par des graphiques (histogramme ou QQ-plot) ou 2) par des tests statistiques (test Shapiro-Wilk ou test de Jarque-Berra).

Récupérons les résidus :

• A partir d'une régression statsmodels

residues = myreg.resid_pearson

• A partir d'une régression sklearn

residues = myreg.resid # Y-myreg.predict(X)

Remarque : il faut ici utiliser les résidus standards : chaque résidu doit être divisé par la racine carrée de la somme des carrés des résidus divisée par N.

residues_std = residues/sqrt(sum(residues**2)/(len(residues)-1))

1) Contrôle graphique par histogramme ou QQ-plot

2) Contrôle par un test statistique

Le test de Shapiro (qui fonctionne pour les petits échantillons inférieurs à 50) permet de contrôler l'hypothèse nulle d'une répartition normale. Si p-value > 0.05 <-> Normalité.

from scipy.stats import shapiro

x, pval = shapiro(residues_std) ; print(pval)

Pour les grands échantillons, on peut aussi faire de test de Jarque-Berra :

from scipy.stats import jarque_bera

x, pval = jarque_bera(residues_std) ; print(pval)

Un autre contrôle permet de vérifier que la variance des résidus est constante (cas d'homoscédasticité).

Là encore, 1) une représentation graphique suffit ou 2) le test statistique de Breusch-Pagan.

On va cette fois-ci raisonner sur les résidus standardisés, la racine carrée des résidus standards.

sqrt_residues_std = sqrt(abs(residues_std)) # ne fonctionne pas si on ne met pas abs

2) Test de Breusch-Pagan

Cet test renvoie une p-value qui interroge la variance des résidus. Si la p-value est inférieure à 0,05, on doit rejeter la régression comme présentant de l'hétéroscédasticité.

from statsmodels.stats.api import het_breuschpagan

lagrande, pval, f, fpval = het_breuschpagan(myreg.resid, myreg.model.exog)

On notera qu'on ne donne pas à avaler les résidus standardisés à cette fonction. Elle se charge de mener les bons calculs.

myreg.model.exog correspond à une array qui regroupe les indices des résidus les plus susceptibles d'induire une hétéroscédasticité.

• En savoir plus sur l'analyse de la pertinence d'une régression avec statsmodels

• Aller encore plus loin

• Pour faire comme en R avec un code détaillé

On vient de voir de nombreux graphiques qui permettent de valider un modèle de régression.

Ce n'est pas toujours simples de les avoir sous python, alors, on peut émuler R pour les avoir grâce au module Rpy2.

Sinon, j'envisage de compiler une fonction valreg() qui pourrait faire toutes les validations de façon automatique. Ce serait pour le moins, pratique.

• En lien, comment émuler R sous python avec le module Rpy2

• En lien, comment interpréter les graphiques qui dont le bilan d'un modèle de Régression : Régression avec R (Aide à l'utilisation de R)

Faire une régression multiple suit la même logique qu'une régression simple.

Etape 1 - Définir les variables X

Option 1 - Définir les variables X que l'on veut dans la dataframe en supprimant celles que l'on ne veut pas avec drop()

list_var=compil.columns.drop(["imc","num","sexe"]) ; X=compil[list_var] ; print(X)

Option 2 - Définir les noms des variables X que l'on désire

X = compil.loc[:,["poids","taille"]] ; print(X)

Option 3 - Définir les variables X que l'on désire en les appelant par son index

X = compil.iloc[:,1:3] ; print(X)

Etape 2 - Définir la/les variables Y

Y = compil.imc

Etape 3 - Faire le modèle de régression

myreg = LinearRegression().fit(X,Y)

print('coef dir =',myreg.coef_, 'ord orgine =',myreg.intercept_)

• Approche basique

import statsmodels.api as sm

mod = sm.OLS(compil.imc, compil.iloc[:,1:3])

# Ajuster le modèle

lmfit = mod.fit() # Résumé des principales

print(lmfit.summary())

• Approche avec des formules

import statsmodels.api as sm

import statsmodels.formula.api as smf

print(compil) # données simulées dans le paragraphe 0 ci-dessus

# Définition du modèle

import statsmodels.formula.api as smf

model = smf.ols(formula='imc ~ poids + taille',data=compil)

myreg = model.fit()

print(myreg .summary())

• Si message d'erreur : copier coller dans la console windows

pip install --upgrade patsy

pip install statsmodels --upgrade

pip install pandas --upgrade

Analyser certains paramètres tels l'AIC et le BIC qui doivent être le plus bas possible ! (cf. summary())

Avec statsmodels

myreg.bic

myreg.aic

Avec sklearn (à contrôler... )

from sklearn.linear_model import LassoLarsIC

model_bic = LassoLarsIC(criterion="bic", normalize=False)

model_bic.fit(X, y)

alpha_bic_ = model_bic.alpha_ # Coefficients de pénalité

model_aic = LassoLarsIC(criterion="aic", normalize=False)

model_aic.fit(X, y)

alpha_aic_ = model_aic.alpha_ # Coefficients de pénalité

Je n'ai pas pu encore fait une comparaison sur la solidité de l'approche sklearn, cette fonction étant adapté normalement à Lasso et Lars, à vérifier donc.

Rubrique en cours de développement. cf. la démarche sous R pour équivalence.

Pour valider un modèle, il existe plusieurs approches :

1. Validation avec échantillon d'apprentissage et échantillon test - on s'assure alors que le MSE soit bas. Et que le modèle marche pour l'échantillon test.

2. Validation croisée : un bootstrap permet de refaire plusieurs fois la méthode ci-dessus.

• Sous sklearn ou statsmodels

Y_pred = myreg.predict(X)

Y_pred = myreg.predict(compil[["poids","taille"]]) # Dans un exemple ou le modèle de prédiction de Y aurait été conçu à partir des variables poids et taille.

Reprenons les données du haut de la page et faisons une régression de l'IMC en fonction de deux variables X1 et X2 : la Taille et le Poids

Etape 1 - Récupérons à partir des deux variables des valeurs couvrant l'ensemble d'une grille

variable1 = compil.poids # à modifier

variable2 = compil.taille # à modifier

res = 40 # résolution à paramétrer soi-même : plus c'est grand, plus c'est joli mais lourd

X1 = np.arange(min(variable1 ),max(variable1 ),((max(variable1 )-min(variable1 ))/res))

X2 = np.arange(min(variable2 ),max(variable2 ),((max(variable2 )-min(variable2 ))/res))

Etape 2 -Je combine les deux étendues de X1 et X2 pour créer une Grille de valeurs possibles

# Une fonction utile

def expand_grid(x, y,name_x1="X1",name_x2="X2"):

x1, x2 = np.meshgrid(x, y) # create the actual grid

x1G = x1.flatten() # make the grid 1D

x2G = x2.flatten() # conversion 1D

return x1, x2, pd.DataFrame({name_x1:x1G,name_x2:x2G})

X1, X2, Grille=expand_grid(X1,X2,name_x1="poids",name_x2="taille")

Etape 3 - J'utilise la commande predict() pour simuler toutes les valeurs d'IMC possible pour les combinaisons de X1 et X2

out = myreg.predict(Grille)

# Je remets en 2D les prédictions

out = np.array(out).reshape(res,res)

Etape 4 -J'affiche le graphique en 3D

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()

ax = fig.gca(projection='3d')

surface = ax.plot_surface(X1,X2,out, rstride=1, cstride=1, cmap="plasma",

linewidth=0, antialiased=False)

ax.set_zlim(17, 30) # A paramétrer

fig.colorbar(surface, shrink=0.5, aspect=5)

plt.show()