Régression linéaire

Pour faire une régression linéaire :

• à partir d'une array X d'observations (en ligne) x paramètres (en colonne) et un vecteur y :

  import statsmodels.regression
  mdl = statsmodels.regression.linear_model.OLS(y, X, hasconst = False)
  res = mdl.fit()
      

  mais par défaut, pas d'ajout de constante (intercept). Si on veut en rajouter une, il faut faire avant la régression : import statsmodels.api; X = statsmodels.api.add_constant(X)
• modèle linéaire avec ordonnée à l'origine (intercept) à partir d'un dataframe pandas (qui a ici au moins les colonnes x1, x2 et y) :

  import pandas
  import numpy
  import statsmodels.formula.api
  df = pandas.DataFrame({'x1': [2, 6, 7, 8, 6, 2], 'x2': [4, 2, 9, 1, 7, 2]})
  df['y'] = df['x1'] * 2 + df['x2'] * 5 + 0.2 * numpy.random.randn(len(df)) + 3
  model = statsmodels.formula.api.ols('y ~ x1 + x2', data = df)
  result = model.fit()
      

• ici, une constante (intercept) est automatiquement rajoutée.
• si on ne veut pas de constante, il faut utiliser la formule : 'y ~ x1 + x2 - 1'
• on peut aussi faire (équivalent) : from statsmodels import regression; model = statsmodels.regression.linear_model.OLS.from_formula('y ~ x1 + x2', data = df)
• result est de type statsmodels.regression.linear_model.RegressionResultsWrapper
• pour avoir les résultats sous forme textuelle, faire result.summary().
• valeurs dans les résultats :
  • result.model : le modèle (de type statsmodels.regression.linear_model.OLS)
  • result.model.formula : la formule du modèle (string).
  • result.model.data.frame : les données d'origine pour le modèle.
  • result.model.exog : la variable à prédire
  • result.model.endog : les variables prédictives sous forme de matrice et incluant en 1ère colonne que des 1 pour l'intercept.
  • result.nobs : le nombre d'observations.
  • result.params : les paramètres du modèle (intercept en tête). C'est une series avec comme noms : Intercept et les noms de colonnes du dataframe (ici, x1 et x2)
  • result.bse : les erreurs standards sur les coefficients de la regression.
  • result.tvalues : les valeurs de la statistique t (result.params / result.bse)
  • result.fittedvalues : les valeurs prédites.
  • result.resid : les résidus (series).
  • result.pvalues : les p values pour chaque paramètre (series).
  • result.f_pvalue : la p value globale.
  • result.rsquared : le R2
  • result.rsquared_adj : le ajusted R2
  • result.df_model : le nombre de degrés de liberté du modèle (nombre de paramètres - 1)
  • result.df_resid : le nombre de degrés de liberté des résidus (nombre de points - nombre de paramètres du modèle).
  • result.aic : l'AIC
  • result.bic : le BIC
  • result.centered_tss : la variance totale (somme des carrés des écarts à la moyenne)
  • result.ess : la variance expliquée (somme des carrés des différences entre valeurs prédites et moyenne)
  • result.ssr : la variance résiduelle (somme des carrés des résidus). centered_tss = ess + ssr.
  • result.mse_model : ess divisé par le nombre de degrés de liberté des paramètres.
  • result.mse_resid : ssr divisé par le nombre de degrés de liberté des résidus.
  • result.mse_total : centered_tss divisé par somme des degrés de liberté des paramètres et des résidus
  • result.fvalue : la statistique F (mse_model / mse_resid)
  • result.conf_int(0.05) : l'intervalle de confiance sur chacun des paramètres au niveau de confiance 0.05 (dataframe à 2 colonnes pour le min et le max).
• test des outliers :

  studentizedResiduals = result.get_influence().resid_studentized_external
  numpy.argwhere(numpy.abs(studentizedResiduals) > scipy.stats.t.ppf(0.975, df = len(df) - parameterNbr - 2))[:,1]
        

  renvoie les indices correspondant aux outliers. On peut l'appliquer de manière récursive (enlever les outliers, et recommencer le calcul du modèle et recalculer les outliers.
• si on a passé un dataframe, on peut avoir les donnees fittées en utilisant result.predict(df) (même résultats que result.fittedvalues).
• on peut alors prédire les valeurs correspondantes à un nouveau dataframe : print(result.predict(pandas.DataFrame({'x1': [2, 1], 'x2': [4, 1]}))) (le résultat est une series).
• statsmodels.api.stats.linear_rainbow(result) : teste l'hypothèse nulle que la relation est bien linéaire.

On peut avoir un intervalle de confiance des valeurs prédites avec :

• import statsmodels.sandbox.regression.predstd
• (stdError, lower, upper) = statsmodels.sandbox.regression.predstd.wls_prediction_std(result) avec stdError l'erreur standard, lower et upper l'intervalle de confiance (par défaut à 0.05)

Regression linéaire robuste aux valeurs extrèmes (outliers) :

• model = statsmodels.robust.robust_linear_model.RLM.from_formula('y ~ x1 + x2', data = df)
• puis, result = model.fit() et l'utilisation de result comme avec la regression linéaire.
• on peut changer la norme utilisée : model = statsmodels.robust.robust_linear_model.RLM.from_formula('y ~ x1 + x2', data = df, M = statsmodels.robust.norms.TrimmedMean()) (le défaut est statsmodels.robust.norms.HuberT(), mais la trimmed mean est souvent utilisée).

df.corr() : permet d'avoir la matrice de corrélation, ce qui donne les variables fortement corrélées, dont il faut éliminer une partie pour ne garder que les variables non corrélées (sinon, regression est instable).

Détermination des multicolinéarités par le VIF (Variance Inflation Factor) : on peut pour cela utiliser la fonction suivante :

import statsmodels.stats.outliers_influence
df = pandas.DataFrame({'x1': x1, 'x2': x2, 'x3': x3, 'y': y})
for i in range(len(['x1', 'x2', 'x3']))])
  print(statsmodels.stats.outliers_influence.variance_inflation_factor(df.loc[:, ['x1', 'x2', 'x3']].values, i))
  

il faut alors éliminer une par une les variables qui donnent une valeur supérieure à 5 (en commençant par la plus grande, puis on refait tourner, etc ...).