Statistická analýza determinantů světového štěstí¶
V tomto projektu budeme analyzovat datasety World Happiness Report a Countries of the World. Budeme hledat souvislosti mezi štěstím populace státu a různými ukazateli jako je gramotnost, hrubý domácí produkt a další. Ze získaných poznatků bude možné zjistit, na jaké ukazatele je třeba se zaměřit, pro zvýšení štěstí obyvatelstva.
Plán analýzy¶
První se seznámíme s daty a provedeme kontroly, že data rámcově obsahují to, co bychom čekali (sanity-check). Následně data pročistíme a připravíme na analýzu. Poté vytvoříme základní model, který bude predikovat štěstí populace na základě gramotnosti populace. Očekáváme, že lidé jsou díky své gramotnosti schopni tvořit lépe fungující společnost a díky tomu žít šťastnější životy. Dále vytvoříme model, který bude dělat predikce na základě více proměnných. Na závěr porovnáme vytvořené modely mezi sebou. Pro samotnou analýzu dat použijeme technologie Python, matplotlib, statsmodels, scikit-learn a další.
Data¶
Dataset týkající se štěstí populace státu obsahuje data z let 2015 až 2019 a je dostupný ve formátu CSV. Pro jednoduchost budeme analýzu provádět pouze pro rok 2019, která mají velikost 9 KB a obsahuje data o 157 zemích. Z dat z roku 2019 použijeme následující sloupce:
- Country or region - země či region, jichž se řádek týká
- Score - Měření subjektivní pohody (SWB) pochází z vydání Gallup World Poll (GWP) z 26. února 2021 a pokrývá roky 2005 až 2020. Pokud není uvedeno jinak, jde o národní průměr odpovědí na otázku týkající se hodnocení života. Překlad znění otázky do Českého jazyka je následující: „Představte si žebřík, jehož stupně jsou očíslovány od 0 na spodním konci po 10 na horním konci. Vrchol žebříku představuje nejlepší možný život pro vás a spodní část žebříku nejhorší možný život pro vás. Na kterém stupni žebříku byste řekli, že se osobně nacházíte v tuto chvíli?“ Toto měření je také označováno jako Cantrilův žebřík života.
Dataset týkající se metrik zemí světa je složen z několika veřejných datasetů a je dostupný ve formátu CSV. Velikost datasetu je 38 KB a popisuje metriky pro 228 zemí. Sloupce obsahující základní geografická data jako je populace, či rozloha obsahují hodnoty pro všechny země v tabulce. Ostatní metriky jako gramotnost nejsou dostupné pro všechny země v tabulce, ale vždy chybí pouze několik jednotek hodnot. Chybějící hodnoty by neměly zásadně ovlivnit naši analýzu. Popisy jednotlivých sloupců bylo potřeba najít v dokumentací zdrojových datasetů. Mnoho z těchto datasetů bohužel popisy sloupců neobsahovaly a proto je popis nahrazen pouhým překladem názvu sloupce. Nedostupnost popisu jednotlivých sloupců značně komplikuje analýzu a snížuje důvěryhodnost datasetu. Následuje výčet sloupců, které budou použity pro naši analýzu:
- Country – Země
- Birthrate – Porodnost
- Infant mortality (per 1000 births) – Úmrtnost kojenců (na 1000 narozených)
- GDP ($ per capita) – HDP (v dolarech na obyvatele)
- Literacy (%) – Gramotnost (%)
- Phones (per 1000) – Telefony (na 1000 obyvatel)
- Agriculture – Zemědělství
Kompletní popis sloupců použitých datasetů je dostupný v příloze na konci dokumentu.
Čištění dat¶
V této podsekci se budeme zabývat čištěním a normalizací dat.
V datasetu s metrikami zemí ve sloupci "Country" končí každá hodnota mezerou. To znemožňuje sloučení s druhým datasetem na základě názvu země. Proto je potřeba z těchto hodnot smazat bílé znaky z okrajů řetězců. Po sloučení tabulek na základě názvu země tímto způsobem bychom získali ve výsledném tabulce 142 řádků. Oba datasety mají však více řádků, čímž se dostáváme k dalšímu problému.
Užité datsety používají různé identifikátory států, např. dataset o štěstí
populace používá název "Trinidad and Tobago" zatímco dataset s metrikami zemí
použá pro stejný stát název "Trinidad & Tobago". Budeme potřebovat
sofistikovanější způsob nežli pouhé porovnání na rovnost. Před samotným
sloučením v obou datasetech sjednotíme názvy zemí použítím Python balíčku
pycountry
, který umí dělat fuzzy hledání v databázi zemí.
import numpy as np
import pandas as pd
import pycountry
def standardize_country_name(country_name: str):
"""
Try fuzzy matching against a DB of country names.
Otherwise at least trim whitespace.
"""
stripped_country_name = country_name.strip()
try:
return pycountry.countries.search_fuzzy(stripped_country_name)[0].name
except LookupError:
return stripped_country_name
countries_of_the_world = pd.read_csv(
"./data/countries_of_the_world/countries_of_the_world.csv",
dtype={"Literacy (%)": np.float32},
decimal=",",
)
countries_of_the_world["Country"] = countries_of_the_world["Country"].apply(
standardize_country_name
)
happiness = pd.read_csv(
"./data/world_happiness/2019.csv", dtype={"Score": np.float32}
)
happiness["Country or region"] = happiness["Country or region"].apply(
standardize_country_name
)
Sloučení tabulek¶
Pročišěné tabulky nyní sloučíme a ověříme, že velikost výsledné tabulky je dostatečná.
print('Velikost tabulky "Countries of the world":', len(countries_of_the_world))
print('Velikost tabulky "World happiness":', len(happiness))
print()
merged = pd.merge(
countries_of_the_world,
happiness[["Score", "Country or region"]],
how="inner",
left_on="Country",
right_on="Country or region",
)
print("Velikost sloučené tabulky: ", len(merged))
Velikost tabulky "Countries of the world": 227 Velikost tabulky "World happiness": 156 Velikost sloučené tabulky: 146
Až na pár jednotek řádků se podařilo tabulky sloučit, což bude pro naši analýzu dostatečné.
Explorační analýza¶
Datasety jsou poměrně velké, takže začneme výběrem sloupců, které mohou být užity pro predikci štěstí populace. Vybereme 5 sloupců s největší korelací se skórem štěstí a s těmi budeme dále pracovat. Pro výpočet korelace použijeme Pearsonův korelační koeficient. Pro jednoduchost budeme počítat korelaci pouze se sloupci, která obsahují numerické hodnoty. Není pro nás důležité zda-li je hodnota korelace pozitivní či negativní, proto budeme pro seřazení uvažovat absolutní hodnotu korelace.
from scipy.stats import pearsonr
# Column to compare with all others
target_column = "Score"
numeric_columns = merged.select_dtypes(include=["number"]).columns
# Calculate correlations
correlations = {}
for column in numeric_columns:
if column != target_column:
valid_data = merged[[target_column, column]].dropna()
corr, _ = pearsonr(valid_data[target_column], valid_data[column])
correlations[column] = corr
# Display results as a DataFrame
corr_col_name = "Correlation with " + target_column
col_name = "Column name"
correlations_df = pd.DataFrame(
correlations.items(), columns=[col_name, corr_col_name]
)
abs_corr_col_name = "Absolute Correlation"
correlations_df[abs_corr_col_name] = correlations_df[corr_col_name].abs()
correlations_df.sort_values(by=abs_corr_col_name, ascending=False, inplace=True)
correlations_df.drop(columns=[abs_corr_col_name], inplace=True)
# Get top n results
top_n = 5
top_correlations = correlations_df.head(top_n)
# Store columns with highest correlation in a variable for later use
top_corr_columns = top_correlations[col_name].to_list()
# Display result without row ID
top_correlations.style.hide()
Column name | Correlation with Score |
---|---|
Phones (per 1000) | 0.753111 |
GDP ($ per capita) | 0.750321 |
Infant mortality (per 1000 births) | -0.700899 |
Birthrate | -0.658775 |
Agriculture | -0.641866 |
Nyní už máme vybrané sloupce, podle které můžeme využít k predikci. Překvapivě zde není gramotnost populace, kterou budeme používat pro tvorbu referenčního modelu.
Nyní se podíváme na histogramy skóre štěstí (predikovaná proměnná), gramotnosti (prediktor v referenčním modelu) a ukazatelů s vysokou korelací (prediktory ve vylepšeném modelu), abychom se s daty seznámili.
%matplotlib inline
import matplotlib.pyplot as plt
import math
import pandas as pd
def plot_histograms(df: pd.DataFrame, colors: dict = {}):
"""
Plot histograms for each column in the provided DataFrame.
colors - dictionary with column names as keys and colors as values
"""
fig, axs = plt.subplots(math.ceil(len(df.columns) / 2), 2, figsize=(10, 10))
for col, column_name in enumerate(df.columns):
row = col // 2
col = col % 2
axs[row, col].hist(
df[column_name],
bins=15,
edgecolor="black",
color=colors.get(column_name),
)
axs[row, col].set(xlabel="Hodnota", ylabel="Frekvence")
axs[row, col].set_title(column_name)
fig.delaxes(axs[3, 1]) # remove unused subplot
plt.tight_layout()
plt.show()
reference_model_predictor_column_name = "Literacy (%)"
predicted_variable_column_name = "Score"
data = merged[
top_corr_columns
+ [reference_model_predictor_column_name, predicted_variable_column_name]
]
plot_histograms(data)
Data vypadají tak, jak bychom očekávali. Žádné zvláštnosti z histogramů nelze vidět. Nyní se podívejme jak data těchto sloupců pokrývají země světa. Pro vykreslení mapy použijeme veřejně dostupný dataset od World Bank Group. Státy vykresleny zelenou barvu obsahují data pro všechny sloupce, které budou použity pro vytvoření modelů. Státy pro, které tato podmínka neplatí, jsou vykresleny červeně.
import geopandas as gpd
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
world_map = gpd.read_file("zip://data/countries.zip!WB_countries_Admin0_10m")
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot()
data_for_model = merged[
top_corr_columns
+ [
"Country",
reference_model_predictor_column_name,
predicted_variable_column_name,
]
].dropna()
# add color for each country based on coverage
color_col_name = "color"
def assign_color(name: str):
if (
name is not None
and data_for_model["Country"].eq(standardize_country_name(name)).any()
):
return "green"
return "red"
world_map[color_col_name] = world_map["FORMAL_EN"].map(assign_color)
world_map.plot(
ax=ax,
color=world_map[color_col_name],
edgecolor="black",
alpha=0.5,
legend=True,
categorical=True,
missing_kwds={
"color": "green",
"label": "Covered",
},
)
# setup legend
color_dict = {"Pokryto": "green", "Nepokryto": "red"}
custom_points = [
Line2D([0], [0], marker="o", linestyle="none", markersize=10, color=color)
for color in color_dict.values()
]
leg_points = ax.legend(custom_points, color_dict.keys())
# turn off axis ticks
ax.set_xticks([])
ax.set_yticks([])
plt.title("Pokrytí zemí světa")
plt.show()
Jednotlivé kontinenty jsou pokryté mírně nerovnoměrně. Špatně je pokrytá zejména Afrika a Austrálie. Překvapivě dobře je pokryta Asie a Jižní Amerika. Data z mnoha chudších zemí sice máme, ale statistické úřady těchto zemí lze jen těžko považovat za spolehlivé. Náš model bude tedy ovlivněn nepřesností či neúplností některých dat. Celkově bychom však měli mít dostatečně kompletní a kvalitní data pro vytvoření funkčního modelu.
Transformace dat¶
Na histogramech, které byly vytvořeny v rámci explorační analýzy, jsme si mohli všimnout, že některé sloupce vypadají, že mají log-normální rozdělení. Pro lepší predikci štěstí populace může být vhodné tyto sloupce převést na logaritmickou škálu. V případě gramotnosti je potřeba před aplikovaním logaritmu otočit hodnoty tak, aby byly větší hodnoty gramotnosti blízko 0 a naopak menší hodnoty gramotnosti blízko 100. Abychom předešli aplikaci logaritmu na nulu, přičteme k hodnotám gramotnosti 1. Tímto způsobem se nám podaří převést hodnoty gramotnosti na logaritmickou škálu. Na některé sloupce tuto transformaci nebudeme aplikovat, jelikož nemají vhodné rozdělení (např. sloupce Agriculture). Na histogramech níže můžeme vidět transformovaná data. Sloupce, které byly transformovány mají zelenou barvu. Sloupce, které byly ponechány v původním stavu mají modrou modrou barvu. Transformovaná data by nyní měla mít normální rozdělení.
transformed_data = merged[top_corr_columns + ["Literacy (%)", "Score"]]
# disable false positive warning
pd.options.mode.chained_assignment = None
# For these columns we apply the default log transformation
column_names_default_transformation = [
"Birthrate",
"GDP ($ per capita)",
"Phones (per 1000)",
"Infant mortality (per 1000 births)",
]
for column_name in column_names_default_transformation:
transformed_data[column_name] = transformed_data[column_name].map(
lambda x: math.log(x)
)
# For literacy we apply a custom transformation
transformed_data["Literacy (%)"] = (
transformed_data["Literacy (%)"]
.map(lambda x: 100 - x)
.map(lambda x: math.log(x + 1))
)
pd.options.mode.chained_assignment = "warn"
colorMap = {
column_name: "green"
for column_name in column_names_default_transformation + ["Literacy (%)"]
}
plot_histograms(transformed_data, colors=colorMap)
Modelování dat¶
První vytvoříme referenční model využívající gramotnost populace. Zobrazme si data v grafu pro nastavení očekávání výkonnosti referečního modelu.
import matplotlib.pyplot as plt
# if happiness score or literacy is unknown, ignore country
filtered = merged[
[reference_model_predictor_column_name, predicted_variable_column_name]
].dropna()
print("Velikost vzorku: ", len(filtered))
plt.scatter(
filtered[reference_model_predictor_column_name],
filtered[predicted_variable_column_name],
)
plt.xlabel("Gramotnost (%)")
plt.ylabel("Skóre štěstí")
plt.xlim(-5, 105)
plt.show()
Velikost vzorku: 144
Na první pohled nevypadá vztah úplně silně, ale na levé polovině grafu (nízká gramotnost) jsou hodnoty skóre menší než 6 a na pravé straně grafu (vysoká gramotnost) jsou hodnoty štěstí i nad 7. Pro vytvoření modelu využijeme lineární regresi.
import statsmodels.api as sm
reference_model_data = merged[
[reference_model_predictor_column_name, predicted_variable_column_name]
].dropna()
ref_predictors = reference_model_data[reference_model_predictor_column_name]
ref_predicted_var = reference_model_data[predicted_variable_column_name]
ref_predictors_with_constant = sm.add_constant(ref_predictors)
ref_model = sm.OLS(ref_predicted_var, ref_predictors_with_constant).fit()
ref_model.summary(title="Referenční model s netransformovanými daty")
Dep. Variable: | Score | R-squared: | 0.350 |
---|---|---|---|
Model: | OLS | Adj. R-squared: | 0.345 |
Method: | Least Squares | F-statistic: | 76.32 |
Date: | Tue, 17 Dec 2024 | Prob (F-statistic): | 6.09e-15 |
Time: | 10:40:31 | Log-Likelihood: | -187.59 |
No. Observations: | 144 | AIC: | 379.2 |
Df Residuals: | 142 | BIC: | 385.1 |
Df Model: | 1 | ||
Covariance Type: | nonrobust |
coef | std err | t | P>|t| | [0.025 | 0.975] | |
---|---|---|---|---|---|---|
const | 2.9475 | 0.297 | 9.933 | 0.000 | 2.361 | 3.534 |
Literacy (%) | 0.0309 | 0.004 | 8.736 | 0.000 | 0.024 | 0.038 |
Omnibus: | 3.296 | Durbin-Watson: | 1.771 |
---|---|---|---|
Prob(Omnibus): | 0.192 | Jarque-Bera (JB): | 2.507 |
Skew: | -0.180 | Prob(JB): | 0.286 |
Kurtosis: | 2.463 | Cond. No. | 333. |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
y_pred = ref_model.predict(ref_predictors_with_constant)
plt.figure(figsize=(8, 6))
plt.scatter(
ref_predictors, ref_predicted_var, color="blue", label="Datové body"
)
plt.plot(ref_predictors, y_pred, color="red", label="Předpověď")
plt.xlabel("Gramotnost (%)")
plt.ylabel("Skóre štěstí")
plt.title("Referenční model")
plt.legend()
plt.show()
Na grafu výše můžeme vidět jak se nám povedlo napasovat přímku na data. Pro vyhodnocení přesnosti predikcí modelu můžeme použít mnoho metrik. Některé z nich lze vidět ve vypsaném shrnutí vlastností modelu. Pro naše porovnávání modelů použijeme $R^2$, která je v případě našeho modelu zhruba $0.35$. Z této hodnoty vyplývá, že pomocí gramotnosti populace jsme schopni alespoň částečně odhadnout štěstí populace. Zkusme referenční model vylepšit tím, že použijeme transformovaná data pomocí logaritmu.
reference_model_transformed_data = transformed_data[
[reference_model_predictor_column_name, predicted_variable_column_name]
].dropna()
ref_trans_predictors = reference_model_transformed_data[
reference_model_predictor_column_name
]
ref_trans_predicted_var = reference_model_data[predicted_variable_column_name]
ref_trans_predictors_with_constant = sm.add_constant(ref_trans_predictors)
ref_trans_model = sm.OLS(
ref_predicted_var, ref_trans_predictors_with_constant
).fit()
ref_trans_model.summary(title="Referenční model s transformovanými daty")
Dep. Variable: | Score | R-squared: | 0.439 |
---|---|---|---|
Model: | OLS | Adj. R-squared: | 0.435 |
Method: | Least Squares | F-statistic: | 111.3 |
Date: | Tue, 17 Dec 2024 | Prob (F-statistic): | 1.44e-19 |
Time: | 10:40:31 | Log-Likelihood: | -176.90 |
No. Observations: | 144 | AIC: | 357.8 |
Df Residuals: | 142 | BIC: | 363.7 |
Df Model: | 1 | ||
Covariance Type: | nonrobust |
coef | std err | t | P>|t| | [0.025 | 0.975] | |
---|---|---|---|---|---|---|
const | 6.7154 | 0.138 | 48.647 | 0.000 | 6.443 | 6.988 |
Literacy (%) | -0.5533 | 0.052 | -10.548 | 0.000 | -0.657 | -0.450 |
Omnibus: | 3.053 | Durbin-Watson: | 2.060 |
---|---|---|---|
Prob(Omnibus): | 0.217 | Jarque-Bera (JB): | 2.562 |
Skew: | -0.215 | Prob(JB): | 0.278 |
Kurtosis: | 2.509 | Cond. No. | 5.82 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Dostali jsme mírné zlepšení. Hodnota $R^2$ se zvýšila na zhruba $0.44$. Nyní vytvoříme vylepšený model použitím dat s vysokými hodnotami korelace s predikovanou proměnnou. Hned poté zkusíme vylepšený model natrénovat na transformovaných datech. Nakonec modely porovnáme.
import statsmodels.api as sm
improved_model_data = merged[
top_corr_columns + [predicted_variable_column_name]
].dropna()
imp_predictors = improved_model_data[top_corr_columns]
imp_predicted_var = improved_model_data[predicted_variable_column_name]
imp_predictors_with_constant = sm.add_constant(imp_predictors)
imp_model = sm.OLS(imp_predicted_var, imp_predictors_with_constant).fit()
imp_model.summary(title="Vylepšený model s netransformovanými daty")
Dep. Variable: | Score | R-squared: | 0.652 |
---|---|---|---|
Model: | OLS | Adj. R-squared: | 0.639 |
Method: | Least Squares | F-statistic: | 51.25 |
Date: | Tue, 17 Dec 2024 | Prob (F-statistic): | 9.99e-30 |
Time: | 10:40:31 | Log-Likelihood: | -142.06 |
No. Observations: | 143 | AIC: | 296.1 |
Df Residuals: | 137 | BIC: | 313.9 |
Df Model: | 5 | ||
Covariance Type: | nonrobust |
coef | std err | t | P>|t| | [0.025 | 0.975] | |
---|---|---|---|---|---|---|
const | 5.1972 | 0.243 | 21.357 | 0.000 | 4.716 | 5.678 |
Phones (per 1000) | 0.0009 | 0.001 | 1.226 | 0.222 | -0.001 | 0.002 |
GDP ($ per capita) | 4.195e-05 | 1.23e-05 | 3.409 | 0.001 | 1.76e-05 | 6.63e-05 |
Infant mortality (per 1000 births) | -0.0098 | 0.004 | -2.681 | 0.008 | -0.017 | -0.003 |
Birthrate | 0.0056 | 0.011 | 0.500 | 0.618 | -0.016 | 0.028 |
Agriculture | -0.4193 | 0.706 | -0.594 | 0.553 | -1.815 | 0.976 |
Omnibus: | 2.764 | Durbin-Watson: | 2.194 |
---|---|---|---|
Prob(Omnibus): | 0.251 | Jarque-Bera (JB): | 2.810 |
Skew: | -0.319 | Prob(JB): | 0.245 |
Kurtosis: | 2.746 | Cond. No. | 1.81e+05 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.81e+05. This might indicate that there are
strong multicollinearity or other numerical problems.
improved_model_transformed_data = transformed_data[
top_corr_columns + [predicted_variable_column_name]
].dropna()
imp_trans_predictors = improved_model_transformed_data[top_corr_columns]
imp_trans_predicted_var = improved_model_transformed_data[
predicted_variable_column_name
]
imp_trans_predictors_with_constant = sm.add_constant(imp_trans_predictors)
imp_trans_model = sm.OLS(
imp_trans_predicted_var, imp_trans_predictors_with_constant
).fit()
imp_trans_model.summary(title="Vylepšený model s transformovanými daty")
Dep. Variable: | Score | R-squared: | 0.681 |
---|---|---|---|
Model: | OLS | Adj. R-squared: | 0.670 |
Method: | Least Squares | F-statistic: | 58.59 |
Date: | Tue, 17 Dec 2024 | Prob (F-statistic): | 2.36e-32 |
Time: | 10:40:31 | Log-Likelihood: | -135.68 |
No. Observations: | 143 | AIC: | 283.4 |
Df Residuals: | 137 | BIC: | 301.1 |
Df Model: | 5 | ||
Covariance Type: | nonrobust |
coef | std err | t | P>|t| | [0.025 | 0.975] | |
---|---|---|---|---|---|---|
const | -0.3342 | 1.640 | -0.204 | 0.839 | -3.578 | 2.909 |
Phones (per 1000) | 0.2057 | 0.088 | 2.347 | 0.020 | 0.032 | 0.379 |
GDP ($ per capita) | 0.4417 | 0.131 | 3.374 | 0.001 | 0.183 | 0.700 |
Infant mortality (per 1000 births) | -0.5139 | 0.133 | -3.876 | 0.000 | -0.776 | -0.252 |
Birthrate | 0.8711 | 0.249 | 3.497 | 0.001 | 0.379 | 1.364 |
Agriculture | 0.8917 | 0.719 | 1.240 | 0.217 | -0.530 | 2.313 |
Omnibus: | 6.481 | Durbin-Watson: | 2.105 |
---|---|---|---|
Prob(Omnibus): | 0.039 | Jarque-Bera (JB): | 6.057 |
Skew: | -0.475 | Prob(JB): | 0.0484 |
Kurtosis: | 3.337 | Cond. No. | 337. |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
$R^2$ se v obou případech zvýšilo. V případě netranformovaných dat na $0.65$, takže jsme dostali lepší model. V případě transformovaných dat máme $R^2 \doteq 0.68$. Transformace dat pro model pracující s více proměnnými velké zlepšení nepřinesla. Pro lineární regresi pracující s více proměnnými se pro ohodnocení výkonnosti modelu můžeme podívat i na hodnotu $\bar{R}^2$ (adjusted R squared), která penalizuje model za užívání "neužitečných" proměnných. V případě našeho vylepšeného modelu natrénovaného na netransformovaných datech je $\bar{R}^2$ zrhuba $0.64$, což je výrazné zlepšení oproti $\bar{R}^2$ referenčního modelu, která je zhruba $0.35$. Model natrénovaný na transformovaných datech dosáhl na hodnotu $\bar{R}^2$ zhruba $0.67$.
Zatím jsme počítali metriky modelu pouze na trénovacích datech. Abychom ověřili,
že model funguje i na jiných datech, použijeme jiný postup. Vzhledem k tomu, že
dataset o zemích světa obsahuje data pouze pro jeden okamžik v čase, nelze dobře
data rozdělit na trénovací a testovací. Proto použijeme křížovou validaci (cross
validation). Abychom tuto validaci nemuseli implmentovat ručně, tak použijeme
knihovnu sklearn
. Protože tato knihovna není přímo kompatibilní s naším
modelem, tak budeme muset vytvořit adapter, do kterého náš model zabalíme.
Následně provedeme opakovaně k-fold křížovou validaci a z výsledků spočítáme
průměr.
import statsmodels.api as sm
from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.model_selection import cross_val_score, RepeatedKFold
from numpy import mean
class SMWrapper(BaseEstimator, RegressorMixin):
"""A universal sklearn-style wrapper for statsmodels regressors"""
def __init__(self, model_class, fit_intercept=True):
self.model_class = model_class
self.fit_intercept = fit_intercept
def fit(self, X, y):
if self.fit_intercept:
X = sm.add_constant(X)
self.model_ = self.model_class(y, X)
self.results_ = self.model_.fit()
return self
def predict(self, X):
if self.fit_intercept:
X = sm.add_constant(X)
return self.results_.predict(X)
random_seed = 52
ref_model_mean_r2 = mean(
cross_val_score(
SMWrapper(sm.OLS),
ref_predicted_var,
ref_predictors,
scoring="r2",
cv=RepeatedKFold(random_state=random_seed),
)
)
ref_trans_model_mean_r2 = mean(
cross_val_score(
SMWrapper(sm.OLS),
ref_trans_predicted_var,
ref_trans_predictors,
scoring="r2",
cv=RepeatedKFold(random_state=random_seed),
)
)
print("Referenční model r^2:", round(ref_model_mean_r2, 2))
print(
"Referenční model s transformovanými daty r^2:",
round(ref_trans_model_mean_r2, 2),
)
print("Rozdíl:", round(ref_trans_model_mean_r2 - ref_model_mean_r2, 2))
print()
ref_trans_model_mean_r2 = mean(
cross_val_score(
SMWrapper(sm.OLS),
ref_trans_predicted_var,
ref_trans_predictors,
scoring="r2",
cv=RepeatedKFold(random_state=random_seed),
)
)
imp_model_mean_r2 = mean(
cross_val_score(
SMWrapper(sm.OLS),
imp_predicted_var,
imp_predictors,
scoring="r2",
cv=RepeatedKFold(random_state=random_seed),
)
)
print(
"Referenční model s tranformovanými daty r^2:",
round(ref_trans_model_mean_r2, 2),
)
print("Vylepšený model r^2:", round(imp_model_mean_r2, 2))
print("Rozdíl:", round(imp_model_mean_r2 - ref_trans_model_mean_r2, 2))
print()
imp_trans_model_mean_r2 = mean(
cross_val_score(
SMWrapper(sm.OLS),
imp_trans_predicted_var,
imp_trans_predictors,
scoring="r2",
cv=RepeatedKFold(random_state=random_seed),
)
)
print("Vylepšený model r^2:", round(imp_model_mean_r2, 2))
print(
"Vylepšený model s transformovanými daty r^2:",
round(imp_trans_model_mean_r2, 2),
)
print("Rozdíl:", round(imp_trans_model_mean_r2 - imp_model_mean_r2, 2))
print()
print(
"Celkové zlepšení:", round(imp_trans_model_mean_r2 - ref_model_mean_r2, 2)
)
Referenční model r^2: 0.31 Referenční model s transformovanými daty r^2: 0.4 Rozdíl: 0.09 Referenční model s tranformovanými daty r^2: 0.4 Vylepšený model r^2: 0.45 Rozdíl: 0.05 Vylepšený model r^2: 0.45 Vylepšený model s transformovanými daty r^2: 0.49 Rozdíl: 0.03 Celkové zlepšení: 0.18
Základní vylepšený model dává lepší hodnoty, ale vidíme, že rozdíl není tak výrazný jako v případě, kdy jsme $R^2$ počítali na trénovacích datech. Vylepšený model natrénovaný na transformovaných datech prokazuje malé zlepšení oproti vylepšenému modelu natrénovanému na netransformovaných datech.
Závěr¶
Zjistili jsme, že na základě gramotnosti populace lze alespoň částečně předpovědět skóre štěstí populace. Přestože gramotnost populace má jistou vypovídací hodnotu sama o sobě, pro přesnější predikce je vhodné data transformovat na logaritmickou škálu. Pro ještě lepší výsledky je vhodné zahrnout více faktorů, jako jsou ekonomické ukazatele, zdravotní péče apod. Nejvyšší hodnoty $R^2$ se podařilo dosáhnout kombinací lineární regrese s více proměnnými a aplikací transformace na logaritmickou škálu na trénovací data. Vylepšený model natrénovaný na transformovaných datech dosáhl $\overline{R}^2 \doteq 0.67$. Při křížové validaci dosáhl $R^2 \doteq 0.45$.
Přílohy¶
Popis sloupců datasetu World Happiness Report¶
Popis sloupců datasetu World Happiness Report, který byl převzat z přílohy 1 reportu, který byl k datasetu vydán:
- Country or region - země či region, jichž se řádek týká
- Score - Měření subjektivní pohody (SWB) pochází z vydání Gallup World Poll (GWP) z 26. února 2021 a pokrývá roky 2005 až 2020. Pokud není uvedeno jinak, jde o národní průměr odpovědí na otázku týkající se hodnocení života. Překlad znění otázky do Českého jazyka je následující: „Představte si žebřík, jehož stupně jsou očíslovány od 0 na spodním konci po 10 na horním konci. Vrchol žebříku představuje nejlepší možný život pro vás a spodní část žebříku nejhorší možný život pro vás. Na kterém stupni žebříku byste řekli, že se osobně nacházíte v tuto chvíli?“ Toto měření je také označováno jako Cantrilův žebřík života.
- Overall rank - pozice při setřízení dle štěstí populace
- GDP per capita - Hrubý domací produkt na obyvatele v paritě kupní síly (PPP) v konstantních mezinárodních cenách dolaru roku 2017 pochází z aktualizace Světových vývojových indikátorů (WDI) z 14. října 2020. Údaje o HDP pro Tchaj-wan, Sýrii, Palestinu, Venezuelu, Džibuti a Jemen pocházejí z Penn World Table 9.1.
- Social support - (nebo mít někoho, na koho se lze spolehnout v těžkých časech) je národní průměr binárních odpovědí (buď 0, nebo 1) na otázku GWP „Kdybyste měli problém, máte příbuzné nebo přátele, na které se můžete kdykoliv obrátit s žádostí o pomoc?“
- Healthy life expectancy - Očekávaná délka zdravého života při narození je založena na datech získaných z datového úložiště Světové zdravotnické organizace (WHO) Global Health Observatory (Poslední aktualizace: 28. září 2020). Data jsou dostupná pro roky 2000, 2005, 2010, 2015 a 2016. Aby odpovídala období sledovanému v této zprávě (2005–2020), byla použita interpolace a extrapolace.
- Freedom to make life choices - je národní průměr odpovědí na otázku GWP „Jste spokojeni nebo nespokojeni se svou svobodou rozhodovat, co dělat se svým životem?“
- Generosity - je reziduál regrese národního průměru odpovědí na otázku GWP „Darovali jste v minulém měsíci peníze na charitu?“ vzhledem k HDP na obyvatele.
- Perceptions of corruption - Toto měření je národní průměr odpovědí v průzkumu na dvě otázky v GWP: „Je korupce rozšířená ve vládě, nebo ne?“ a „Je korupce rozšířená v podnicích, nebo ne?“ Celkové vnímání je průměrem těchto dvou odpovědí (0 nebo 1). V případě, že chybí údaje o korupci ve vládě, používáme údaje o korupci v podnicích jako celkové vnímání. Vnímání korupce na národní úrovni je jednoduše průměrem individuálního vnímání.
Popis sloupců datasetu Countries of the World¶
Dataset je složen z několika veřejně dostupných datasetů. Popisy jednotlivých sloupců bylo potřeba najít v dokumentací zdrojových datasetů. Mnoho z těchto datasetů bohužel popisy sloupců neobsahovaly a proto je popis nahrazen pouhým překladem názvu sloupce. Nedostupnost popisu jednotlivých sloupců značně komplikuje analýzu a snížuje důvěryhodnost datasetu.
- Country – Země
- Birthrate – Porodnost
- Infant mortality (per 1000 births) – Úmrtnost kojenců (na 1000 narozených)
- GDP ($ per capita) – HDP (v dolarech na obyvatele)
- Literacy (%) – Gramotnost (%)
- Phones (per 1000) – Telefony (na 1000 obyvatel)
- Agriculture – Zemědělství
- Region – Region
- Area (sq. mi.) – Rozloha (čtvereční míle)
- Pop. Density (per sq. mi.) – Hustota obyvatelstva (na čtvereční míli)
- Coastline (coast/area ratio) – Pobřeží (poměr pobřeží k rozloze)
- Net migration – Čistá migrace
- Arable (%) – Orná půda (%)
- Crops (%) – Zemědělské plodiny (%)
- Other (%) – Ostatní (%)
- Climate – Hodnota popisuje podnebí. 1 znamená suché tropické nebo tundra a led, klasifikace B a E. 2 znamená vlhké tropické, klasifikace A. 3 znamená mírně vlhké, subtropické a mirné kontinentální, klasifikace Cfa, Cwa a D. 4 znamená suché a teplé léta a vlkhé zimy.
- Deathrate – Úmrtnost
- Industry – Průmysl
- Service – Služby
- Population – Populace