Régression linéaire

Python works with packages. A large number of packages can be dowloaded all together in the Anaconda environment. Start downloading here (version >3) https://www.anaconda.com/download/ Anaconda comes with

  • Spyder, a GUI for python
  • Jupyter Notebook, how to create great documents

Enjoy!

Bee careful :

  • Python is sensible to spaces in the beginning of lines,
  • Python is sensible to number type : $3$ and $3.$ are totally different: integer and float types.
In [1]:
import warnings
warnings.filterwarnings('ignore')

import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
import pandas as pd

Les données

Le fichier Supra_data_s.txt contient les données relevées lors d’une étude prospective (SUPRA - Supervision of the PRogression of Atherosclerosis in HIV-1 infected patients) parmi des participants de la cohorte Aquitaine. L’étude portait sur l’évaluation des facteurs de risques associés aux maladies cardiovasculaires chez les patients infectés par le VIH, en utilisant le marqueur IMT - intima media thickness of the carotid arteries (colonne mesure). L’étude comportait le relevé d’autres marqueurs connus comme étant des facteurs de risques des maladies cardiovasculaires : sexe (SEXE = 1 si homme, 2 si femme), age (AGE en années), taille (en cm), poids (en kg), habitudes liées au sport (0=non sportif, 1=sportif), tabac et prise d’alcool (0=consommation nulle, 1=consommation occasionnelle, 2=consommation régulière).

In [2]:
data = np.loadtxt("Supra_data_s.txt",skiprows=1)
header = ["SEXE","AGE","taille","poids","tabac","SPORT","mesure","alcool"]
print(data[0:10,:])

[[  1.    33.   170.    70.     1.     0.     0.52   1.  ]
 [  2.    33.   177.    67.     2.     0.     0.42   1.  ]
 [  2.    53.   164.    63.     1.     0.     0.65   0.  ]
 [  2.    42.   169.    76.     1.     1.     0.48   1.  ]
 [  2.    53.   152.    54.     0.     0.     0.45   1.  ]
 [  2.    50.   162.    53.     2.     0.     0.49   1.  ]
 [  1.    22.   181.    72.     0.     1.     0.42   1.  ]
 [  2.    26.   162.    61.     1.     1.     0.45   1.  ]
 [  1.    50.   170.    80.     0.     0.     0.65   2.  ]
 [  2.    47.   178.    73.     0.     0.     0.52   1.  ]]

Creation du BMI

Créer la variable BMI Indice de Masse Corporelle telle que $$ BMI=\frac{Poids}{Taille^2}, $$ et l'ajouter au dataset.

In [3]:
bmi = data[:,3:4]/(data[:,2:3]/100)**2
data = np.concatenate((data,bmi),axis=1)
header = header + ["bmi"]
print(data.shape)

(110, 9)

La régression linéaire

QUESTION Visualiser :

  • IMT en fonction de BMI
  • IMT en fonction de age
  • IMT en fonction de taille
  • IMT en fonction de poids
In [25]:
plt.figure(1)
plt.subplot(221)
plt.scatter(data[:,8],data[:,6])
plt.xlabel("BMI")
plt.ylabel("IMT")

plt.subplot(222)
plt.scatter(data[:,1],data[:,6])
plt.xlabel("age")
plt.ylabel("IMT")

plt.subplot(223)
plt.scatter(data[:,2],data[:,6])
plt.xlabel("taille")
plt.ylabel("IMT")

plt.subplot(224)
plt.scatter(data[:,3],data[:,6])
plt.xlabel("poids")
plt.ylabel("IMT")

plt.subplots_adjust(top=0.92, bottom=-0.3, left=0.10, right=0.95, hspace=0.25,wspace=0.35)

QESTION Quel est selon vous la variable qui permette le mieux de prédire, avec un modèle linéaire la variable IMT ?

QUESTION Pour chacun des couples vu précédemment, calculer le coefficient de corrélation de Pearson. Et commenter en fonction de ce qui a été vu dans la question précédente.

In [5]:
for i in [8,1,2,3]:
    print(np.corrcoef(data[:,6],data[:,i])[0,1])

0.30717287224763695
0.552366859681787
-0.06865902433560454
0.22952003414114933

QUESTION Utiliser la fonction OLS de la librairie statsmodels.api afin d'estimer le modèle de régression prédisant l'IMT en fonction des autres variables du problème. Le modèles sera noté model, la matrice des coefficients de régression sera notée $B$ et l'intercept (=ordonnée à l'origine) $B_0$. On fera attention à inclure l'intercept dans le modèle.

In [6]:
import statsmodels.api as sm

# Delete 6th column of data and put it in the matrix X
# sm.add_constant permits to include the computation of the intercept in the model :) do not forget it!
X = sm.add_constant(np.delete(data,6,1))
Y = data[:,6:7]

model = sm.OLS(Y, X).fit()

B_0 = model.params[0]
B = np.delete(np.array([model.params]).T,0,0)

QUESTION Observer le détail du modèle linéaire construit et interprêter les différentes p-values.

In [7]:
print(model.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.398
Model:                            OLS   Adj. R-squared:                  0.350
Method:                 Least Squares   F-statistic:                     8.337
Date:                Wed, 14 Nov 2018   Prob (F-statistic):           1.26e-08
Time:                        14:04:36   Log-Likelihood:                 141.35
No. Observations:                 110   AIC:                            -264.7
Df Residuals:                     101   BIC:                            -240.4
Df Model:                           8                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.8571      0.794      2.339      0.021       0.282       3.432
x1             0.0078      0.019      0.420      0.675      -0.029       0.045
x2             0.0041      0.001      5.571      0.000       0.003       0.006
x3            -0.0097      0.005     -2.049      0.043      -0.019      -0.000
x4             0.0137      0.006      2.317      0.023       0.002       0.025
x5             0.0040      0.009      0.439      0.662      -0.014       0.022
x6            -0.0007      0.014     -0.052      0.958      -0.028       0.027
x7             0.0060      0.012      0.481      0.631      -0.019       0.031
x8            -0.0334      0.016     -2.033      0.045      -0.066      -0.001
==============================================================================
Omnibus:                       13.948   Durbin-Watson:                   1.962
Prob(Omnibus):                  0.001   Jarque-Bera (JB):               17.717
Skew:                           0.681   Prob(JB):                     0.000142
Kurtosis:                       4.418   Cond. No.                     2.24e+04
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.24e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

QUESTION Quelles sont les variables significatives ?

Ce sont :

  • const soit l'intercept
  • x2 soit AGE
  • x3 soit taille
  • x4 soit poids
  • x8 soit IMC

QUESTION Créer un nouveau modèle avec ces variables et observer les éventuelles modifications.

In [26]:
X = sm.add_constant(data[:,[1,2,3,8]])
Y = data[:,6:7]

model_sparse = sm.OLS(Y, X).fit()

B_0 = model_sparse.params[0]
B = np.delete(np.array([model_sparse.params]).T,0,0)

print(model_sparse.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.395
Model:                            OLS   Adj. R-squared:                  0.371
Method:                 Least Squares   F-statistic:                     17.10
Date:                Wed, 14 Nov 2018   Prob (F-statistic):           7.89e-11
Time:                        14:14:07   Log-Likelihood:                 141.05
No. Observations:                 110   AIC:                            -272.1
Df Residuals:                     105   BIC:                            -258.6
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.8393      0.771      2.386      0.019       0.311       3.368
x1             0.0043      0.001      6.676      0.000       0.003       0.006
x2            -0.0095      0.005     -2.068      0.041      -0.019      -0.000
x3             0.0133      0.006      2.344      0.021       0.002       0.025
x4            -0.0325      0.016     -2.046      0.043      -0.064      -0.001
==============================================================================
Omnibus:                       12.679   Durbin-Watson:                   1.958
Prob(Omnibus):                  0.002   Jarque-Bera (JB):               15.569
Skew:                           0.640   Prob(JB):                     0.000416
Kurtosis:                       4.326   Cond. No.                     2.21e+04
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.21e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

QUESTION Représenter les valeurs prédites en fonction des valeurs réelles.

In [9]:
plt.scatter(Y,model_sparse.predict(X))
plt.xlabel("IMT réelle")
plt.ylabel("IMT prédit")
Out[9]:
Text(0,0.5,'IMT prédit')

QUESTION Tracer les résidus en fonction des valeurs estimées et juger de l'hypothèse d'homoscédasticité. Qu'observe-t-on ?

In [10]:
plt.scatter(model_sparse.predict(X),Y.T-model_sparse.predict(X))
plt.xlabel("IMT prédit")
plt.ylabel("Résidus")
Out[10]:
Text(0,0.5,'Résidus')

QUESTIONS FINALES

  • Lequel des 2 modèles construits est-il le plus précis ?
  • Lequel le clinicien va-t-il privilégier ?
  • Lequel a le plus d'intérêt ?

QUESTION EN PLUS

Si le temps le permet regarder les résidus en fonction de toutes les variables explicatives et dire si' l'éventuelle hétéroscédasticité peut venir de telle ou telle variable.