Retrouvez l'énoncé et les détails de la mission traitée en cliquant ici. ¶
Sommaire ¶
I. Nettoyage des données
Mission 3: Génération des individus
I. Coefficients d'élasticité
I. Sélection diversifiée des pays à analyser
I. Préliminaires des analyses à venir
Nous importons dans un premier temps l'ensemble des librairies nécessaires pour l'ensemble du projet. ¶
import pandas as pd ; import numpy as np ; import random ; from IPython.display import Image
import matplotlib.pyplot as plt ; from matplotlib.pylab import arange, plot ; from matplotlib.pyplot import figure
from sklearn import preprocessing ; from sklearn.preprocessing import StandardScaler ; sc = StandardScaler()
import scipy ; import scipy.stats as st ; from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
from scipy.stats import t, shapiro, ks_2samp ; from collections import Counter
import statsmodels.formula.api as smf ; import statsmodels.api as sm
from statsmodels.formula.api import ols ; from statsmodels.stats.outliers_influence import variance_inflation_factor
# Ci-dessous sont importées les fonctions données par l'énoncé du projet 7 pour la mission 3, enregistrés dans un dossier .py:
from codes_mission_3 import* ; from functions import*
# La fonction dim() nous permettra de définir à chaque fois les dimensions des figures que l'on affichera:
def dim(x, y): return(plt.figure(figsize=(x, y)))
# Outils de génération des cartes du monde colorisées:
import cartopy ; import cartopy.io.shapereader as shpreader ; import cartopy.crs as ccrs
shpfilename = shpreader.natural_earth(resolution="110m", category="cultural", name="admin_0_countries")
reader = shpreader.Reader(shpfilename) ; from world_cartography import *
from IPython.core.display import HTML ; import warnings ; warnings.filterwarnings('ignore')
HTML("""<style> .output_png {display: table-cell; text-align: center; vertical-align: middle} </style>""")
Invoquons notre DataFrame de base - que nous appellerons "df" - issu de la World Income Distribution qui contient les variables suivantes:
- [country_code]: Le code ISO d'un pays donné.
- [year_survey]: L'année où ont été relevées les données.
- [quantile]: Pour chaque pays, on dispose de 100 quantiles, chacun représentant 1% de la population. Ils sont numérotés de 1 à 100 où quantile numéro 1 représente la moyenne des revenus la plus faible et le quantile numéro 100 la moyenne la plus élevée. Cette colonne indique alors le numéro du quantile en question.
- [nb_quantile]: Renvoie le nombre de quantiles pour un pays donné. En l'occurence, ce nombre est toujours égal à 100.
- [income]: La moyenne des revenus pour un quantile donné dans un pays donné au cours de l'année en question, converti en dollards.
- [gdpppp]: Le PIB par habitant à parité de pouvoir d'achats.
df=pd.read_csv("data/bank/data_p7.csv").rename(columns={"country": "country_code", "year_survey": "year"}) ; df
Examinons dans un premiers temps le type de chaque colonne de notre base.
print(df.info())
On constate que les colonnes "income" et "gdpppp" représentent des nombres mais sont interprétés comme des objets. En effet, certaines de ces données sont sous forme de string. Nous convertirons donc ces deux colonnes en nombres flottants, et pour une meilleure visibilité dans toute la suite du sujet, nous arrondirons ces valeurs au cents (centime américain) près.
df["income"]=[float(str(i).replace(",",".")) for i in df.income]
df["gdpppp"]=[float(str(i).replace(",",".")) for i in df.gdpppp]
df[["income", "gdpppp"]].info()
Constatons que le PIB par habitant des Îles Fidji paraissait anormalement élevé par rapport aux autres pays.
print("PIB/hbts de la France:", df.loc[df["country_code"]=="FRA", "gdpppp"].unique()[0])
print("PIB/hbts du Luxembourg:", df.loc[df["country_code"]=="LUX", "gdpppp"].unique()[0])
print("PIB/hbts des Îles Fidji:", df.loc[df["country_code"]=="FJI", "gdpppp"].unique()[0])
À titre comparatif avec la France et le Luxembourg, tous deux présentants cet indice parmi les plus élevés, notre base de données indique que celui des Îles Fidji serait environ 50 fois plus grand que ces deux derniers. D'après la World Bank - source https://data.worldbank.org/indicator/NY.GDP.PCAP.PP.CD?locations=FJ, cet indice s'élèverait en réalité à environ 7777.7 dollars PPA. Nous le remplacerons donc dans notre base de données.
df.loc[df["country_code"]=="FJI", "gdpppp"]=7777.7
print("PIB/hbts corrigé des Îles Fidji:", df.loc[df["country_code"]=="FJI", "gdpppp"].unique()[0])
Intéressons-nous au nombre de pays que contient df et à ses dimensions.
print("Nombre de lignes de df: %d" %len(df))
print("Nombre de pays contenu dans df: %s" %len(df["country_code"].unique()))
Chaque pays est supposé être représenté par 100 quantiles (centiles). Donc les lignes de notre tableau sont censée être exactement au nombre de 100 fois le nombre de pays présents. Les valeurs données ci-dessus nous indiquent donc qu'un centile est absent pour un pays donné. Nous allons donc l'identifier ce pays ainsi que son quantile manquant et ajouter la ligne manquante à df.
for cc in df["country_code"].unique(): # Identification du pays concerné
if len(df.loc[df["country_code"]==cc, "quantile"])<100:
country=cc
print("Code ISO du pays concerné:", cc)
for i in range(1, 101): # Identification du quantile manquant
if i not in df.loc[df["country_code"]==country,"quantile"].values: print("Il lui manque le %d ème quantile!" %i)
Il s'agit donc de la Lituanie, et le 41ème quantile est donc absent. Nous allons donc nous baser sur les revenus moyens des centiles 40 et 42 de la Lituanie puis prendre leur moyenne pour estimer le revenu moyen du centile 41. Puis nous ajoutons cette nouvelle ligne à df.
filtre1, filtre2 = (df["country_code"]=="LTU") & (df["quantile"]==40), (df["country_code"]=="LTU") & (df["quantile"]==42)
rm = (df.loc[filtre1, "income"].values[0] + df.loc[filtre2, "income"].values[0])/2
df.loc[len(df)]=["LTU", 2008, 41, 100, rm, 17571.0]
# On peut observer que le nombre de lignes de df est bien passé de 11599 à 16000:
df=df.sort_values(by=["country_code", "quantile"]).reset_index(drop=True) ; df
Repérons maintenant les éventuelles colonnes présentant au moins une valeur manquante.
df.isnull().sum()
On observe une 200 valeurs manquantes dans la colonne "gdpppp". Repérons alors les pays concernés.
[i for i in df.loc[df["gdpppp"].isnull()==True]["country_code"].unique()]
À titre comparatif avec la France et le Luxembourg, tous deux présentants cet indice parmi les plus élevés, notre base de données indique que celui des Îles Fidji serait environ 50 fois plus grand que ces deux derniers. D'après la World Bank, cet indice s'élèverait en réalité à environ 7777.7 dollars PPA. Nous le remplacerons donc dans notre base de données.
df.loc[df["country_code"]=="XKX"]=df.loc[df["country_code"]=="XKX"].fillna(7249)
df.loc[df["country_code"]=="PSE"]=df.loc[df["country_code"]=="PSE"].fillna(3816.20)
# On peut vérifier que l'on en a terminé avec les valeurs manquantes dans le tableau de base:
df.isnull().sum()
Récupérons dans un premier temps le nombre d'habitants de chaque pays présents dans notre DataFrame pour l'année en question. Ici, nous avons récupérer sur la faostats la liste des populations de chaque pays de 2004 à 2011 (années où s'étendent nos données de base). Nous appellerons "pop" ce DataFrame. Nous en profiterons pour récupérer également le nom des pays indiqués jusqu'ici par leur code ISO.
pop=pd.read_csv("data/bank/gini/population.csv").drop(
columns=["Code Domaine", "Domaine", "Code Élément", "Élément", "Produit", "Code Produit", "Symbole",
"Description du Symbole", "Code année", "Unité"]).rename(columns={"Code zone": "country_code",
"Zone": "country", "Année": "year",
"Valeur": "population"})
pop["population"]=1000*pop.population
pop.loc[pop["country"]=="Royaume-Uni de Grande-Bretagne et d'Irlande du Nord", "country"]="Royaume-Uni"
df=pd.merge(df, pop, how="left").reindex(columns=["country_code", "country", "year", "quantile", "nb_quantiles", "income",
"gdpppp", "population"]) ; df.tail(4)
On récupère également les données manquantes à la liste trouvée sur la faostats, à savoir les populations et noms des pays ci-dessous, d'après le site: http://datatopics.worldbank.org/world-development-indicators/.
df.loc[df["country_code"]=="BTN", "country"]="Bhoutan"
df.loc[df["country_code"]=="XKX", "country"]="Kosovo"
df.loc[df["country_code"]=="SDN", "country"]="Soudan"
df.loc[df["country_code"]=="SYR", "country"]="Syrie"
df.loc[df["country_code"]=="PSE", "country"]="Palestine"
df.loc[df["country_code"]=="COD", "country"]="République du Congo"
df.loc[df["country_code"]=="BTN", "population"]=664876
df.loc[df["country_code"]=="XKX", "population"]=1747000
df.loc[df["country_code"]=="SDN", "population"]=33780000
df.loc[df["country_code"]=="SYR", "population"]=17830000
df.loc[df["country_code"]=="PSE", "population"]=3689000
df.loc[df["country_code"]=="COD", "population"]=60410000 ; df
Nous allons rapportés la population de chaque pays présents dans la base à 100 habitants, chacun représenté par un centile. Et nous considérerons le revenu moyen attribué à un centile dans la colonne "income" comme le revenu de ce centile pour calculer l'indice de Gini de chaque pays présents dans "df", puis nous ajouterons ces indices à celui-ci.
# On crée ici la liste P de chaque pays présent dans notre base:
P = df.country_code.unique()
# On crée maintenant la liste G qui contiendra l'indice de Gini de chaque pays présent dans P:
G=[]
for i in P:
revenus = df.loc[df["country_code"]==i].income.values
lorenz = np.append([0], np.cumsum(revenus)/sum(revenus))
gini=2*(0.5-lorenz.sum()/len(revenus))
G=np.append(G, gini)
# Créons le DataFrame contenant chaque pays présents dans notre base accompagné de son indice de Gini:
gini_list=pd.DataFrame({"country_code": P, "gini": G})
# Puis on l'ajoute à "df" par une jointure gauche:
df=pd.merge(df, gini_list, how="left") ; df
• Nous avons alors une liste de 116 pays dont les données sont relevées chacun pour une seule année donnée entre 2004 et 2011 incluses, exceptée 2005.
• La colonne "gdpppp" (Gross Domestic Product based on Purchasing Power Parity) représente le PIB à parité de pouvoir d'achat calculé divisé par le nombre d'habitants d'un pays donné.
• La population de chaque pays est représentée par des centiles, nous avons divisé par 100 chaque population, et regroupé les habitants selon leur revenu. Plus le revenu est élevé, plus l'indice du quantile dans lequel se situe un habitant sera proche de 100, puis nous avons donné le revenu moyen des habitants de chaque quantile dans la colonne "income".
• Échantillonner une population en utilisant des quantiles est, selon moi, une bonne méthode car cela permet de minimiser le besoin en récoltes des données, et nous permet d'éviter les cas particuliers. Cela nous a permis de donner un indice de Gini à chaque pays dans la colonne "gini" en ne considérants que les revenus moyens des centiles donnés.
La plupart des données que nous traitons ayant été relevées en 2008. Considérons alors la population mondiale comme la population mondiale de 2008 + les valeurs ajoutées dans la partie précédente par la World Bank, et calculons le pourcentage de la population mondiale couverte par notre analyse.
XX=["BTN", "XKX", "SDN", "SYR", "PSE", "COD"]
eps=df.loc[df["country_code"].isin(XX)==True, "population"].unique().sum()
per=round(100*sum(df.drop_duplicates(["country_code"]).population)/(sum(pop.loc[pop["year"]==2008].population)+eps), 2)
CARTO="L'analyse couvre donc environ {}% de la population mondiale et les pays ci-dessous".format(per)
ciso=mapzone(df[["country_code", "country"]].drop_duplicates().set_index("country"))
dim(22, 7); ax = plt.axes(projection=ccrs.PlateCarree()) ; ax = carte(ax) ; cartography([ciso], ["crimson"], ax)
plt.title(CARTO, fontsize=24) ; plt.show()
Il nous faut sélectionner entre 5 à 10 pays de notre base de données pour montrer la diversité mondiale en terme de distribution de revenus. Pour les sélectionner de manière la plus objective possible, une classification hiérarchique ascendante non-supervisée de ces pays sur df selon les variables "income" et "gdpppp", en tracerons un dendrogramme pour déterminer le nombre de groupes distincts raisonnables à étudier, et sélectionnerons un pays donné par groupe obtenu.
var=df[["country_code", "income", "gdpppp"]]
var=pd.DataFrame(var.groupby("country_code").mean())
X=var.values
std_scaler = preprocessing.StandardScaler().fit(X)
X_scaled = std_scaler.transform(X)
df_scaled=pd.DataFrame(X_scaled, index=var.index, columns=var.columns) ; df_scaled
# Préparation des données pour le clustering:
names = df_scaled.index
# Clustering hiérarchique:
Z = linkage(X, 'ward')
plot_dendrogram(Z, names) ; plt.show()
En analysant le dendrogramme, il nous laisse suggérer que réaliser 8 clusters nous permettrait une bonne répartition des diversités.
Y = df_scaled.values
Z = linkage(Y, 'ward')
# On ajoute la colonne "cluster" à "df" qui indique le cluster attribué à un pays pays donné selon notre classification:
clusters = fcluster(Z, 8, criterion='maxclust')
var["cluster"]=clusters
varc=var[["cluster"]].reset_index()
df=pd.merge(df, varc, how="left") ; df
Nous sélectionnons alors un pays par cluster. On aura pris soin de choisir pour chaque cluster un pays où les indices de Gini de 2004 à 2011 sont pour le mieux donnés par la Banque Mondiale. Nous allons donc construire un DataFrame - que nous appellerons "selection" - contenant la liste des 8 pays choisis, leur indice de gini pour l'année donnée dans notre base et les numéros de cluster auquel ils ont été attribués lors de la classification hierarchique.
# Ci-dessous, la liste S des codes ISO des 8 pays sélectionnés:
S=["FRA", "ESP", "ISL", "LUX", "ARM", "RUS", "UKR", "ARG"]
P=[] # La liste P nous affichera les noms des pays en question.
G=[] # La liste G nous affichera, elle, l'indice de gini attribué à chaque pays.
C=[] # Enfin, La liste C nous indiquera le cluster auquel appartient chaque pays.
for i in S:
filtre = df["country_code"]==i
country = df.loc[filtre, "country"].unique()[0]
clust = df.loc[filtre, "cluster"].unique()[0]
ind_gini = df.loc[filtre, "gini"].unique()[0]
P=np.append(P, country)
C=np.append(C, clust)
G=np.append(G, ind_gini)
# Résumons alors les pays sélectionnés et leurs données dans un DataFrame:
selection=pd.DataFrame({"country_code": S, "country": P, "gini": G, "cluster": range(1, 9)}) ; selection
À présent, nous allons représenter graphiquement le revenu moyen en fonction de chacune des classes de revenus pour les 8 pays sélectionnés ci-dessus pour montrer la diversité des cas.
figure(num=None, figsize=(8, 6), dpi=80, facecolor="w", edgecolor="k")
for n, i in enumerate(selection["country_code"]):
filtre = df["country_code"]==i
plt.plot(df.loc[filtre, "quantile"], df.loc[filtre, "income"], label="Cluster %d" %(n+1))
plt.title("Diversité des classes de revenus dans le monde", fontsize=16) ; plt.legend(fontsize=10)
plt.xlabel("Quantiles", fontsize=14); plt.ylabel("Revenus Moyens",fontsize=14); plt.yscale("symlog",linthreshx=1); plt.show()
On constate tout d'abord que les étendues des revenus pour chaque pays sont relativement disjointes. Par exemple, le plus haut revenu moyen du cluster 5 est inférieur au plus bas revenu du cluster 4. Cependant, les courbes semblent relativement superposables. On peut donc estimer que les inégalités de revenus entre les différents pays sont plus marquées au niveau international qu'au niveau national.
Nous allons donc afficher la courbe de Lorenz pour chaque pays sélectionné. Nous indiquerons chaque cluster dans le titre de chaque graphique.
figure(num=None, figsize=(8, 5), dpi=80, facecolor="w", edgecolor="k")
for p,i in enumerate(selection["country"]):
data=df.loc[df["country"]==i]
revenus = data.income.values
n = len(revenus)
lorenz = np.append([0], np.cumsum(revenus)/sum(revenus))
gini=selection["gini"][p]
plt.axes().axis("equal")
plt.plot(np.linspace(0-1/n, 1+1/n, n+1),
lorenz, label=" %s - Cluster %d - Indice de Gini: %f" % (i, p+1, gini))
X = arange(0,1,0.01) ; Y=X
plot(X, Y)
plt.title("Courbes de Lorenz", fontsize=20) ; plt.legend(fontsize=9.5) ; plt.tight_layout(); plt.show()
Les indices de Gini calculés se situent entre 0,24 et 0,44. On constatera tout de même de bonnes inégalités en terme de distribution de revenus dans chaque pays. Tous éloignés de 0, on en déduit que les revenus ne sont pas équitablement répartis dans aucun état.
Le DataFrame "gin" nous indique les indices de gini que la Banque Mondiale a pu estimer pour chaque pays au fil des ans. Nous ne conserverons que la période donnée dans notre base, à savoir les années 2004 à 2011, et uniquement les pays présents dans notre sélection.
gin=pd.read_csv("data/bank/gini/api.csv", skiprows=[0,1,2])
# On supprime les colonnes inutiles et on renomme les deux premières:
gin=gin.drop(columns="Unnamed: 64").rename(columns={"Country Code": "country_code"}).set_index("country_code")
# On ne garde que les années 2004 à 2011:
gin=gin.iloc[:, 47:]
gin=gin.iloc[:, 0:8]
gin=gin.reset_index()
# On ne conserve que les pays présents dans notre sélection:
gin=gin.loc[gin["country_code"].isin(selection["country_code"])==True]
# On le réordonne comme le DataFrame "selection":
sel=selection[["country_code"]].reset_index()
gin=pd.merge(gin, sel).sort_values(by="index").set_index("index").reset_index().drop(columns=["index"])
# Le DataFrame renvoie les indices de Gini en pourcentage, nous les divisons donc par 100:
gin.iloc[:, 1:]=gin.iloc[:, 1:]/100 ; gin
Nous allons donc maintenant représenter graphiquement pour ces 8 pays l'évolution de l'indice de Gini de 2004 à 2011.
figure(num=None, figsize=(10, 5), dpi=80, facecolor="w", edgecolor="k")
for i in gin.index:
ans=pd.DataFrame(gin.iloc[i, 1:9]).reset_index()
ans.columns=["year", "gini"]
pays=selection["country"][i]
groupe=selection["cluster"][i]
plt.plot(ans.year, ans.gini, label=("%s - Cluster %d" %(pays, groupe)))
plt.legend(loc="upper left", fontsize=12, bbox_to_anchor=(1, 1)) ; plt.ylabel("Indice de Gini", fontsize=20)
plt.title("Évolution de l'indice de Gini au fil des ans", fontsize=25) ; plt.show()
On notera généralement une tendance à la baisse des indices de gini proches de 0,5 au fil des ans, tandis que les indices de gini les plus bas ont présentent une légère tendance à la hausse ou bien semblent osciller autour d'une certaine constante.
Pour classer les pays choisis selon leur indice de Gini, il nous suffit de reprendre le DataFrame "selection" et de l'ordonner selon la variable "gini". On crée un DataFrame qui renvoie cette liste ordonnée que nous appellerons "classement".
classement=selection[["country", "gini", "cluster"]].sort_values(by="gini").reset_index().drop(columns="index") ; classement
Nous allons maintenant calculer la moyenne de ces indices de gini, ainsi que la moyenne des 4 pays présentant l'indice le plus bas, puis la moyenne des 4 présentant les plus haut.
low_classement , top_classement = classement.head(4) , classement.tail(4)
print("Moyenne des 8 pays:", classement["gini"].mean())
print("Moyenne des 4 pays présentant les indices de gini les plus bas:", low_classement["gini"].mean())
print("Moyenne des 4 pays présentant les indices de gini le plus hauts:", top_classement["gini"].mean())
La France semble présenter une distribution de revenus légèrement moins équitable que la plupart des autres pays du monde, mais son indice de gini n'est pas très éloigné de la moyenne générale.
L'objectif de cette première partie de la mission 3 est de récupérer les coefficients d'élasticité propres à chaque pays présents dans notre base de données et de les ajouter à df. Le DataFrame suivant - qu'on appellera "gdi" - issu de la World Bank, nous indiquera ce coefficient pour une première partie de nos pays.
gdi=pd.read_csv("data/bank/gdim.csv")[["wbcode", "year", "IGEincome"]].rename(columns={"wbcode": "country_code",
"IGEincome": "elasticity"})
gdi=gdi.sort_values(by=["country_code", "year"]) ; gdi.sort_values(by="elasticity").drop_duplicates()
Nous allons nous débarrasser des doublons de ce DataFrame, et conserver l'élasticité la plus récente pour les pays où elle est indiquée. Puis nous ajouterons une nouvelle colonne "elasticity" à df qui nous renverra déjà l'élasticité des pays où elle est indiquée dans gdi, et renverra un NaN pour ceux où elle ne l'est pas pour le moment.
# Cration de la table "ela" qui représentera "gdi" sans les doublons avec, si elle est indiquée, l'élasticité la plus récente
# pour un pays donné:
ela=pd.DataFrame({"country_code": [], "year": [], "elasticity": []})
for i in gdi["country_code"].unique():
loca=gdi["country_code"]==i
if len(gdi.loc[loca].dropna())>0:
ok=gdi.loc[loca].dropna().drop_duplicates()
ela=pd.concat([ela, ok])
else:
ko=gdi.loc[loca].drop_duplicates(subset=["elasticity"])
ela=pd.concat([ela, ko])
ela=ela.drop(columns="year")
# Ajout de la colonne "elasticity" à df:
df=pd.merge(df, ela, how="left") ; df
On observe que nous n'avons pas encore l'ensemble des indices d'élasticité de tous les pays observés, ce qui fera l'objet de la sous-partie suivante.
df["elasticity"].isnull().sum()
Nous allons donc maintenant repérer les pays de notre base de données pour lesquels l'élasticité n'a pas été relevée par gdi. Le DataFrame suivant, "estim_elasticity", propose une estimation de l'élasticité de chaque pays selon leur région continentale. Nous identifierons alors la liste des pays de df où l'élasticité n'est pas encore indiquée ainsi que la région continentale à laquelle appartient chacun de ces pays, puis nous leur attribuerons l'élasticité moyenne indiquée dans la colonne "Base case" du DataFrame "estim_eslasticity" définit ci-dessous.
estim_elasticity=pd.read_csv("data/bank/elasticity.csv").set_index("Unnamed: 0").rename_axis(None) ; estim_elasticity
Identifions les pays de df dont l'élasticité reste à déterminer.
B=[i for i in df.loc[df["elasticity"].isnull()==True]["country"].unique()] ; print(B)
Créons maintenant de la liste des élasticités attribués dans "estim_elasticity" aux régions territoriales auxquelles appartiennent chacun de ces pays.
[a, b, c, d, e] = estim_elasticity["Base case"]
K=[e, d, d, e, b, d, e, e,
e, d, e, b, c, b, e, b, d,
d, d, a, d, d, b, d, e,
d, b, b, e, b, e, e, e, e,
d, b, e, b, e, e, b, e, d, d,
d, d, b, b, e, e, d, d, e]
# Ajout des valeurs manquantes dans la colonne "elasticity" de "df":
CC=[i for i in df.loc[df["elasticity"].isnull()==True]["country_code"].unique()]
ajout=pd.DataFrame({"country_code": CC, "elasticity": K}) ; ajout.head(10)
Nous pouvons donc enfin ajouter les indices d'élasticités manquants dans df afin de compléter le DataFrame.
for i in df.index:
a=df["country_code"][i]
if (a in list(ajout["country_code"]))==True:
df["elasticity"][i]=list(ajout.loc[ajout["country_code"]==a]["elasticity"])[0]
df
Nous allons générer une population dans le cadre de la construction d'un tableau des distributions conditionnelles. Dans la cellule [1] des importations de librairies ont été, entre autres, importées les fonctions données par l'énoncé d'Openclassrooms dans le cadre de la mission 3. Commençons par générer les revenus d'une population et ceux de leurs parents. Nous prendrons pour cette partie un indice d'élasticité de 0,4.
# Nombre de quantiles:
nb_quantiles = 100
rg = range(1, nb_quantiles+1)
# Nombre d'individus générés et coefficient d'élasticité choisi:
n, pj = 1000*nb_quantiles, 0.9
# Génération aléatoire des revenus des enfants et de leurs parents:
y_child, y_parent = generate_incomes(n, pj)
y_child, y_parent = pd.Series(y_child), pd.Series(y_parent)
y_child.head()
On détermine maintenant la classes de revenu de chaque enfant et de ses parents.
c_i_child = quantiles(y_child, 100) # Classes de revenu des enfants.
c_i_parent = quantiles(y_parent, 100) # Classes de revenu des parents.
# Aperçut des classes de revenus:
sample=compute_quantiles(y_child, y_parent, nb_quantiles) ; sample
Établissons maintenant le tableau de distribution pour un indice d'élasticité de 0,4. Le DataFrame suivant comprendra en index la classe de revenu d'un enfant, en colonne la classes de revenu d'un parent et dans chaque cellule, la probabilité qu'un parent appartienne à la classe indiquée en index sachant que son enfant appartiennent à la classe de revenu indiquée en colonne et que l'élasticité est de 0,4.
cd = conditional_distributions(sample, nb_quantiles)
dis_rev = pd.DataFrame(cd, columns=range(1, 101), index=range(1, 101)) ; dis_rev
À titre d'exemple, calculons la probabilité qu'un parent appartienne à la classe de revenus 47 sachant que son enfant appartient à la classe 50 et que l'indice d'élasticité est de 0,4.
print("\nP(c_i_parent = {} | c_i_child = {}, pj = {}) = {}".format(47, 50, 0.4, dis_rev.loc[50, 47]))
Cette seconde sous-partie est une anticipation de la sous-partie suivante qui conclura la mission 3. Nous y présentons notamment la fonction qui nous permettra plus tard de construire les classes de revenu des parents des individus de la World Income Distribution. On généralise la sous-partie précédante en créant la fonction "tableau_distributivite" qui prendra en argument un indice d'élasticité pj, et qui nous renverra le tableau des distributions de revenus selon pj pour une population générée aléatoirement selon une loi normale, et un échantillonnage de 100 quantiles.
def tableau_distributivite(pj):
nb_quantiles=100
n=500*nb_quantiles
rg=range(1, nb_quantiles+1)
y_child, y_parent = generate_incomes(n, pj)
y_child, y_parent = pd.Series(y_child), pd.Series(y_parent)
c_i_child=quantiles(y_child, nb_quantiles)
c_i_parent=quantiles(y_parent, nb_quantiles)
sample=compute_quantiles(y_child, y_parent, nb_quantiles)
cd=conditional_distributions(sample, nb_quantiles)
data=pd.DataFrame(cd, columns=rg, index=rg)
return(data)
Pour une question pratique, nous utiliserons aussi la fonction "effectif" qui prend une série, et renvoie une nouvelle série qui comprend i fois le nombre k si l'index de k est i.
def effectif(S):
W=[]
for i in S.index:
nb=int(S[i])
W+=[i]*nb
return(W)
Retirons dans un premier temps les colonnes de df qui ne nous servirons plus jusqu'à la fin du projet, et renommons, conformément à l'énoncé, la colonne "gdpppp" en "mj", que nous considérerons dorénavant comme le revenu moyen du pays, ainsi que la colonne "income" en "y_child", qui représentera maintenant le revenu des individus, et la colonne "quantile" en "c_i_child" qui représentera la classe de revenus de ces individus.
df=df.drop(columns=["country", "nb_quantiles", "year", "cluster"]).rename(columns={"gdpppp": "mj", "income": "y_child",
"quantile": "c_i_child"}) ; df
Reprenons maintenant df, et tirons-en le DataFrame "individus" qui représentera notre échantillon, en multipliant chaque ligne de df par 500.
individus=pd.concat([df]*500)
individus=individus.sort_values(["country_code", "c_i_child"]).reset_index(drop=True)
individus=individus.reindex(columns=["country_code", "c_i_child", "y_child", "mj", "gini", "elasticity", "population"])
individus
Conformément à ce qui nous est demandé dans l'énoncé de la mission 3, nous allons ajouter au DataFrame "individus" une nouvelle colonne "c_i_parents" représentant la distribution des classes de revenus des parents, dont le taux d'aparitions d'une même classe de revenu pour un parent dont l'enfant présente une classe donné dans un pays donné (donc pour un coefficient d'élasticité donné) sera proportionnel à la probabilié d'appartenance de cette classe sachant ces derniers paramètres.
df.to_csv("data/bank/transferts/df_before.csv", index=None)
individus.to_csv("data/bank/transferts/individus_before.csv", index=None)
L'exécution de l'instruction pour cet ajout étant trop volumineux pour mon matériel, la partie suivante a donc été réalisée sur Collaboratory. C'est pourquoi ci-dessus ont été enregistrés "df" et "individus", et est joint ci-dessous une capture d'écran du code exécuté sur Colaboratory pour obtenir "individus" doté de sa nouvelle colonne "c_i_parents".
# Génération de la colonne "c_i_parents":
Image(filename='data/bank/algo_capture.png')
Ci-dessous, on rappelle alors le DataFrame "individus" avec la modification exécutée sur Collaboratory, à savoir l'ajout de la colonne "c_i_parents" qui renvoie pour chaque quantile d'un pays donné une liste de 500 classes de revenu des parents, conformément aux conditions demandées dans l'énoncé du projet.
individus=pd.read_csv("data/bank/transferts/individus_after_m3.csv")
df=pd.read_csv("data/bank/transferts/df_before.csv")
individus=individus.drop(columns=["mj", "country"]).rename(columns={"gdpppp": "mj"})
individus=individus.reindex(columns=["country_code", "c_i_child", "y_child", "mj", "gini", "c_i_parents",
"elasticity", "population"])
# On peut visualiser la nouvelle colonne "c_i_parents" dans le DataFrame "individus":
individus
Terminons cette mission 3 en supprimant la colonne ["c_i_child] et en ajoutant la colonne ["y_parents"]: Le revenu des parents basé sur le quantile "c_i_parents" du pays donné dans le tableau de base.
# Suppression de la colonne "c_i_child":
individus=individus.drop(columns="c_i_child")
# Ajout de la colonne "y_parents":
income_p=df[["country_code", "c_i_child", "y_child"]]
income_p.columns=["country_code", "c_i_parents", "y_parents"]
individus=pd.merge(individus, income_p)
# Enfin, on réorganise les colonnes de "individus":
individus=individus.reindex(columns=["country_code", "y_child", "y_parents", "c_i_parents", "mj", "gini", "elasticity",
"population"]) ; individus
On notera que ce nouveau DataFrame contient finalement donc bien les 3 variables demandées:
- mj, le revenu moyen du pays j indiqué dans la colonne [mj].
- Gj, l'indice de Gini du pays j indiqué dans la colonne [gini].
- c_i_parent, la classe de revenu des parents d'un individu dans la colonne [c_i_parents].
On rappelle pour la mission 4 que nous avons conservé "df" avant la multiplication par 500 des individus.
Nous réaliserons chaque test statistique au niveau de 5%. On rappelle dans la cellule suivante notre échantillon "individus" décrit à la fin de la mission précédente, et on en profite pour fixer notre niveau de test nt pour toute modélisation qui s'exécutera dans toute la suite du projet.
# Fixation du niveau de test:
nt=0.05
# Chacun des tests sera basé sur notre nouvelle base de données "individus":
individus
Dans le cadre des modélisations à venir, nous ajoutons à "individus" les colonnes "log_y_child", "log_y_parents" et "log_mj", respectivements les données de "y_child", "y_parents" et "mj" exprimées en logarithme. Nous ajoutons également "log_y_child" à df car cela va nous être utile dans la partie suivante.
df["log_y_child"]=np.log(df["y_child"])
individus["log_y_child"]=np.log(individus["y_child"])
individus["log_mj"]=np.log(individus["mj"])
individus["log_y_parents"]=np.log(individus["y_parents"]) ; individus
En vue des analyses répétées que nous exécuterons sur les résidus des modélisations à venir, nous créons ici les fonctions décrites ci-dessous qui prendront toutes en entrée une modélisation donnée pour renvoyer un certain résultat sur ses résidus. À noter que les différents tests d'adéquation à une loi gausienne que nous exécuterons testeront donc l'adéquation d'une certaine variable X à la loi normale de paramètres (mean(X), std(X)).
- 1°) "residus_henry": qui nous superposera les résidus d'une modélisation donnée avec la droite de Henry.
- 2°) "nuage_residuel": qui nous renverra le nuage de la variance résiduelle.
- 3°) "analyse_des_residus": qui nous renverra la droite de Henry ET le nuage résiduel dans un subplots horizontal.
- 4°) "test_shapiro": qui nous renverra les résultats d'un test de Shapiro sur les résidus (peu probable qu'on s'en serve vu la taille des échantillons) .
- 5°) "test_residus_KS": qui nous renverra les résultats d'un test de Kolmogorov-Smirnoff sur les résidus.
- 6°) "test_residus_BP": qui nous renverra les résultats d'un test de Breusch-Pagan sur les résidus en vue de l'analyse de leur homoscédasticité.
Pour chaque analyse des résidus d'un modèle de régression, nous observerons son éventuelle adéquation à la loi normale sur la droite de Henry, et son homoscédasticité sur le nuage de variance résiduelle.
# 1°) Droite de Henry:
def residus_henry(MODEL, nom, couleur, ax):
sm.qqplot(MODEL.resid, ax=ax, line="45", fit=True, color=couleur)
ax.set_title("Droite de Henry - Résidus de %s" %nom, fontsize=18)
ax.set_xlabel("Quantiles théoriques" ,fontsize=16), ax.set_ylabel("Quantiles observés", fontsize=16)
# 2°) Nuage de la variance résiduelle:
def nuage_residuel(MODEL, nom, couleur, ax):
ax=plt.plot(MODEL.fittedvalues, MODEL.resid, ".", color=couleur, alpha=0.3)
plt.title("Nuage de la variance résiduelle de %s" %nom, fontsize=18)
plt.xlabel("GWh", fontsize=16), plt.ylabel("Résidus", fontsize=16)
# 3°) Les deux en même temps:
def analyse_des_residus(MODEL, nom, couleur):
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))
residus_henry(MODEL, nom, couleur, ax1)
nuage_residuel(MODEL, nom, couleur, ax2)
# 4°) Test de Shapiro:
def test_shapiro(MODEL):
SH=shapiro(MODEL.resid)
print("TEST DE SHAPIRO:", " statistic =", SH[0], " pvalue =", SH[1])
# 5°) Test de Kolmogorov-Smirnov:
def test_residus_KS(MODEL):
MR=MODEL.resid
KS=ks_2samp(MR,list(np.random.normal(MR.mean(), MR.std(), len(MR))))
print("TEST DE KOLMOGOROV-SMIRNOV:", " statistic =", KS[0], " pvalue =", KS[1])
# 6°) Test de Breusch-Pagan:
def test_residus_BP(MODEL):
BP=sm.stats.diagnostic.het_breuschpagan(MODEL.resid, MODEL.model.exog)
print("TEST DE BREUSCH-PAGAN:", " statistic =", BP[0], " pvalue =", BP[1])
Les analyses de la mission 4 se porteront notamment sur diverses modélisations de données. Pour chaque régression linéaire multiple que nous exécuterons, il nous faudra analyser la multicolinéarité des variables explicatives incluses dans le modèle. La fonction suivante "VIF" nous prendra en argument un modèle de régression linéaire multiple donné et la liste des variables explicatives du dit-modèle pour nous renvoyer un DataFrame indiquant le facteur d'influence de la variance (VIF en anglais) de chacune de ces variables explicatives au cours de la modélisation.
def VIF(RLM, variables):
MEX=RLM.model.exog
return(pd.DataFrame({"VIF de la variable": [variance_inflation_factor(MEX, i) for i in np.arange(1, MEX.shape[1])]},
index=variables))
Nous créons également, ici, la fonction "decomposition_des_variables" qui, pour une régression linéaire multiple donnée, nous renverra directement un DataFrame affichant pour chaque variable explicative le pourcentage d'explication de la variance du revenu des individus.
def decomposition_des_variables(RLM):
data=sm.stats.anova_lm(RLM, typ=1)
ssq=data[["sum_sq"]]
ssq=np.round(100*ssq/ssq.sum(), 2)
ssq["sum_sq"]=ssq.sum_sq
ssq.columns=["Pourcentage de l'explication de la variance"]
return(ssq)
Dans cette première partie, on se demande si le pays d'origine a une influence sur le revenu de l'individu. En vue d'une ANOVA sur le revenu en ne considérant que le pays d'origine comme variable explicative, nous allons dans un premier temps rappeler un aperçu de la diversité des revenus dans les pays étudiés. Étant donné que les variables "country_code" et "y_child" de "individus" sont celles de "df" répétées 500 fois, nous utiliserons dans cette première partie df au lieu de "individus" en dépit du poids de ce dernier DataFrame, cela ne changera rien à nos résultats.
Nous observerons que les boxplots suivantes indiquent déjà qu'il y a bien un effet "pays" sur les revenus.
fig, ax = plt.subplots(figsize=(14,7))
df.boxplot(column="y_child",by="country_code", ax=ax)
plt.title("Revenus par pays", fontsize=28) ; plt.ylabel("Revenus", fontsize=23) ; plt.xlabel("")
plt.xticks(rotation=90, fontsize=7) ; plt.show()
On modélise alors ici notre ANOVA que nous appellerons "AN". Le DataFrame "individus" étant trop volumineux pour cette exécution, nous la ferons sur df (qui est toujours notre modèle de base avant d'avoir multiplité chaque ligne par 500 au cours de la mission 3), ce qui ne changera rien à nos résultats étant donnée que les salaires sont identiquement répétés 500 fois pour chaque pays dans "individus".
AN = ols("y_child ~ country_code", data=df).fit(alpha=nt) #On rappelle que nt=0,05.
AN.summary().tables[0]
Le coefficient de détermination donné par R-squared dans le modèle indique la performance du modèle. On notera toutefois qu'il diffère, ici, du coefficient de détermination ajusté.
print("Coefficient de détermination:", AN.rsquared)
print("Le pays d'origine explique environ", 100*round(AN.rsquared, 3), "% de la variance de y_child")
Nous allons exécuter un test de Fisher sur nos données. On émet alors les deux hypothèses suivantes:
- H0: Le pays d'origine n'a aucune influence sur le revenu de l'individu.
- H1: Le pays d'origine conditione le revenu de l'individu.
La validation de H0 se confirmera si les différences entre chaque revenu ne sont pas significatives, à contrario de H1 qui sera validé si elles le sont.
sm.stats.anova_lm(AN, typ=1)
La pvalue étant infinitésimalement proche de 0, on peut au niveau de test de 5% rejeter H0 et admettre H1: ainsi on en conclu, selon ce test, que le pays d'origine influe bien sur les revenus.
Distribution des résidus: Nous allons dans un premier temps vérifier l'adéquation de la distribution des résidus de l'ANOVA à la loi normale de paramètres la moyenne et l'écart-type de cette variable, l'hypothèse nulle H0 étant que les résidus ne suivent pas une loi normale, effectuons un test de Kolmogorov-Smirnov sur les résidus.
test_residus_KS(AN)
Le test de Kolmogorov-Smirnov affichant une pvalue infinitésimalement nulle, il ne permet pas le rejet de l'hypothèse nulle. Observons alors la distribution des résidus sur une droite de Henry.
fig, ax = plt.subplots(1, 1, figsize=(6, 5))
residus_henry(AN, "AN", "darkslategrey", ax) ; plt.show()
À la vue de la droite de Henry, les points ne sont pas alignés et l'on remarque un grand nombre de points éloignés de la moyenne, ce qui laisse fortement supposer que les résidus ne suivent pas une loi normale, pas même si l'on supprimait quelques pays à forte influence de nos observations.
Homoscédasticité: Nous allons maintenant observer la variance des résidus. Avec comme hypothèse nulle l'hétéroscédasticité des résidus, un test de Breusch-Pagan sur les résidus nous indiquera ou non son rejet.
test_residus_BP(AN)
Le test de Breusch-Pagan ne nous permet pas le rejet de H0. Observons alors les variances des résidus sur le nuage de variance résiduelle.
fig, ax = plt.subplots(1, 1, figsize=(6, 5))
nuage_residuel(AN, "AN", "darkslategrey", ax) ; plt.show()
Le nuage de la variance résiduelle montre clairement que la variance des résidus n'est pas constante et le test de Breusch-Pagan affichant une pvalue proche de 0, on rejette alors l'homoscédasticité des résidus.
Nous allons ajouter à df la variable "log_y_child" qui sera la variable "y_child" exprimée en logarithme, puis, les résultats de la première ANOVA laissant quelque peu à désirer en dépit de l'analyse des résidus. Nous réaliserons à nouveau chacune des étape précédente pour cette nouvelle ANOVA, à commencer par les boîtes à moustaches.
fig, ax = plt.subplots(figsize=(14,7))
df.boxplot(column="log_y_child",by="country_code", ax=ax)
plt.title("Logarithme des revenus par pays", fontsize=28) ; plt.ylabel("Logarithme des revenus", fontsize=23) ; plt.xlabel("")
plt.xticks(rotation=90, fontsize=7) ; plt.show()
Le passage au logarithme permet effectivement de diminuer le volume de grands nombres représentants une série. On observe toujours une irrégularité de la variance des revenus, mais déjà une distribution plus serrée. Avec comme variable explicative "country_code" et comme varialbe à expliquer "log_y_child", on réalise maintenant l'ANOVA dont nous appellerons "AN2" la modélisation.
AN2 = ols("log_y_child ~ country_code", data=df).fit(alpha=nt) ; AN2.summary().tables[0]
On relève dans le sommaire de la modélisation le R-squared qui évalue la performance du modèle en nous donnant le pourcentage d'explication de la variance des revenus des individus.
print("Coefficient de détermination:", AN2.rsquared)
print("Le pays d'origine explique environ", round(100*AN2.rsquared, 2), "% de la variance de log_y_child")
On peut d'ores et déjà admettre que le passage au logarithme des revenus nous permet de modéliser un modèle bien plus performant que le précédent, avec 72,92% d'explication de variance, contre 49,6% pour le précédent. Soulignons l'effet du pays sur les revenus malgré la réduction de leur variance lors du passage au logarithme par le test de Fisher.
sm.stats.anova_lm(AN2, typ=1)
La pvalue étant nulle, à nouveau, nous concluons également à une influence significative du pays d'origine sur les revenus des individus malgré leur variance réduite, ce qui au passage, appui davantage cette même conclusion pour le modèle précédent. Nous allons maintenant analyser les résidus de cette seconde ANOVA, en testant cette fois-ci simultanément leur distribution par un test de Kolmogorov-Smirnov et leur homoscédasticité par un test de Breusch-Pagan.
test_residus_KS(AN2) # Test de Kolmogorov-Smirnov:
test_residus_BP(AN2) # Test de Breusch-Pagan:
Les pvalues sont effectivement quasi-nulles. On ne peut pour le moment ni rejeter l'idée que la distribution soit gaussienne, ni rejeter son hétéroscédasticité. Néanmoins, observons tout de même encore la droite de Henry et le nuage de la variance résiduelle.
analyse_des_residus(AN2, "AN2", "goldenrod") ; plt.show()
On observe donc sur la droite de Henry un bien meilleur alignement des résidus que pour le modèle précédent. En analysant cette droite, les données semblent suivrent une loi normale excepté pour une poignée d'outliers influants sur le modèle. De plus, on observe sur le nuage résiduel une tendance à la constance des variances. Nous étudierons plus tard dans la mission les individus faussant les résultats des tests statistiques précédents, mais nous pouvons, a priori, nous fier à cette nouvelle ANOVA selon les deux graphiques ci-dessus. Ce modèle est clairement plus fiable que le précédent, en plus d'être plus significatif à la vue des pourcentages d'explication de la variance.
On suppose qu'il existe une relation linéaire entre la variable à expliquer y_child et les variables explicatives mj et gini, donc qu'il existe des paramètres b1, b2 et r tels que: y_child = b1 x mj + b2 x gini + r.
On émet deux à deux les quatre hypothèses de test:
- H0: b1=0 - le paramètre de la variable mj n'est pas significatif.
- H1: b1≠0 - le paramètre de la variable mj est significatif.
- G0: b2=0 - le paramètre de la variable gini n'est pas significatif.
- G1: b2≠0 - le paramètre de la variable gini est significatif.
On effectue maintenant le test de régression linéaire sur y_child selon les variables mj et gini. On appellera "RL1" cette modélisation.
RL1=ols("y_child ~ mj+gini", data=individus).fit(alpha=nt) ; RL1.summary()
Nous obtenons que les 2 pvalues sont de 0 pour chacune des deux variables, ce qui nous conduit à:
- rejeter H0 et admettre H1.
- rejeter G0 et admettre G1.
Les deux variables mj et gini sont alors significatives au niveau de test de 5%. Testons néanmoins leur éventuelle colinéarité.
VIF(RL1, ["mj", "gini"])
Les deux facteurs d'influence sur la variance sont inférieurs à 10. Il n'y a donc pas de problème de colinéarité dans notre modèle.
On suppose qu'il existe une relation linéaire entre la variable à expliquer log_y_child et les variables explicatives log_mj et gini, donc qu'il existe des paramètres b1, b2 et r tels que: log_y_child = b1 x log_mj + b2 x gini + r.
On émet deux à deux les quatre hypothèses de test:
- H0: b1=0 - le paramètre de la variable log_mj n'est pas significatif.
- H1: b1≠0 - le paramètre de la variable log_mj est significatif.
- G0: b2=0 - le paramètre de la variable gini n'est pas significatif.
- G1: b2≠0 - le paramètre de la variable gini est significatif.
On effectue maintenant le test de régression linéaire sur log_y_child selon les variables log_mj et gini. On appelle ce second modèle "RL2", et "srl_2" son sommaire.
RL2=ols("log_y_child ~ log_mj+gini", data=individus).fit(alpha=nt) ; RL2.summary()
Nous obtenons que les 2 pvalues sont de 0 pour chacune des deux variables, donc b1 et b2 sont significativement différents de 0, ce qui nous conduit à:
- rejeter H0 et admettre H1.
- rejeter G0 et admettre G1.
Les deux variables log_mj et gini sont alors significatives au niveau de test de 5%. Testons également leur colinéarité.
VIF(RL2, ["log_mj", "gini"])
Les deux facteurs d'influence sur la variance sont inférieurs à 10. Il n'y a donc pas de problème de colinéarité dans notre modèle non plus.
Nous allons analyser le pourcentage de la variance expliquée des revenus par chacun des deux derniers modèles de régressions linéraires RL1 où les données sont inchangées et RL2 où les données sont exprimées en logarithme. Les droites de Henry sur les résidus de chacune des deux régressions linéaires précédentes nous indiqueront lesquels sont le plus en adéquation avec la loi gaussienne.
# Analyse de la première régression linéaire sans passer au logarithme:
print("COEFFICIENT DE DÉTERMINATION DE RL1 - SANS PASSER AU LOGARITHME:", round(RL1.rsquared, 4))
analyse_des_residus(RL1, "RL1", "purple") ; plt.show()
print("_"*60) ; print(" ")
# Analyse de la seconde régression linéaire en passant au logarithme:
print("COEFFICIENT DE DÉTERMINATION DE RL2 - EN PASSANT AU LOGARITHME:", round(RL2.rsquared, 4))
analyse_des_residus(RL2, "RL2", "mediumblue") ; plt.show()
Le pourcentage de la variance expliquée par le premier modèle est d'environ 44,77% tandis que celui du second modèle est proche de 65,18%, donc un coefficient de détermination plus élevé pour RL2 que pour RL1.
Les droites de Henry indiquent clairement que les résidus du second modèle (en passant au logarithme) sont bien plus en adéquation avec une loi gaussienne que ne le sont les résidus du premier modèle:
- Pour le premier modèle, sans passer au logarithme, nous observons un éparpillement total des résidus, comme pour la dernière modélisation, aucunement une adéquation à la loi gaussienne, et nul échantillonnage raisonnable de nos individus ne semblerait y remédier.
- Pour la deuxième régression linéaire, en revanche, où l'on passe au logarithme, la distribution empirique des résidus semble suivre une loi gaussienne selon la droite de Henry, excepté pour les valeurs faibles des résidus.
Les nuages de variance résiduelle indiquent clairement que les résidus du second modèle (en passant au logarithme) sont bien plus proches:
- Pour le premier modèle, sans passer au logarithme, les variances diffèrent significativement.
- Pour la deuxième régression linéaire, en revanche, où l'on passe au logarithme, les variances semblent toutes plus ou moins égales entre elles.
C'est donc indégnablement ce second modèle que nous considérerons pour la suite, le premier n'étant pas viable. Le passage au logarithme nous a donc bien permis un modèle de régression linéaire plus fiable que le premier réalisé.
On s'intéresse maintenant au pourcentage de la variance de "log_y_child" expliquée par chacune des deux variables "log_mj" et "gini". Pour les obtenir, nous modéliserons la table de la régression linéaire de RL2 que nous appellerons "table_RL2". Nous appliquons ici sur RL2 la fonction "decomposition_des_variables" définie dans la partie I de la mission 4 pour observer le pourcentage de l'explication de la variance de "log_y_child" de chaque variable du modèle.
table_RL2=sm.stats.anova_lm(RL2, typ=2)
decomposition_des_variables(RL2)
Ainsi, log_mj explique 64,33% de la variance de log_y_child tandis que la variable gini n'apporte qu'une contribution de 0,85%. On observe sur les résidus que cette variance est expliquée à 34,82% par d'autres facteurs.
Notre dernier modèle de régression linéaire explique environ 65,18% de la variance. Voyons si la variable "c_i_parents" contribue à nous l'expliquer davantage. On modélise donc une nouvelle régression linéaire "RL3", qui reprend "RL2" en y ajoutant les classes de revenu des parents des individus.
RL3=ols("log_y_child ~ log_mj+gini+c_i_parents", data=individus).fit(alpha=nt) ; RL3.summary()
Le sommaire de la régression linéaire affiche un coefficient de détermination de 0,703, nous parvenons donc maintenant à expliquer environ 70,3% de la variance grâce à l'ajout de cette variable. Notons que le coefficient de détermination est, ici, ajusté. Analysons maintenant la décomposition de l'explication de la variance du revenu pour RL3.
decomposition_des_variables(RL3)
En effet, la classe de revenus des parents explique 5,12% de la variance des revenus des enfants selon RL3. Analysons maintenant les facteurs d'influence de la variance de chacune des 3 variables.
VIF(RL3, ["log_mj", "gini", "c_i_parents"])
Les VIF étant tous inférieurs à 10, nous ne détectons donc aucune problème de colinéarité. Observons maintenant la droite de Henry liée aux résidus de cette dernière régression linéaire.
fig, ax = plt.subplots(1, 1, figsize=(6, 5))
residus_henry(RL3, "RL3", "darkorange", ax) ; plt.show()
Les résidus semblent bien suivrent une loi gaussienne selon la droite de Henry, excepté encore une fois pour quelques individus à faible résidus. On remarquera une adéquation à la loi normale encore meilleure que pour RL2. Analysont maintenant l'homoscédasticité de RL3 avec le nuage de variance résiduelle.
fig, ax = plt.subplots(1, 1, figsize=(6, 5))
nuage_residuel(RL3, "RL3", "darkorange", ax) ; plt.show()
Les variances des résidus sont significativement égales selon ce nuage, on en déduit l'homoscédasticité de RL3 et attestons ainsi, de par ce résultat et les précédents, de la fiabilité du modèle de régression linéaire RL3.
Nous allons maintenant vérifier si d'autres variables présentes dans nos données de base influent sur la variance des revenus des individus.
RL4=ols("log_y_child ~ log_mj+gini+c_i_parents+log_y_parents+elasticity+population", data=individus).fit(alpha=nt)
print(RL4.summary()) ; print("") ; print("")
print("CE NOUVEAU MODELE EXPLIQUE ENVIRON", round(100*RL4.rsquared, 2), "% DE LA VARIANCE DE log_y_child")
La pvalue de chaque variable étant égale à 0, elles semblent alors toutes significatives pour le modèle, et ce nouveau modèle explique environ 78,42% de la variance des revenus, soit 8,12% de plus que RL3. Avant d'analyser la décomposition des variables, vérifions d'abord leurs éventuelles colinéarités.
VIF(RL4, ["log_mj", "gini", "c_i_parents", "log_y_parents", "elasticity", "population"])
Nous obtenons des facteurs d'influence plus grands, cette fois-ci, dont certains sont supérieurs à 5, mais aucun n'est inférieur à 10. On constate alors une certaine colinéarité pour quelques variables mais on peut admettre qu'elles n'influent pas sur le modèle. On peut donc observer la décomposition de la variance des revenus expliquée par ces variables.
decomposition_des_variables(RL4)
Il est intéressant de noter que le revenu même des parents vient s'ajouter à l'explication de la variance des revenus, avec un score de 8,12% dans notre dernier modèle, et que nous pouvons mettre de côté les variables "population" et "elasticity" selon ce modèle. On peut donc reconsidérer RL4 sans ces deux dernières variables.
RL4=ols("log_y_child ~ log_mj+gini+c_i_parents+log_y_parents", data=individus).fit(alpha=nt)
print("Et ce modèle explique toujours", round(100*RL4.rsquared, 2), "% de la variance des revenus.")
decomposition_des_variables(RL4)
Analysons les résidus pour nous assurer de la validité de cette dernière régression linéaire multiple.
analyse_des_residus(RL4, "RL4", "chartreuse") ; plt.show()
Les résidus semblent moins confondues avec la droite de Henry mais présentent tout de même une distribution sensiblement proche de la loi gaussienne, et les variances restent sensiblement proches sur le nuage pour ce qui est de l'homoscédasticité. On considérera alors ce dernier modèle comme étant fiable et nous nous reposerons dessus pour la conclusion.
En se basant sur notre dernière régression linéaire, on établi ci-dessous la liste descendante des variables explicant la variance des revenus des individus selon leur pourcentage de contribution.
decomposition_des_variables(RL4).head(4).sort_values(by="Pourcentage de l'explication de la variance", ascending=False)
Et nous avons vu que ce modèle explique environ 78,42% de cette variance, ce qui signifie qu'environ 21,58% de cette variance s'explique par des facteurs extérieurs (efforts, chances, etc...).
En passant au logarithme, nous avons pu expliquer jusqu'à 78,42% de la variance des revenus avec les données présentes dans notre base. En observant les résultats de la mission 4, on peut admettre que le revenu d'un individu est conditionné en priorité par le revenu moyen de son pays d'origine, représentant plus de 70% de cette variance. Cette dernière s'explique au second plan par les conditions financières de leurs parents, d'abord selon leur revenu même, puis de leur classe de revenu au sein de leur pays. Enfin, l'observation de l'indice de Gini est un facteur faible mais qui, après analyse des données précédentes, peut être déterminant. Enfin, tout facteur extérieur à ceux étudiés dans ce projet (effort, chance, potentiel...) contribue à 21,58% de la variance des revenus.
Dans le scénario du projet, nous travaillons dans une banque où l'on a accès à chacune des quatre variables explicatives définies pour chacun de nos clients: le revenu moyen du pays où ils habitent, leur revenu, leur classe de revenu au sein de leur pays, et l'indice de Gini de leur pays. Voici donc un algorithme qui repose sur notre base de données qui, pour un client donné, prendrait en argument le revenu moyen de son pays, son revenu, sa classe de revenu et l'indice de gini de son pays et qui nous renvoie une estimation du revenu de ses enfants.
# On se base uniquement sur les 4 variables en vigueurs présentes dans notre base:
donnees=individus[["y_child", "mj", "y_parents", "c_i_parents", "gini"]]
# Algorithme préliminaire pour l'algorithme de prédiction:
def proche(x, L):
K=(abs(x-L))
a=min(K)
dt=pd.concat([L, K], axis=1)
dt.columns=["vari", "abso"]
dt.drop_duplicates(inplace=True)
i=list(dt.loc[dt["abso"]==a]["vari"])[0]
return(i)
# Algorithme de prédiction:
def prediction(revenu_moyen_pays, revenu_parent, classe_parent, gini):
rmp=proche(revenu_moyen_pays, donnees["mj"])
data=donnees.loc[donnees["mj"]==rmp]
rp=proche(revenu_parent, data["y_parents"])
data=data.loc[data["y_parents"]==rp]
cp=proche(classe_parent, data["c_i_parents"])
data=data.loc[data["c_i_parents"]==cp]
g=proche(gini, data["gini"])
data=data.loc[data["gini"]==g]
return(round(data["y_child"].mean(), 2))
Imaginons qu'un client gagne environ 3000 dollars ppa par mois, appartienne au 64ème centile, vit dans un pays où le revenu moyen est de 2500 dollars ppa, et où l'indice de gini est de 0,35. On propose une estimation du revenu de son enfant.
prediction(2500, 3000, 60, 0.35)