P7_04_Mission 3 & 4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4c9e163e",
   "metadata": {},
   "source": [
    "# Import des librairies, dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b33afda2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as st\n",
    "import plotly.express as px\n",
    "import plotly\n",
    "import plotly.graph_objs as go\n",
    "from statsmodels.stats.outliers_influence import variance_inflation_factor\n",
    "\n",
    "import pylab \n",
    "import math\n",
    "\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.api as sm\n",
    "import seaborn as sns\n",
    "import scipy.stats as stats\n",
    "from statsmodels.stats.stattools import durbin_watson\n",
    "\n",
    "# pour afficher dynamiquement dans le notebook\n",
    "from IPython.display import clear_output\n",
    "from IPython.display import display, Math, Markdown\n",
    "import time\n",
    "\n",
    "from datetime import datetime\n",
    "\n",
    "pd.set_option('display.max_columns', None) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77e42346",
   "metadata": {},
   "outputs": [],
   "source": [
    "echantillon = pd.read_csv(\"données/echantillon.csv\")\n",
    "echantillon = echantillon.drop(columns=['Unnamed: 0'])\n",
    "data_ok = pd.read_csv(\"données/data_ok.csv\")\n",
    "data_ok = data_ok.drop(columns=['Unnamed: 0'])\n",
    "gini_4 = pd.read_csv(\"données/gini_4.csv\")\n",
    "dist_revenus_ok = pd.read_csv(\"données/dist_revenus_ok.csv\")\n",
    "dist_revenus_ok = dist_revenus_ok.drop(columns=['Unnamed: 0'])\n",
    "dist_revenus_ok = dist_revenus_ok.rename(columns = {\"country\" : \"iso3\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b05ccd93",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_rev_moy = dist_revenus_ok.groupby(by = \"iso3\").mean()\n",
    "dist_rev_moy = dist_rev_moy.rename(columns = {\"income\": \"income_mean\"}).reset_index()\n",
    "dist_rev_moy = dist_rev_moy[[\"iso3\", \"income_mean\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da36bba1",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_revenus_ok = pd.merge(left = dist_rev_moy, right = dist_revenus_ok, on = \"iso3\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59db0580",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_revenus_ok = dist_revenus_ok[[\"iso3\", \"income_mean\", \"quantile\", \"income\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ece1ffe",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_revenus_ok"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4b75862",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Merge du dataset\n",
    "#dist_revenus_ok = dist_revenus_ok.rename(columns = {\"country\" : \"iso3\"})\n",
    "données_completes = pd.merge(left = data_ok, right = dist_revenus_ok, on = \"iso3\", how = \"outer\")\n",
    "données_completes = données_completes[données_completes['country'].notna()]\n",
    "WID = données_completes[[\"country\", \"quantile\", \"income\", \"income_mean\", \"moy_gini\", \"elasticite_ok\"]]\n",
    "WID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "018ee787",
   "metadata": {},
   "outputs": [],
   "source": [
    "WID.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "442941f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(WID['country'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6860e28",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(WID.isnull().any())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3232fb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(WID.duplicated().any())\n",
    "print(WID.duplicated().sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa0da594",
   "metadata": {},
   "source": [
    "# Mission 3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "723f726a",
   "metadata": {},
   "source": [
    "## Définition des fonctions utiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3de875e9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "\n",
    "def generate_incomes(n, pj):\n",
    "    # On génère les revenus des parents (exprimés en logs) selon une loi normale.\n",
    "    # La moyenne et variance n'ont aucune incidence sur le résultat final (ie. sur le caclul de la classe de revenu)\n",
    "    ln_y_parent = st.norm(0,1).rvs(size=n)\n",
    "    # Génération d'une réalisation du terme d'erreur epsilon\n",
    "    residues = st.norm(0,1).rvs(size=n)\n",
    "    return np.exp(pj*ln_y_parent + residues), np.exp(ln_y_parent)\n",
    "\n",
    "def quantiles(l, nb_quantiles):\n",
    "    size = len(l)\n",
    "    l_sorted = l.copy()\n",
    "    l_sorted = l_sorted.sort_values()\n",
    "    quantiles = np.round(np.arange(1, nb_quantiles+1, nb_quantiles/size) -0.5 +1./size)\n",
    "    q_dict = {a:int(b) for a,b in zip(l_sorted,quantiles)}\n",
    "    return pd.Series([q_dict[e] for e in l])\n",
    "\n",
    "def compute_quantiles(y_child, y_parents, nb_quantiles):\n",
    "    y_child = pd.Series(y_child)\n",
    "    y_parents = pd.Series(y_parents)\n",
    "    c_i_child = quantiles(y_child, nb_quantiles)\n",
    "    c_i_parent = quantiles(y_parents, nb_quantiles)\n",
    "    sample = pd.concat([y_child, y_parents, c_i_child, c_i_parent], axis=1)\n",
    "    sample.columns = [\"y_child\", \"y_parents\", \"c_i_child\",\"c_i_parent\"]\n",
    "    return sample\n",
    "\n",
    "def distribution(counts, nb_quantiles):\n",
    "    distrib = []\n",
    "    total = counts[\"counts\"].sum()\n",
    "    if total == 0 :\n",
    "        return [0] * nb_quantiles\n",
    "    for q_p in range(1, nb_quantiles+1):\n",
    "        subset = counts[counts.c_i_parent == q_p]\n",
    "        if len(subset):\n",
    "            nb = subset[\"counts\"].values[0]\n",
    "            distrib += [nb / total]\n",
    "        else:\n",
    "            distrib += [0]\n",
    "    return distrib   \n",
    "\n",
    "def conditional_distributions(sample, nb_quantiles):\n",
    "    counts = sample.groupby([\"c_i_child\",\"c_i_parent\"]).apply(len)\n",
    "    counts = counts.reset_index()\n",
    "    counts.columns = [\"c_i_child\",\"c_i_parent\",\"counts\"]\n",
    "    mat = []\n",
    "    for child_quantile in np.arange(nb_quantiles)+1:\n",
    "        subset = counts[counts.c_i_child == child_quantile]\n",
    "        mat += [distribution(subset, nb_quantiles)]\n",
    "    return np.array(mat) \n",
    "\n",
    "def plot_conditional_distributions(p, cd, nb_quantiles):\n",
    "    #plt.figure(figsize=(10,10))\n",
    "    # La ligne suivante sert à afficher un graphique en \"stack bars\", sur ce modèle : https://matplotlib.org/gallery/lines_bars_and_markers/bar_stacked.html\n",
    "    cumul = np.array([0] * nb_quantiles)\n",
    "    for i, child_quantile in enumerate(cd):\n",
    "        plt.bar(np.arange(nb_quantiles)+1, child_quantile, bottom=cumul, width=0.95, label = str(i+1) +\"e\")\n",
    "        cumul = cumul + np.array(child_quantile)\n",
    "    plt.axis([.5, nb_quantiles*1.3 ,0 ,1])\n",
    "    plt.title(\"p=\" + str(p))\n",
    "    plt.legend(ncol=1)\n",
    "    plt.xlabel(\"quantile parents\")\n",
    "    plt.ylabel(\"probabilité du quantile enfant\")\n",
    "    plt.show()\n",
    "    \n",
    "def proba_cond(c_i_parent, c_i_child, mat):\n",
    "    return mat[c_i_child, c_i_parent]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3445d27e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def smooth(x,y, box_percent=0.05,res=50,median=True):\n",
    "    surface = max(x)-min(x)\n",
    "    my_pas = np.arange(min(x),max(x),surface/res)\n",
    "    box = surface*box_percent\n",
    "    demi_box = box/2\n",
    "    y_sortie = np.array([])\n",
    "    x_sortie = np.array([])\n",
    "    for myx in my_pas :\n",
    "        temp = [y[i] for i in range(len(x)) if ((x[i]>=(myx-demi_box))and(x[i]<=(myx+demi_box)))]\n",
    "        if median==True :\n",
    "            temp_y = np.median(temp)\n",
    "        else :\n",
    "            temp_y = np.mean(temp)\n",
    "        #print(temp_y)\n",
    "        y_sortie = np.append(y_sortie,temp_y)\n",
    "        #print(y_sortie)\n",
    "        x_sortie = np.append(x_sortie,myx)\n",
    "    return x_sortie, y_sortie"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50296514",
   "metadata": {},
   "source": [
    "### Question 1 à 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbc1fce6",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "#Question 1 à 2 :\n",
    "pj = 0.9\n",
    "nb_quantiles = 10\n",
    "n = 1000*nb_quantiles\n",
    "\n",
    "#Question 3 : \n",
    "y_child, y_parents = generate_incomes(n, pj)\n",
    "\n",
    "#Question 4 :\n",
    "sample = compute_quantiles(y_child, y_parents, nb_quantiles)\n",
    "\n",
    "#Question 5 :\n",
    "cd = conditional_distributions(sample, nb_quantiles)\n",
    "\n",
    "#plot_conditional_distributions(pj, cd, nb_quantiles) #prend beaucoup de temps\n",
    "\n",
    "c_i_child = 5\n",
    "c_i_parent = 8\n",
    "p = proba_cond(c_i_parent, c_i_child, cd)\n",
    "print(\"\\nP(c_i_parent = {} | c_i_child = {}, pj = {}) = {}\".format(c_i_parent, \n",
    "                                                                   c_i_child,pj, p))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96cb5194",
   "metadata": {},
   "source": [
    "### Question 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac6f6cc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "pj_2 = 0.9 \n",
    "nb_quantiles_2 = 10\n",
    "n_2 = 1000*nb_quantiles_2\n",
    "y_child_2, y_parents_2 = generate_incomes(n_2, pj_2)\n",
    "sample_2 = compute_quantiles(y_child_2, y_parents_2, nb_quantiles_2)\n",
    "cd_2 = conditional_distributions(sample_2, nb_quantiles_2)\n",
    "plot_conditional_distributions(pj_2, cd_2, nb_quantiles_2) # Cette instruction prendra du temps si nb_quantiles > 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d69bebcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "pj_3 = 0.1 \n",
    "nb_quantiles_3 = 10\n",
    "n_3 = 1000*nb_quantiles_3\n",
    "y_child_3, y_parents_3 = generate_incomes(n_3, pj_3)\n",
    "sample_3 = compute_quantiles(y_child_3, y_parents_3, nb_quantiles_3)\n",
    "cd_3 = conditional_distributions(sample_3, nb_quantiles_3)\n",
    "plot_conditional_distributions(pj_3, cd_3, nb_quantiles_3) # Cette instruction prendra du temps si nb_quantiles > 10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dd9b580",
   "metadata": {},
   "source": [
    "### Question 7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84cc8a44",
   "metadata": {},
   "outputs": [],
   "source": [
    "#On supprime les individus créés, on ne garde que les distributions conditionnelles \"cd\"\n",
    "del pj, nb_quantiles, n, y_child, y_parents, sample, c_i_child, c_i_parent, p, cd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcb3302f",
   "metadata": {},
   "source": [
    "### Question 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e08f8865",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Multiplication par 500\n",
    "data_cloned = pd.concat([WID]*500, ignore_index=True)\n",
    "\n",
    "print('WID shape :', WID.shape)\n",
    "print('data_cloned shape :', data_cloned.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f3cd81b",
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in data_cloned.columns:\n",
    "    print(col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89f110b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cloned.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d92205a",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cloned = data_cloned[[\"country\", \"quantile\", \"income\", \"moy_gini\", \"elasticite_ok\"]]\n",
    "data_cloned.rename(columns={'quantile': 'c_i_child', 'income': 'y_child', 'moy_gini': 'G_j', 'elasticite_ok': 'p_j'}, inplace=True)\n",
    "\n",
    "data_cloned.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90bfa7f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cloned.sample(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc23e97d",
   "metadata": {},
   "source": [
    "### Question 9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5dcb371",
   "metadata": {},
   "outputs": [],
   "source": [
    "country_list = data_cloned['country'].unique()\n",
    "len(country_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4af45876",
   "metadata": {},
   "outputs": [],
   "source": [
    "# creation d'un variable quantile parents, à remplir\n",
    "list_proba = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3799130",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "#Pour chaque pays dans la liste\n",
    "for country in country_list :\n",
    "    \n",
    "    #coef d'elasticité pour chaque pays de data_cloned dans la liste country_list, premiere ligne\n",
    "    pj = data_cloned.loc[data_cloned['country'] == country,'p_j'].iloc[0] \n",
    "    \n",
    "    #nombre de quantiles (nombre de classes de revenu)\n",
    "    nb_quantiles = 100 \n",
    "    \n",
    "    #taille de l'échantillon\n",
    "    n = 50000 \n",
    "    \n",
    "    #Génération de revenus selons une loi normale\n",
    "    y_child, y_parents = generate_incomes(n, pj) \n",
    "    \n",
    "    #Retourne un df avec y_child, y_parents, c_i_child, c_i_parents\n",
    "    sample = compute_quantiles(y_child, y_parents, nb_quantiles) \n",
    "    \n",
    "    #Calcul & attribution des probabilités conditonnelles \n",
    "    cd = conditional_distributions(sample, nb_quantiles)\n",
    "    \n",
    "    #On compte chaque combinaison c_i_child, c_i_parent\n",
    "    for c_i_child in range(100):\n",
    "        for c_i_parent in range(100):\n",
    "            p = proba_cond(c_i_parent, c_i_child, cd)\n",
    "            #print(\"\\nP(c_i_parent = {} | c_i_child = {}, pj = {}) = {}\".format(c_i_parent, c_i_child,pj,p))\n",
    "            \n",
    "            #Association des probas conditionelles aux individus\n",
    "            list_proba.extend([c_i_parent+1]*(int(p*500)))\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06c3d2d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb286238",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83bd82cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mesure de la mobilité\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5569c9af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Je check que ma liste à le meme nombre de ligne que mon data_cloned\n",
    "len(list_proba)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97ced9d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Je créé une colonne dans mon data_cloned pour la classe parents\n",
    "data_cloned['proba'] = list_proba\n",
    "data_cloned = data_cloned.rename({'proba':'c_i_parent'},axis=1)\n",
    "data_cloned.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5be569f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Insertion du revenu moyen par pays\n",
    "mean_income = data_cloned.groupby(by='country').mean()\n",
    "mean_income.reset_index(inplace=True)\n",
    "mean_income = mean_income[['country', 'y_child']]\n",
    "mean_income.rename(columns={'y_child': 'm_j'}, inplace=True)\n",
    "\n",
    "mean_income.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a553c31a",
   "metadata": {},
   "source": [
    " #Insertion du revenu moyen par pays\n",
    "mean_income = WID[[\"country\", \"income_mean\"]].drop_duplicates().reset_index()\n",
    "mean_income.rename(columns={'income_mean': 'm_j'}, inplace=True)\n",
    "mean_income = mean_income[[\"country\", \"m_j\"]]\n",
    "mean_income.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1dd9658",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Jointure de mean_income avec data_cloned & suppression de c_i_child\n",
    "data_cloned_b = pd.merge(data_cloned, mean_income, on='country')\n",
    "data_cloned_b = data_cloned_b[[\"country\", \"y_child\", \"G_j\", \"p_j\", \"c_i_parent\", \"m_j\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed52d774",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Insertion de colonnes log de y_child et m_j\n",
    "data_cloned_b[\"income_log\"] = np.log(data_cloned_b[\"y_child\"])\n",
    "data_cloned_b[\"m_j_log\"] = np.log(data_cloned_b[\"m_j\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b8c3e5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cloned_b.loc[data_cloned_b[\"country\"] == \"France\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14fe9bcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cloned_b.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c1d9674",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sauvegarde du dataframe\n",
    "#data_cloned_b.to_csv(\"données/data_cloned_b.csv\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d482b270",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cloned_b.sample(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68a3e56d",
   "metadata": {},
   "source": [
    "# Mission 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "796515fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "echantillon_m4 = data_cloned_b.loc[(data_cloned_b[\"country\"] == \"France\")\n",
    "                     | (data_cloned_b[\"country\"] == \"Chile\")\n",
    "                     | (data_cloned_b[\"country\"] == \"Iceland\") \n",
    "                     | (data_cloned_b[\"country\"] == \"Paraguay\")\n",
    "                     | (data_cloned_b[\"country\"] == \"Vietnam\") \n",
    "                     | (data_cloned_b[\"country\"] == \"Congo, Dem. Rep.\")\n",
    "                    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2ebdcaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "echantillon_m4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8beb84a",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = echantillon_m4\n",
    "fig = px.box(df, x=\"y_child\", y=\"country\", color=\"country\")\n",
    "\n",
    "fig.show()\n",
    "\n",
    "print(\"Les distributions des revenus des individus varient fortement entre les pays dans l'échantillon choisi\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c53fc50",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = echantillon_m4\n",
    "fig = px.box(df, x=\"income_log\", y=\"country\", color=\"country\")\n",
    "\n",
    "fig.show()\n",
    "\n",
    "print(\"Les distributions des revenus en log des individus varient fortement entre les pays dans l'échantillon choisi\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12e9c6e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(30,5))\n",
    "rs = plt.plot(data_cloned_b.groupby('country').y_child.mean(), 'o')\n",
    "rs = plt.xticks(rotation=90)\n",
    "plt.xlabel('pays')\n",
    "plt.ylabel('revenu moyen des enfants')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f89c33e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(30,5))\n",
    "rs = plt.plot(data_cloned_b.groupby('country').income_log.mean(), 'o')\n",
    "rs = plt.xticks(rotation=90)\n",
    "plt.xlabel('pays')\n",
    "plt.ylabel('revenu moyen des enfants')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fa83f80",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Agrégation pour gagner en temps de calcul\n",
    "data_cloned_c = data_cloned_b.groupby(by=['country',\n",
    "                                  'y_child',\n",
    "                                  'm_j',\n",
    "                                  'G_j',\n",
    "                                  'p_j',\n",
    "                                  'm_j_log',\n",
    "                                  'income_log']).mean()\n",
    "data_cloned_c.reset_index(inplace=True)\n",
    "#data_projet7_m4_2.drop(columns=['c_i_parent'], inplace=True) # variable non necessaire pour la suite\n",
    "\n",
    "data_cloned_c.c_i_parent = data_cloned_c.c_i_parent.round()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f95a3604",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('data_cloned_c shape :', data_cloned_c.shape)\n",
    "print('data_cloned_b shape :', data_cloned_b.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10fe324f",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cloned_c.sample(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2dd4ca9",
   "metadata": {},
   "source": [
    "## Anova"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5f2fd1d",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "anova_pays_rar = smf.ols('y_child~country', data=data_cloned_c).fit()\n",
    "print(anova_pays_rar.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b84504f",
   "metadata": {},
   "source": [
    "La pvalue du F test (test de Fisher) étant proche de 0, on peut au niveau de test de 5% rejeter H0 et admettre H1: ainsi on peut conclure que le pays d'origine influe bien sur les revenus. </br>\n",
    "Par ailleurs, d'après ce modèle, le pays d'origine explique **50%** de la variance du revenu."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f96a1656",
   "metadata": {},
   "source": [
    "### Tests"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c76c71c",
   "metadata": {},
   "source": [
    "**Test 1 : Distribution des résidus** </br>\n",
    "On vérifie l'adéquation de la distribution des résidus de l'ANOVA à une loi normale. l'hypothèse nulle H0 étant que les résidus suivent une loi normale, on procède à un test de Kolmogorov-Smirnov sur les résidus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b417051e",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = anova_pays_rar.resid\n",
    "stats.kstest(x, 'norm')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbd7a111",
   "metadata": {},
   "source": [
    "HO : les données suivent une loi normale </br>\n",
    "H1 : les données ne suivent pas une loi normale </br>\n",
    "pvalue < 0,05 : on rejette H0, les données ne suivent pas une loi normale"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2844bc56",
   "metadata": {},
   "source": [
    "Le KS test affichant une pvalue nulle, les données ne suivent pas une loi normale. </br>\n",
    "On observe la distribution des résidus sur une droite de Henry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "653ab4cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.probplot(x, dist=\"norm\", plot=pylab)\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d204f66",
   "metadata": {},
   "source": [
    "On a graphiquement une part élevée de résidus s'éloignant de la droite théorique, donc on peut dire que les résidus ne suivent vraisemblablement pas une loi normale."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75188a0d",
   "metadata": {},
   "source": [
    "**Test 2 : Homoscédasticité** </br>\n",
    "On observe la variance des résidus.\n",
    "HO l'hypothèse nulle d'hétéroscédasticité des résidus, on réalise un test de Breusch-Pagan sur les résidus.\n",
    "\n",
    "H0 : homoscédasticité\n",
    "H1 : hétéroscédasticité"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9c428a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sm.stats.diagnostic.het_breuschpagan(anova_pays_rar.resid, anova_pays_rar.model.exog))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "908b3b45",
   "metadata": {},
   "source": [
    "Avec une p value proche de 0 (seconde valeur), on admet l'hétéroscédasticité. </br>\n",
    "On peut observer les variances des résidus sur le nuage de variance résiduelle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "909115d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=plt.plot(anova_pays_rar.fittedvalues, anova_pays_rar.resid , \".\",  alpha=0.3)\n",
    "plt.title(\"Nuage de la variance résiduelle\", fontsize=18)\n",
    "#plt.xlabel(\"GWh\", fontsize=16), plt.ylabel(\"Résidus\", fontsize=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abdb5506",
   "metadata": {},
   "source": [
    "On peut confirmer la non linéarité avec un **test de Rainbow**, qui vérifie H0 : la représentation statistique est bien linéaire. </br> \n",
    "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é. </br>\n",
    "Ici, la p value étant inférieure au seuil, on rejette H0. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a528492c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.stats.diagnostic import linear_rainbow\n",
    "Ftest, pval = linear_rainbow(anova_pays_rar)\n",
    "print(pval)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a34fd45",
   "metadata": {},
   "source": [
    "En conclusion, si ce modèle explique bien 50% de la variance du revenu des individus avec le pays d'origine, il ne semble pas être robuste aux tests usuels et ne peut donc être utilisé pour une prédiction pertinente."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be705ffa",
   "metadata": {},
   "source": [
    "## Régression linéaire avec le revenu moyen et l'indice de gini"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96c0b605",
   "metadata": {},
   "source": [
    "### Modèle 1 : revenu & gini"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12aa922b",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Variables explicatives : revenu moyen du pays de l’individu et l’indice de Gini du pays de l’individu\n",
    "reg_v1_1 = smf.ols('y_child ~ G_j+m_j', data=data_cloned_c).fit()\n",
    "print(reg_v1_1.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a660681",
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_v1_1_sct = reg_v1_1.centered_tss\n",
    "reg_v1_1_sce = reg_v1_1.ess\n",
    "reg_v1_1_scr = reg_v1_1.ssr\n",
    "reg_v1_1_r2 = reg_v1_1.rsquared\n",
    "\n",
    "display(Markdown(f\"\"\"\n",
    "Décomposition de la variance : <br>\n",
    "\n",
    "**SCT = SCE + SCR**\n",
    "- SCT = {reg_v1_1_sct:.0f} <br>\n",
    "- SCE = {reg_v1_1_sce:.0f} <br>\n",
    "- SCR = {reg_v1_1_scr:.0f} <br>\n",
    "\n",
    "$R^2 = {reg_v1_1_r2:.2f}$\n",
    "\n",
    "Le modèle explique **{reg_v1_1_r2*100:.0f}% de la variance totale**.\n",
    "\"\"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62648dc8",
   "metadata": {},
   "source": [
    "Selon ce modèle, le pays de naissance (ie. le revenu moyen et l’indice de Gini) explique donc 50% de la variance totale, tandis que les autres facteurs non considérés dans le modèle (efforts, chance, etc.) représentent l'autre moitié.\n",
    "Par ailleurs, l'indice de Gini n'est pas statistiquement significatif au seuil de 0,05."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d39c9bce",
   "metadata": {},
   "source": [
    "### Modèle 2 : revenus en log & indice de gini"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "376b2500",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Variables explicatives : revenu moyen du pays de l’individu en log et l’indice de Gini du pays de l’individu\n",
    "reg_v1_2 = smf.ols('income_log ~ G_j+m_j_log', data=data_cloned_c).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a88ef0ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_v1_2_sct = reg_v1_2.centered_tss\n",
    "reg_v1_2_sce = reg_v1_2.ess\n",
    "reg_v1_2_scr = reg_v1_2.ssr\n",
    "reg_v1_2_r2 = reg_v1_2.rsquared\n",
    "\n",
    "display(Markdown(f\"\"\"\n",
    "Décomposition de la variance : <br>\n",
    "\n",
    "**SCT = SCE + SCR**\n",
    "- SCT = {reg_v1_2_sct:.0f} <br>\n",
    "- SCE = {reg_v1_2_sce:.0f} <br>\n",
    "- SCR = {reg_v1_2_scr:.0f} <br>\n",
    "\n",
    "$R^2 = {reg_v1_2_r2:.2f}$\n",
    "\n",
    "Le modèle explique **{reg_v1_2_r2*100:.0f}% de la variance totale**.\n",
    "\"\"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c44873a",
   "metadata": {},
   "source": [
    "Selon ce modèle, le pays de naissance (ie. le revenu moyen et l’indice de Gini) explique donc **73%** de la variance totale, tandis que les autres facteurs non considérés dans le modèle (efforts, chance, etc.) représentent les **27%** restant.\n",
    "Le passage en log des variables de revenus améliore le pouvoir explicatif du modèle, et l'indice de Gini devient significatif. Cependant, l'interprétation est rendue plus compliquée du fait des logs (pas de \"lecture naturelle\" posible)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78be0356",
   "metadata": {},
   "source": [
    "## Régression linéaire avec le revenu moyen, l'indice de gini et la classe de revenu des parents"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a6f8a4c",
   "metadata": {},
   "source": [
    "### Modèle 3 : revenu, gini & classe de revenus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04829604",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Variables explicatives : revenu moyen du pays de l’individu indice de Gini du pays de l’individu, classe de revenu des parents\n",
    "reg_v2_1 = smf.ols('y_child ~ G_j+m_j+c_i_parent', data=data_cloned_c).fit()\n",
    "print(reg_v2_1.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1376e75d",
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_v2_1_sct = reg_v2_1.centered_tss\n",
    "reg_v2_1_sce = reg_v2_1.ess\n",
    "reg_v2_1_scr = reg_v2_1.ssr\n",
    "reg_v2_1_r2 = reg_v2_1.rsquared\n",
    "\n",
    "display(Markdown(f\"\"\"\n",
    "Décomposition de la variance : <br>\n",
    "\n",
    "**SCT = SCE + SCR**\n",
    "- SCT = {reg_v2_1_sct:.0f} <br>\n",
    "- SCE = {reg_v2_1_sce:.0f} <br>\n",
    "- SCR = {reg_v2_1_scr:.0f} <br>\n",
    "\n",
    "$R^2 = {reg_v2_1_r2:.2f}$\n",
    "\n",
    "Le modèle explique **{reg_v2_1_r2*100:.0f}% de la variance totale**.\n",
    "\"\"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42702d40",
   "metadata": {},
   "source": [
    "Selon ce modèle, le pays de naissance (ie. le revenu moyen, la classe de revenus des aprents, et l’indice de Gini) explique donc **65%** de la variance totale, tandis que les autres facteurs non considérés dans le modèle (efforts, chance, etc.) représentent les **35%** restant. </br>\n",
    "Cependant, et comme pour le premier modèle, l'indice de gini n'est pas statistiquement significatif au seuil de 0,05."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5a55058",
   "metadata": {},
   "source": [
    "### Modèle 4 : revenus en log, gini & classes de revenus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cbc0b96",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Variables explicatives : revenu moyen du pays de l’individu en log, \n",
    " # l'indice de Gini du pays de l’individu, classe de revenu des parents\n",
    "reg_v2_2 = smf.ols('income_log ~ G_j+m_j_log+c_i_parent', data=data_cloned_c).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9e038e0",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "reg_v2_2_sct = reg_v2_2.centered_tss\n",
    "reg_v2_2_sce = reg_v2_2.ess\n",
    "reg_v2_2_scr = reg_v2_2.ssr\n",
    "reg_v2_2_r2 = reg_v2_2.rsquared\n",
    "\n",
    "display(Markdown(f\"\"\"\n",
    "Décomposition de la variance : <br>\n",
    "\n",
    "**SCT = SCE + SCR**\n",
    "- SCT = {reg_v2_2_sct:.0f} <br>\n",
    "- SCE = {reg_v2_2_sce:.0f} <br>\n",
    "- SCR = {reg_v2_2_scr:.0f} <br>\n",
    "\n",
    "$R^2 = {reg_v2_2_r2:.4f}$\n",
    "\n",
    "Le modèle explique **{reg_v2_2_r2*100:.0f}% de la variance totale**.\n",
    "\"\"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f937f05c",
   "metadata": {},
   "source": [
    "Le modèle le plus performant est celui expliquant le **revenu des enfants en log** par le **revenu du pays en log**, **l'indice de gini** et la **classe de revenus des parents**, avec un **$R^2 = 0,96$**.<br>\n",
    "On risque cependant de faire face à un problème de **sur-ajustement du modèle**, et une complication de l'interprétation du fait de l'utilisation des log."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40ad033e",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Le R^2 ajusté est similaire, avec une valeur de \", round(reg_v2_2.rsquared_adj, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96346f19",
   "metadata": {},
   "source": [
    "## Tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bded79b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Récapitulatif :\n",
    "R2 = [reg_v1_1_r2, reg_v1_2_r2,\n",
    "     reg_v2_1_r2, reg_v2_2_r2]#,reg_v3_1_r2, reg_v3_2_r2]\n",
    "reg_name = [\"Modèle 1\", \"Modèle 1_2\",\n",
    "           \"Modèle 2\", \"Modèle 2_2\"]#,\"Modèle 3\", \"Modèle 3_2\"]\n",
    "log = [\"Non\", \"Oui\", \"Non\", \"Oui\"]#, \"Oui\", \"Oui\"]\n",
    "var = ['y_child ~ G_j+m_j',\n",
    "       'income_log ~ G_j+m_j_log',\n",
    "       'y_child ~ G_j+m_j+c_i_parent',\n",
    "       'income_log ~ G_j+m_j_log+c_i_parent']#,'income_log ~ m_j_log + G_j','income_log ~ m_j_log + G_j+c_i_parent']\n",
    "list_of_tuples = list(zip(reg_name, var, R2, log))\n",
    "list_of_tuples \n",
    "\n",
    "df = pd.DataFrame(list_of_tuples,\n",
    "                  columns = ['Modèle', 'Variables', 'R2', 'Logarithme'])\n",
    "df.R2 = df.R2.round(2)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec77b3b2",
   "metadata": {},
   "source": [
    "On va conserver le modèle 2_2, utilisant à la fois les revenus, l'indice de Gini et la classe de revenu des parents, avec des variables non log, qui permettent une interprétation facilitée.\n",
    "On pourra également aborder le modèle."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f0b19e5",
   "metadata": {},
   "source": [
    "### Calcul des leviers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6a53eb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paramètres de l'étude\n",
    "n = data_cloned_c.shape[0] # échantillon\n",
    "p = 4 # nombre de variables\n",
    "\n",
    "# Seuil levier selon Belsey\n",
    "seuil_levier = 2 * p/n\n",
    "seuil_levier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56d69830",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ajout des leviers\n",
    "data_cloned_c['Leviers'] = reg_v2_2.get_influence().hat_matrix_diag\n",
    "\n",
    "Leviers_sup = data_cloned_c.loc[data_cloned_c['Leviers'] > seuil_levier, :]\n",
    "Leviers_inf = data_cloned_c.loc[data_cloned_c['Leviers'] <= seuil_levier, :]\n",
    "\n",
    "Leviers_sup.head().style.format({\"Leviers\": \"{:,.7f}\"}) # permet d'afficher pour la colonne choisie le nombre de chiffre apres la virgules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b09e40a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(Leviers_sup['country'].unique())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d17004b",
   "metadata": {},
   "source": [
    "On a donc 9 pays présentant des observations au dessus du seuil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ab5d94c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# Représentation des leviers\n",
    "plt.figure(figsize=(10,6))\n",
    "\n",
    "# Individus sous le seuil\n",
    "plt.bar(Leviers_inf['country'],Leviers_inf['Leviers'])\n",
    "\n",
    "# Individus au dessus du seuil\n",
    "plt.bar(Leviers_sup['country'],Leviers_sup['Leviers'])\n",
    "\n",
    "# Décoration et annotations\n",
    "plt.title('Représentation des leviers', fontsize=22)\n",
    "plt.xlabel('Pays', fontsize=18)\n",
    "plt.xticks('', fontsize=16)\n",
    "plt.ylabel('Leviers', fontsize=18)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.axhline(y=seuil_levier, linestyle='-')\n",
    "plt.text(50, 0.00085 , 'Seuil levier', fontsize = '14', color='red')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4ed00ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "country_lev = Leviers_sup['country'].unique()\n",
    "country_lev"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "712aa8d3",
   "metadata": {},
   "source": [
    "### Résidus studentisés : outliers"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b995d34d",
   "metadata": {},
   "source": [
    "Le seuil pour les résidus studentisés est une **loi de Student** à **n-p-1 degrés de liberté**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb6a3f54",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import t, shapiro\n",
    "alpha = 0.05\n",
    "data_cloned_c['rstudent'] = reg_v2_2.get_influence().resid_studentized_internal\n",
    "seuil_rstudent = t.ppf(1-alpha/2,n-p-1)\n",
    "seuil_rstudent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20ace274",
   "metadata": {},
   "outputs": [],
   "source": [
    "# statistique de test par observation\n",
    "data_cloned_c.sort_values(by='rstudent').head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fb7cdcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cloned_c.sort_values(by='rstudent').tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fadccaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "rstudent_in = data_cloned_c.loc[(data_cloned_c['rstudent'] <= seuil_rstudent) & (data_cloned_c['rstudent'] >= -seuil_rstudent), :]\n",
    "rstudent_sup = data_cloned_c.loc[(data_cloned_c['rstudent'] > seuil_rstudent), :]\n",
    "rstudent_inf = data_cloned_c.loc[(data_cloned_c['rstudent'] < -seuil_rstudent), :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d960eb06",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# Représentation des résidus studentisés\n",
    "plt.figure(figsize=(10,10))\n",
    "\n",
    "# Individus entre le seuil mini et le seuil maxi\n",
    "plt.bar(rstudent_in['country'],rstudent_in['rstudent'], color='orange')\n",
    "\n",
    "# Individus au dessus du seuil maxi\n",
    "plt.bar(rstudent_sup['country'],rstudent_sup['rstudent'], color='steelblue')\n",
    "\n",
    "# Individus en dessous du seuil maxi\n",
    "plt.bar(rstudent_inf['country'],rstudent_inf['rstudent'], color='steelblue')\n",
    "\n",
    "# Annotations\n",
    "plt.title('Représentation des résidus studentisés', fontsize=22)\n",
    "plt.xlabel('Pays', fontsize=18)\n",
    "plt.xticks('', fontsize=16)\n",
    "plt.ylim(-22,9)\n",
    "plt.ylabel('Résidus studentisés', fontsize=18)\n",
    "plt.yticks(fontsize=16)\n",
    "\n",
    "plt.axhline(y=seuil_rstudent, color='steelblue', linestyle='-')\n",
    "plt.text(50, 6.2 , 'Seuil rstudent', fontsize = '18', color='steelblue')\n",
    "\n",
    "plt.axhline(y=-seuil_rstudent, color='steelblue', linestyle='-')\n",
    "plt.text(50, -6.2 , '-Seuil rstudent', fontsize = '18', color='steelblue')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1b1fb7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Nombre de valeurs atypiques sur les variables à expliquer\n",
    "res_stu_ln = data_cloned_c.loc[(data_cloned_c['rstudent'] > seuil_rstudent) | (data_cloned_c['rstudent'] < -seuil_rstudent)]\n",
    "\n",
    "len(res_stu_ln)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b51b6d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pays presentants des valeurs atypiques\n",
    "outliers_country = _country = data_cloned_c.loc[data_cloned_c.index.isin(res_stu_ln.index)].groupby(by='country').count().sort_values(by='y_child',ascending=False)\n",
    "outliers_country "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "120f87d8",
   "metadata": {},
   "source": [
    "### Distance de cook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71c93a91",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Création du dataframe avec tous les resultats d'influences disponibles\n",
    "influence_ln = reg_v2_2.get_influence().summary_frame()\n",
    "influence_ln.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "302e3ae2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ajout de la colonne distance de Cook à notre dataframe analyses_ln\n",
    "data_cloned_c['dcooks'] = influence_ln['cooks_d']\n",
    "# Seuil d'influence selon Cook\n",
    "seuil_dcook = 4/(n-p)\n",
    "data_cloned_c.sort_values(by='dcooks').head().style.format({\"dcooks\": \"{:,.13f}\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f3daac6",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cloned_c.sort_values(by='dcooks').tail().style.format({\"dcooks\": \"{:,.7f}\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03da6210",
   "metadata": {},
   "outputs": [],
   "source": [
    "cooks_sup = data_cloned_c.loc[(data_cloned_c['dcooks'] > seuil_dcook), :]\n",
    "cooks_inf = data_cloned_c.loc[(data_cloned_c['dcooks'] <= seuil_dcook), :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2d5d9b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# Représentation de la distances de Cooks\n",
    "plt.figure(figsize=(10,10))\n",
    "\n",
    "# Individus sous le seuil\n",
    "plt.bar(cooks_inf['country'],cooks_inf['dcooks'], color='blue')\n",
    "\n",
    "# Individus au dessus le seuil\n",
    "plt.bar(cooks_sup['country'],cooks_sup['dcooks'], color='orange')\n",
    "\n",
    "# Décoration et annotations\n",
    "plt.title('Représentation de la distance de Cook', fontsize=22)\n",
    "plt.ylabel('Distance de Cook', fontsize=18)\n",
    "plt.xticks('', fontsize=16)\n",
    "plt.xticks(fontsize=16)\n",
    "\n",
    "plt.text(-5, (seuil_dcook + 0.00300), 'Seuil de Cook', fontsize = '18', color='steelblue')\n",
    "plt.axhline(y=seuil_dcook, color='steelblue', linestyle='-')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f4da8cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Nombre d'observations influentes\n",
    "dco_ln = data_cloned_c.loc[data_cloned_c['dcooks'] > seuil_dcook]\n",
    "len(dco_ln)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1512a875",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pays présentants des obervations au dela du seuil\n",
    "data_cloned_c.loc[data_cloned_c.index.isin(dco_ln.index)].groupby(by='country').count().sort_values(by='y_child',ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0006f18e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Obeservations atypiques et influentes\n",
    "ind_aty_infl_ln = data_cloned_c.loc[((data_cloned_c['dcooks'] > seuil_dcook) &\n",
    "                                  (data_cloned_c['rstudent'] > seuil_rstudent) &\n",
    "                                  (data_cloned_c['Leviers'] > seuil_levier)) |\n",
    "                                  ((data_cloned_c['dcooks'] > seuil_dcook) &\n",
    "                                  (data_cloned_c['rstudent'] < -seuil_rstudent) &\n",
    "                                  (data_cloned_c['Leviers'] > seuil_levier))]\n",
    "len(ind_aty_infl_ln)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15e05213",
   "metadata": {},
   "outputs": [],
   "source": [
    "susp_pt = data_cloned_c.loc[data_cloned_c.index.isin(ind_aty_infl_ln.index)]\n",
    "susp_pt.reset_index(inplace=True)\n",
    "susp_pt.groupby('country').count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc9aab43",
   "metadata": {},
   "outputs": [],
   "source": [
    "susp_pt.groupby('country').count().sort_values(by = \"index\",ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fb4656e",
   "metadata": {},
   "source": [
    "### Valeurs influentes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "002bf754",
   "metadata": {},
   "source": [
    "from statsmodels.graphics.regressionplots import *\n",
    "influence_plot(reg_v2_2) # myreg doit être un model de statsmodels\n",
    "\n",
    "#plt.xlim(0,0.002) # paramétrage manuel\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04a7467d",
   "metadata": {},
   "source": [
    "### Colinéarité des variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "828ceee6",
   "metadata": {},
   "outputs": [],
   "source": [
    "variables = reg_v2_2.model.exog\n",
    "[variance_inflation_factor(variables, i) for i in np.arange(1,variables.shape[1])]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56fd0028",
   "metadata": {},
   "source": [
    "### Homoscédasticité des résidus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b69b057",
   "metadata": {},
   "outputs": [],
   "source": [
    "_, pval, __, f_pval = sm.stats.diagnostic.het_breuschpagan(reg_v2_2.resid, variables)\n",
    "print('p value test Breusch Pagan:', pval)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b39c88a",
   "metadata": {},
   "source": [
    "### Normalité des résidus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1a12561",
   "metadata": {},
   "outputs": [],
   "source": [
    "shapiro(reg_v2_2.resid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0687cdef",
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.kstest(reg_v2_2.resid, 'norm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1c1d9e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.probplot(reg_v2_2.resid, dist=\"norm\", plot=pylab)\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2ae991d",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(reg_v2_2.resid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc514219",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Distribution des résidus\n",
    "fig, ax = plt.subplots(1, 2, figsize=(20,10))\n",
    "\n",
    "plt.hist(reg_v2_2.resid, density=True)\n",
    "\n",
    "model_norm_residuals = reg_v2_2.get_influence().resid_studentized_internal\n",
    "QQ = sm.ProbPlot(model_norm_residuals)\n",
    "QQ.qqplot(line='45', alpha=0.5, color='#4C72B0', ax=ax[0])\n",
    "\n",
    "ax[0].set_title('Q-Q Plot')\n",
    "ax[1].set_title('Histogramme des résidus')\n",
    "ax[1].set_xlabel('Valeurs résiduelles')\n",
    "ax[1].set_ylabel('Nombre de résidus')\n",
    "                \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46638451",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Homoscédasticité\n",
    "\n",
    "print(sm.stats.diagnostic.het_breuschpagan(reg_v2_2.resid, reg_v2_2.model.exog))\n",
    "\n",
    "ax=plt.plot(reg_v2_2.fittedvalues, reg_v2_2.resid, \".\",  alpha=0.3)\n",
    "plt.title(\"Nuage de la variance résiduelle\", fontsize=18)\n",
    "#plt.xlabel(\"GWh\", fontsize=16), plt.ylabel(\"Résidus\", fontsize=16)\n",
    "\n",
    "from statsmodels.stats.diagnostic import linear_rainbow\n",
    "Ftest, pval = linear_rainbow(reg_v2_2)\n",
    "print(pval)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e8522c3",
   "metadata": {},
   "source": [
    "### Prédiction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70100f06",
   "metadata": {},
   "outputs": [],
   "source": [
    "country = data_cloned_c['country'].unique()\n",
    "country"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77603f0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "pays_selectionné = 'France'\n",
    "c_i_parent_target = 50.0\n",
    "\n",
    "print('Le pays selectionné est :', pays_selectionné)\n",
    "\n",
    "revenu_moyen_pays = data_cloned_c.loc[(data_cloned_c['country'] == pays_selectionné), 'm_j'].iloc[0]\n",
    "print('Le revenu moyen du pays selectionné est de', round(revenu_moyen_pays, 2),\"$\")\n",
    "\n",
    "revenu_moyen_pays_log = np.log(revenu_moyen_pays)\n",
    "print('Le revenu moyen en log du pays selectionné est de', round(revenu_moyen_pays_log, 3))\n",
    "\n",
    "indice_gini_pays = data_cloned_c.loc[(data_cloned_c['country'] == pays_selectionné), 'G_j'].iloc[0]\n",
    "print(\"L'indice de gini du pays selectionné est de\", round(indice_gini_pays))\n",
    "\n",
    "a_prevoir = pd.DataFrame({'m_j_log':[revenu_moyen_pays_log], 'G_j':[indice_gini_pays],'c_i_parent':[c_i_parent_target]})\n",
    "Revenu_enfant = reg_v2_2.predict(a_prevoir)\n",
    "\n",
    "# Ajouter * devant une variable lors d'un print permet de ne pas afficher le Dtypes ni l'index\n",
    "Revenu_enfant_calculé = np.exp(*Revenu_enfant)\n",
    "Revenu_enfant_calculé\n",
    "\n",
    "print(\"Le revenu d'un individu en\",pays_selectionné,\n",
    "      \"dont les parents ont une classe de revenu égale à\",c_i_parent_target,\n",
    "      \"est de\",round(Revenu_enfant_calculé, 2),\"$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0b456dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = data_cloned_c.loc[(data_cloned_c['country'] == pays_selectionné)]\n",
    "g = f.loc[f[\"c_i_parent\"] == c_i_parent_target]\n",
    "h = g.y_child.mean()\n",
    "print(round(h, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f5a34ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "%whos DataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "668ced76",
   "metadata": {},
   "outputs": [],
   "source": [
    "del data_cloned\n",
    "del data_cloned_b\n",
    "del data_cloned_c\n",
    "del WID\n",
    "del Leviers_inf\n",
    "del dist_revenus_ok\n",
    "del echantillon_m4\n",
    "del influence_ln\n",
    "del rstudent_in\n",
    "del sample\n",
    "del sample_2\n",
    "del sample_3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94990197",
   "metadata": {},
   "source": [
    "# FIN ICI"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}