--- title: Essais et développement de la méthode de calage sur marges keywords: fastai sidebar: home_sidebar nb_path: "notebooks/calmar/calmar.ipynb" ---
import unittest
tc = unittest.TestCase()
import pandas as pd
import seaborn as sns
from leximpact_socio_fisca_simu_etat.config import Configuration
from matplotlib import pyplot as plt
config = Configuration(project_folder="leximpact-prepare-data")
erfs_03 = pd.read_hdf(
config.get("DATA_OUT") + "03_erfs_rfr_cal_ind" + config.get("YEAR_ERFS") + ".h5"
)
erfs_03.columns
erfs_03.tail()
dataset = erfs_03[
[
"idfoy",
"quifoy",
"wprm",
"salaire_de_base",
"chomage_brut",
"retraite_brute",
"f4ba",
"rfr",
]
]
# Pour tester sur une partie de la base
# dataset = dataset.sample(n=1000, random_state=53)
print("On a ", len(dataset), "individus dans cet échantillon")
# Ici, wprm sont les poids initiaux dits "de sondage"
dataset.rename(columns={"wprm": "wprm_init"}, inplace=True)
wprm_init = dataset["wprm_init"]
dataset.head()
margins_avant = {}
margins_avant.update(
{
"salaire_de_base": (dataset["wprm_init"] * dataset["salaire_de_base"]).sum(),
"chomage_brut": (dataset["wprm_init"] * dataset["chomage_brut"]).sum(),
"retraite_brute": (dataset["wprm_init"] * dataset["retraite_brute"]).sum(),
"f4ba": (dataset["wprm_init"] * dataset["f4ba"]).sum(),
}
)
margins_avant
tc.assertEqual(dataset.isna().sum().sum(), 0)
# Si c'est le cas, on pourra les enlever:
# dataset = dataset.fillna(0) ou dataset = dataset.dropna(0)
sns.pairplot(
dataset[["salaire_de_base", "chomage_brut", "retraite_brute", "rfr"]],
diag_kind="kde",
)
TODO: voir si on fait une normalization
margins = {}
margins.update(
{
"salaire_de_base": 650855163531.0, # Salaire imposable POTE 2019
"chomage_brut": 0.97
* (
dataset["wprm_init"] * dataset["chomage_brut"]
).sum(), # On ne connait pas le chiffre, donc on met un ecart pour calibrer
"retraite_brute": 307254581479.0, # Retraite POTE 2019
"f4ba": 1.34 * (dataset["wprm_init"] * dataset["f4ba"]).sum(),
}
) # On ne connait pas le chiffre, donc on met un ecart pour calibrer
# On ajoute une limite de population
margins.update({"total_population": 38549926})
margins
Calibrates weights to satisfy margins constraints
Args:
total_population with value the total population. If absent it is initialized to the actual total populationpondfin_out, lambdasol, margins_new_dict = calmar(
dataset,
margins,
initial_weight="wprm_init",
method="raking ratio",
lo=0.33,
up=3,
use_proportions=False,
xtol=1.49012e-08,
maxfev=256,
)
# pondfin_out
# lambdasol
# margins_new_dict
variables = ["salaire_de_base", "chomage_brut", "retraite_brute"]
def compute_dataset_pondere(dataset, wprm, type_cal):
data_pond = dataset.copy()
for variable in dataset.columns:
data_pond[variable] = data_pond[variable] * wprm
data_pond["wprm"] = wprm
data_pond["type"] = type_cal
return data_pond
def analyse_calage(dataset, pondfin_out):
# Resultat sur les marges
check_calmar(
dataset,
margins,
initial_weight="wprm_init",
pondfin_out=pondfin_out,
lambdasol=lambdasol,
margins_new_dict=margins_new_dict,
)
# Données pondérées après calage
data_cal_pond = compute_dataset_pondere(dataset, pondfin_out, "cal")
# Analyse
sns.pairplot(
data_cal_pond,
diag_kind="kde",
)
plt.show()
analyse_calage(dataset[variables], pondfin_out)
Calibrates weights to satisfy margins constraints
Args:
total_population with value the total population. If absent it is initialized to the actual total populationpondfin_lin, lambdasol, margins_new_dict = calmar(
dataset,
margins,
initial_weight="wprm_init",
method="linear",
lo=0,
up=3,
use_proportions=False,
xtol=1.49012e-08,
maxfev=256,
)
analyse_calage(dataset[variables], pondfin_lin)
pondfin_rr, lambdasol, margins_new_dict = calmar(
dataset,
margins,
initial_weight="wprm_init",
method="raking ratio",
lo=0,
up=3,
use_proportions=False,
xtol=1.49012e-08,
maxfev=256,
)
analyse_calage(dataset[variables], pondfin_rr)
pondfin_log, lambdasol, margins_new_dict = calmar(
dataset,
margins,
initial_weight="wprm_init",
method="logit",
lo=0,
up=3,
use_proportions=False,
xtol=1.49012e-08,
maxfev=256,
)
analyse_calage(dataset[variables], pondfin_log)
def matrice_violin(dataset, results):
# Initial
data_init_pond = compute_dataset_pondere(
dataset[variables], dataset["wprm_init"], "init"
)
matrix = data_init_pond.copy()
# Pour chaque calibration on concatene dans une matrice
for cal in results:
data_cal_pond = compute_dataset_pondere(dataset[variables], results[cal], cal)
matrix = pd.concat([matrix, data_cal_pond]) # .reset_index()
return matrix
results = {}
results.update(
{"linear": pondfin_lin, "raking_ratio": pondfin_rr, "logit": pondfin_log}
)
len(dataset)
matrix = matrice_violin(dataset, results)
len(matrix)
matrix.head()
matrix.tail()
for var in variables:
sns.violinplot(data=dataset, y=var)
plt.show()
for var in variables:
ax = sns.violinplot(data=matrix, y=var, x="type", split=True)
plt.show()
On voit bien qu'elle détonne. Je suis étonnée des résultats négatifs...
min(results["linear"])
results2 = results.copy()
results2.pop("linear")
results2
matrix2 = matrice_violin(dataset, results2)
# On trace les distributions APRES
for var in variables:
ax = sns.violinplot(data=matrix2, y=var, x="type", split=True)
plt.show()
margins
for var in variables:
print(var)
tot_rr = (dataset[var] * results2["raking_ratio"]).sum()
tot_log = (dataset[var] * results2["logit"]).sum()
error_rr = abs(tot_rr - margins[var]) / margins[var]
error_log = abs(tot_log - margins[var]) / margins[var]
print("Total initial", (dataset[var] * dataset["wprm_init"]).sum())
print("Total Raking Ratio", tot_rr, ", erreur : ", error_rr)
print("Total Logit", tot_log, ", erreur : ", error_log, "\n")
tot_rr = (dataset["rfr"] * results2["raking_ratio"]).sum()
tot_log = (dataset["rfr"] * results2["logit"]).sum()
print("Total initial", (dataset["rfr"] * dataset["wprm_init"]).sum())
print("Total Raking Ratio", tot_rr)
print("Total Logit", tot_log)
Ici on met en place pour tracer les distributions pondérées car les précédentes ne représentent rien:
[x, wprm] = [5, 10] -> aujourd'hui on plot une personne de salaire 50, alors qu'on devrait plotter 10 personnes de salaire 5 !
Pour représenter ça on a 2 options:
1 - On génère une matrice avec 'wprm_k' personnes pour chaque ligne k (précis)
2 - On fait des histogrames en adaptant la population totale (rapide)
dataset.head()
# results.pop("linear")
def calmar_hist(dataset, variables, wprm_init, results, min_=None, max_=None):
plot_dict = {}
for var in variables:
if min_ is None:
min_ = 0
if max_ is None:
max_ = max(dataset[var])
fig, ax = plt.subplots(figsize=(15, 5))
plt.hist(
dataset[var],
bins=100,
weights=wprm_init,
range=(min_, max_),
histtype="step",
alpha=0.5,
linewidth=2,
)
legend = ["Initial"]
for method in results.keys():
plt.hist(
dataset[var],
bins=100,
weights=results.get(method),
range=(min_, max_),
histtype="step",
alpha=0.5,
)
legend.append(method)
# Mise en forme
plt.legend(legend)
plt.title("Comparaison des méthodes de calage sur marges pour le " + var)
plt.xlabel(var)
plt.ylabel("Total de population")
plot_dict.update({f"Calage_{var}": fig})
plt.savefig(f"Calage_{var}.png")
return plot_dict
plot_dict = calmar_hist(dataset, variables, wprm_init, results)
plot_dict_cent = calmar_hist(dataset, variables, wprm_init, results, min_=500_000)
plot_dict_hr = calmar_hist(dataset, variables, wprm_init, results, min_=2_000_000)
def normalizer(dataset, variables):
norm = {}
for var in variables:
mean_ = dataset[var].mean()
std_ = dataset[var].std()
dataset[var] = (dataset[var] - mean_) / std_
norm.update({f"norm_{var}": [mean_, std_]})
print(var, dataset[var].min(), dataset[var].max())
return norm, dataset
norm, datanorm = normalizer(dataset, variables)
plot_dict = calmar_hist(datanorm, variables, wprm_init, results, min_=0)
plot_dict = calmar_hist(datanorm, variables, wprm_init, results, min_=90)
var_to_check = "salaire_de_base"
analyse = dataset[[var_to_check, "wprm_init"]]
analyse["logit"] = pondfin_log
analyse["raking_ration"] = pondfin_rr
analyse["linear"] = pondfin_lin
analyse.sort_values(by=var_to_check, inplace=True)
analyse
dataset
Next steps:
Checks calmar:
Calage ERFS 2018 - POTE 2018 -> attention, ici on tire forcement la base dans une direction Vieillissement 2019 -> est-ce que c'est pertinent ? Insertion donnees Rk 2019 Calage avec POTE 2019 -> ca permet de tester la pertinence de nos procedures de vieillissement Vieillissement -> 2021
Vieillissement: il faut adapter la distribution: parfois on a une nouvelle distribution (on adapte les var de l'inflation, le total de pop, mais aussi le nb de gens qui ont/n'ont pas le chomage/la retraite) -> on peut ajouter une cible (=indicatrice) "au chomage" en binaire pour caler sur cette nouvelle marge
Separer l'exercice:
Etat de l'art (IPP)
Pour faire du calage en utilisant les buckets de POTE:
TO DO