ANOVA

Import des outils / jeu de données

import statistics

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import prince
import scipy.stats as stats
import seaborn as sns
import statsmodels.api as sm
from scipy.stats import bartlett, shapiro
from statsmodels.formula.api import ols

from src.config import data_folder
from src.constants import var_categoriques, var_numeriques
from src.utils import init_notebook

init_notebook()

df = pd.read_csv(
    f"{data_folder}/data-cleaned-feature-engineering.csv",
    sep=",",
    index_col="ID",
    parse_dates=True,
)

df_transforme = pd.read_csv(
    f"{data_folder}/data-transformed.csv",
    sep=",",
    index_col="ID",
    parse_dates=True,
)

Variables globales

ANOVA

Problématique

Nous allons tester l’indépendance entre la variable NumStorePurchaes (catégorique, désignant le nombre d’achat en magasin sur un mois pour un client donné) et Income (quantitative continue).
Nous souhaitons répondre à la question : “le revenu influence-t-il effectivement la propension d’un client à acheter plus dans notre magasin ?”

Les boites à moustaches superposées ci-dessous nous donne l’intuition que la réponse est OUI.
Nous allons le vérifier par une ANOVA.

sns.boxplot(df, x="NumStorePurchases", y="Income")

<Axes: xlabel='NumStorePurchases', ylabel='Income'>

On cherche à déterminer si les moyennes des groupes sont significativement différentes.
L’ANOVA fait le rapport des variances interclasse et intraclasse.

Somme des Carrés des Ecarts (SCE) interclasse

$$ SCEinterclasse = \sum\limits_{k=1}^{N}(\bar{Y_k}\ - \bar{Y})^2 * u_k $$

où $ \bar{Y_k} $ est la valeur de la moyenne du groupe k, $ \bar{Y} $ est la moyenne de la population totale et $ u_k $ est le poids du groupe k et N est le nombre de groupe.

SCE intraclasse

$$ SCEintraclasse = \sum\limits_{k=1}^{N} \sum\limits_{i=1}^{n_k} (X_i^k - \bar{Y_k})^2 $$

où $n_k$ est le nombre d’individus dans le groupe k

On a alors $ SCEtotale = SCEintraclasse + SCEinterclasse $

On pose donc :

$H_0$ : Les moyennes de chaque groupe sont égales si la p-value $> 5%$ $H_1$ : Les moyennes de chaque groupe ne sont pas toutes égales si la p-value $< 5%$

Hypothèses à vérifier

1. la variable d’intérêt est qualitative, la variable explicative est quantitative

2. l’indépendance entre les échantillons de chaque groupe

3. l’égalité des variances entre les groupes (homoscédasticité)

4. la normalité des résidus avec un test de Shapiro.

Indépendance

L’indépendance se vérifie en pratique dans la façon dont sont mesurées les données. Ici nos clients ne se connaissent pas, ne sont pas liés entre eux, et agissent indépendemment les uns des autres.

Egalité des variances

df.groupby("NumStorePurchases")["Income"].std()

NumStorePurchases
0     69462.809107
1     77664.598861
2     10728.326529
3     11467.043327
4     17431.118586
5     16340.719749
6     14434.521840
7     13430.040677
8     13379.455717
9     12497.483886
10    11031.807336
11    11274.439006
12    10936.604548
13    11594.032675
Name: Income, dtype: float64

Nous allons effectuer un test de Levene pour vérifier l’égalité statistique des variances.

$H_0$ : Les variances de chaque groupe sont égales si p-value $> 5%$ $H_1$ : Les variances de chaque groupe ne sont pas toutes égales $< 5%$

group_by = [
    df[df["NumStorePurchases"] == i]["Income"]
    for i in range(max(df["NumStorePurchases"]) + 1)
]

stats.levene(
    group_by[0],
    group_by[1],
    group_by[2],
    group_by[3],
    group_by[4],
    group_by[5],
    group_by[6],
    group_by[7],
    group_by[8],
    group_by[9],
    group_by[10],
    group_by[11],
    group_by[12],
    group_by[13],
)

LeveneResult(statistic=41.001633926324985, pvalue=3.3516207987564547e-93)

Notre p-value est largement inférieure à 5%, donc les variances ne sont pas toutes égales.
Nous allons voir ce problème de plus près.

%%capture

group_std = pd.DataFrame(columns=["std"])

for i in range(len(group_by)):
    group_std = group_std.append(
        pd.DataFrame({"std": [statistics.stdev(group_by[i])]}, index=[i])
    )

plt.plot(group_std, "bo")
plt.title("Ecart type des revenus en fonction du nombre d'achat en magasin")

Text(0.5, 1.0, "Ecart type des revenus en fonction du nombre d'achat en magasin")

On constate que les variances sont en effet dispersées, notamment à cause des groupes aux nombres d’achats 0 et 1.
Vérifions que la taille de chacun des groupes est suffisamment grande pour être significatif.

for i in range(len(group_by)):
    print(len(group_by[i]))

print(
    "Les groupes aux nombres d'achats 0 et 1 sont trop petits pour être significatifs, on les retire de l'analyse."
)

15
7
206
447
295
197
161
132
136
99
118
71
93
75
Les groupes aux nombres d'achats 0 et 1 sont trop petits pour être significatifs, on les retire de l'analyse.

del group_by[0:2]
group_std = group_std.drop(labels=[0, 1], axis="index")

plt.plot(group_std, "bo")
plt.title("Ecart type des revenus en fonction du nombre d'achats - données nettoyées")

Text(0.5, 1.0, "Ecart type des revenus en fonction du nombre d'achats - données nettoyées")

On créer un dataframe pour l’ANOVA où l’on conserve uniquement les données qui nous intéressent.

df_anova = df[df["NumStorePurchases"] > 1][["NumStorePurchases", "Income"]]
print(df_anova)

       NumStorePurchases   Income
ID                               
5524                   4  58138.0
2174                   2  46344.0
4141                  10  71613.0
6182                   4  26646.0
5324                   6  58293.0
...                  ...      ...
8080                   3  26816.0
10870                  4  61223.0
7270                  13  56981.0
8235                  10  69245.0
9405                   4  52869.0

[2030 rows x 2 columns]

Normalité des distributions des groupes

group_by_reshape = np.reshape(group_by, (3, 4))

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[15], line 1
----> 1 group_by_reshape = np.reshape(group_by, (3, 4))

File ~/.cache/pypoetry/virtualenvs/customer-base-analysis-F-W2gxNr-py3.10/lib/python3.10/site-packages/numpy/core/fromnumeric.py:285, in reshape(a, newshape, order)
    200 @array_function_dispatch(_reshape_dispatcher)
    201 def reshape(a, newshape, order='C'):
    202     """
    203     Gives a new shape to an array without changing its data.
    204 
   (...)
    283            [5, 6]])
    284     """
--> 285     return _wrapfunc(a, 'reshape', newshape, order=order)

File ~/.cache/pypoetry/virtualenvs/customer-base-analysis-F-W2gxNr-py3.10/lib/python3.10/site-packages/numpy/core/fromnumeric.py:56, in _wrapfunc(obj, method, *args, **kwds)
     54 bound = getattr(obj, method, None)
     55 if bound is None:
---> 56     return _wrapit(obj, method, *args, **kwds)
     58 try:
     59     return bound(*args, **kwds)

File ~/.cache/pypoetry/virtualenvs/customer-base-analysis-F-W2gxNr-py3.10/lib/python3.10/site-packages/numpy/core/fromnumeric.py:45, in _wrapit(obj, method, *args, **kwds)
     43 except AttributeError:
     44     wrap = None
---> 45 result = getattr(asarray(obj), method)(*args, **kwds)
     46 if wrap:
     47     if not isinstance(result, mu.ndarray):

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (12,) + inhomogeneous part.

fig, axs = plt.subplots(3, 4, figsize=(20, 15))

for i in range(3):
    for j in range(4):
        sns.histplot(group_by_reshape[i, j], kde=True, ax=axs[i, j])

Nous utilisons le test de Shapiro-Wilk pour vérifier la normalité des résidus du modèle linéaire OLS.

H0 : Les résidus suivent une loi normale si p-value > 5%
H1 : Les résidus ne suivent pas une loi normale si p-value < 5%

model = ols("NumStorePurchases ~ Income", data=df).fit()
shapiro(model.resid)

La p-value du test de Shapiro-Wilk est inférieure à 5%.
Pour autant, le QQ plot ci-dessous nous indique que nos résidus suivent relativement bien une loi normale.

plt.figure(figsize=(10, 10))
stats.probplot(model.resid, plot=plt, rvalue=True)
plt.title("Probability plot of model's residuals")

Bilan des hypothèses

Hypothèse	Résultat	Remarque
Indépendance	✅	Vérifiée par la nature des données
Egalité des variances	✅	Nécessité de supprimer les groupes à 0 et 1 achat (petit effectif)
Normalité des distributions	✅	Test non concluant mais QQ plot concluant

Test d’ANOVA

H0 : Les moyennes de chaque groupe sont égales si p-value > 5% H1 : Les moyennes de chaque groupe ne sont pas toutes égales < 5%

table = sm.stats.anova_lm(model)
print(table)

Nous obtenons une p-value, pour l’hypothèse nulle “les moyennes des groupes sont égales”, très proche de 0. Nous pouvons donc affirmer que le nombre d’achat en magasin dépend bien du niveau de revenu.

plt.scatter(df_anova["NumStorePurchases"], df_anova["Income"])
plt.title

## TODO: afficher les moyennes
## plt.scatter(df_anova["NumStorePurchases"].unique(), df_anova.groupby(by="NumStorePurchases").mean()["Income"], color="red")