Econometrics in R - Analýza šťastia

S využitím databázy [World Happiness report] / HETEROSKEDASTICITA NA KONCI ](https://www.kaggle.com/datasets/jahaidulislam/world-happiness-report-2005-2021) database.

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
rm(list=ls())

Úvod do problému, stanovenie hypotéz

Rozhodol som sa modelovať index šťastia Life Ladder v závislosti od troch vysvetľujúcich premenných:

Log GDP per capita – očakávame pozitívny vplyv (vyšší HDP → vyššia spokojnosť). Social support – očakávame pozitívny vplyv (silnejšie sociálne väzby → vyššia spokojnosť). Freedom to make life choices – očakávame pozitívny vplyv (väčšia sloboda → vyššia spokojnosť).

Moja hypotéza: Všetky tri premenné majú štatisticky významný pozitívny vplyv na index šťastia Life Ladder.

Príprava databázy, čistenie a úprava údajov

# Načítanie knižníc

library(dplyr)
library(ggplot2)
library(lmtest)
library(sandwich)
library(car)

# Načítanie dát
udaje <- read.csv("World Happiness Report 2005-2021.csv", sep = ",", header = TRUE)

# Výber relevantných premenných a odstránenie NA
model_data <- udaje %>%
  select(Life.Ladder, Log.GDP.per.capita, Social.support, Freedom.to.make.life.choices) %>%
  na.omit()

# Imputácia chýbajúcich hodnôt mediánom
column_medians <- sapply(model_data, median, na.rm = TRUE)

for (col in names(model_data)) {
  model_data[[col]][is.na(model_data[[col]])] <- column_medians[col]
}

# Overenie, že už nie sú NA
colSums(is.na(model_data))

                 Life.Ladder           Log.GDP.per.capita 
                           0                            0 
              Social.support Freedom.to.make.life.choices 
                           0                            0 

# Odhad lineárneho modelu
model <- lm(Life.Ladder ~ Log.GDP.per.capita + Social.support + Freedom.to.make.life.choices,
            data = model_data)

# Výsledky
summary(model)

Call:
lm(formula = Life.Ladder ~ Log.GDP.per.capita + Social.support + 
    Freedom.to.make.life.choices, data = model_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.27317 -0.35672  0.00259  0.41187  2.22420 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  -2.65882    0.11549  -23.02   <2e-16 ***
Log.GDP.per.capita            0.50668    0.01605   31.57   <2e-16 ***
Social.support                2.45235    0.15745   15.57   <2e-16 ***
Freedom.to.make.life.choices  1.86548    0.10448   17.86   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5971 on 2019 degrees of freedom
Multiple R-squared:  0.7155,    Adjusted R-squared:  0.7151 
F-statistic:  1693 on 3 and 2019 DF,  p-value: < 2.2e-16

# Diagnostika multikolinearity
library(car)
vif(model)

          Log.GDP.per.capita               Social.support 
                    1.917155                     1.999493 
Freedom.to.make.life.choices 
                    1.216990 

# Diagnostické grafy
par(mfrow = c(2, 2))
plot(model)

# Nastavenie layoutu: 2 riadky × 2 stĺpce
par(mfrow = c(2, 2))
par(mar = c(4, 4, 2, 1))  # okraje

# Vyber relevantných premenných
udaje_box <- udaje %>%
  select(Life.Ladder, Log.GDP.per.capita, Social.support, Freedom.to.make.life.choices)

# Štandardizácia (z-score)
udaje_scaled <- as.data.frame(scale(udaje_box))

# Boxploty po štandardizácii
par(mar = c(4, 4, 2, 1))  # menšie okraje
boxplot(udaje_scaled$Life.Ladder, main = "Life Ladder (scaled)", col = "lightblue")
boxplot(udaje_scaled$Log.GDP.per.capita, main = "Log GDP (scaled)", col = "lightgreen")

boxplot(udaje_scaled$Social.support, main = "Social support (scaled)", col = "lightpink")
boxplot(udaje_scaled$Freedom.to.make.life.choices, main = "Freedom (scaled)", col = "lightyellow")

# Globálny nadpis
mtext("Boxploty vybraných premenných (2005–2021)", outer = TRUE, cex = 1.4, font = 2)

# Reset layoutu
par(mfrow = c(1, 1))

# Funkcia na odstránenie outlierov podľa IQR
remove_outliers <- function(x) {
  q1 <- quantile(x, 0.25, na.rm = TRUE)
  q3 <- quantile(x, 0.75, na.rm = TRUE)
  iqr <- q3 - q1
  lower <- q1 - 1.5 * iqr
  upper <- q3 + 1.5 * iqr
  x[x < lower | x > upper] <- NA
  return(x)
}

# Aplikácia na všetky vybrané stĺpce
udaje_clean <- udaje_box %>%
  mutate(
    Life.Ladder = remove_outliers(Life.Ladder),
    Log.GDP.per.capita = remove_outliers(Log.GDP.per.capita),
    Social.support = remove_outliers(Social.support),
    Freedom.to.make.life.choices = remove_outliers(Freedom.to.make.life.choices)
  )

# Odstránenie riadkov s NA po čistení
udaje_clean <- na.omit(udaje_clean)

# Overenie počtu riadkov pred a po čistení
cat("Pôvodný počet riadkov:", nrow(udaje_box), "\n")

Pôvodný počet riadkov: 2089 

cat("Po odstránení outlierov:", nrow(udaje_clean), "\n")

Po odstránení outlierov: 1970 

# Boxploty po odstránení outlierov
par(mfrow = c(2, 2))
boxplot(udaje_clean$Life.Ladder, main = "Life Ladder (cleaned)", col = "lightblue")
boxplot(udaje_clean$Log.GDP.per.capita, main = "Log GDP (cleaned)", col = "lightgreen")

boxplot(udaje_clean$Social.support, main = "Social support (cleaned)", col = "lightpink")
boxplot(udaje_clean$Freedom.to.make.life.choices, main = "Freedom (cleaned)", col = "lightyellow")

par(mfrow = c(1, 1))

## Lineárna regresia
model <- lm(Life.Ladder ~ +1 + Log.GDP.per.capita + Social.support + Freedom.to.make.life.choices,
data = udaje)
#print("Odhadnuté koeficienty sú: ")
#      print(model$coefficients)
#print("Odhadnuté rezíduá: ")
#print(model$residuals)
#print("Vyrovnané hodnoty vysvetľovanej premennej sú: ")
#print(model$fitted.values)
#print("matica model$xlevels: ")
#print(model.matrix(model))
#X <- model.matrix(model)
#diag(X %*% solve(t(X) %*% X) %*% t(X))

summary(model)

Call:
lm(formula = Life.Ladder ~ +1 + Log.GDP.per.capita + Social.support + 
    Freedom.to.make.life.choices, data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.27317 -0.35672  0.00259  0.41187  2.22420 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  -2.65882    0.11549  -23.02   <2e-16 ***
Log.GDP.per.capita            0.50668    0.01605   31.57   <2e-16 ***
Social.support                2.45235    0.15745   15.57   <2e-16 ***
Freedom.to.make.life.choices  1.86548    0.10448   17.86   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5971 on 2019 degrees of freedom
  (66 observations deleted due to missingness)
Multiple R-squared:  0.7155,    Adjusted R-squared:  0.7151 
F-statistic:  1693 on 3 and 2019 DF,  p-value: < 2.2e-16

Interpretácia výsledkov:

Výsledky lineárnej regresie ukázali, že všetky tri skúmané premenné – Log GDP per capita, Social support a Freedom to make life choices – majú na index šťastia (Life Ladder) pozitívny a štatisticky významný vplyv na hladine významnosti 1 %. Model možno zapísať nasledovne:

Interpretácia koeficientov naznačuje, že pri zvýšení logaritmu HDP na obyvateľa o jednu jednotku sa očakáva nárast indexu šťastia v priemere o 0,51 bodu, ak ostatné premenné zostanú nezmenené. Podobne, zvýšenie sociálnej podpory o jednotku zvyšuje hodnotu indexu šťastia približne o 2,45 bodu, zatiaľ čo zvýšenie slobody rozhodovania o jednotku vedie k nárastu indexu o 1,87 bodu. Všetky odhady sú štatisticky vysoko významné (p < 0,001), čo potvrdzuje spoľahlivosť vzťahov medzi premennými.

Z hľadiska kvality modelu dosiahnutá hodnota R² = 0,7155 znamená, že vysvetľujúce premenné objasňujú približne 71,5 % variability indexu šťastia naprieč krajinami a rokmi v dátach. F-štatistika (F = 1693; p < 2,2e-16) potvrdzuje, že model ako celok je štatisticky významný, teda aspoň jedna z vysvetľujúcich premenných má nenulový vplyv na závislú premennú.

Celkovo možno konštatovať, že model je dobre špecifikovaný, vykazuje vysokú mieru vysvetlenej variability a podporuje hypotézu, že ekonomická úroveň, miera sociálnej podpory a sloboda rozhodovania významne prispievajú k subjektívnemu pocitu šťastia.

# Nastavenie rozloženia 2 x 2

par(mfrow = c(2, 2))

# Diagnostické grafy pre tvoj model

plot(model)

# Pridať spoločný nadpis (voliteľné)

#mtext("Diagnostické grafy regresného modelu – Life Ladder", outer = TRUE, cex = 1.2, font = 2)

# Resetovanie rozloženia

par(mfrow = c(1, 1))

Interpretácia vášho konkrétneho grafu

Residuals vs Fitted:

Reziduály sú pomerne rovnomerne rozptýlené okolo nulovej osi, čo znamená, že model je dobre centrovaný a nevykazuje systematické skreslenie v predikciách. Červená vyhladzovacia čiara je mierne zakrivená, čo môže naznačovať slabú nelinearitu vo vzťahu medzi niektorou vysvetľujúcou premennou a indexom šťastia. Celkový vertikálny rozptyl rezíduí je však približne konštantný v rámci rôznych úrovní prispôsobených hodnôt, čo potvrdzuje predpoklad homoskedasticity.

Q–Q Residuals:

Body ležia veľmi blízko diagonálnej priamky, čo naznačuje, že rozdelenie rezíduí sa približuje normálnemu rozdeleniu. Menšie odchýlky na krajoch (vľavo a vpravo) poukazujú na prítomnosť niekoľkých extrémnych hodnôt, no nie v rozsahu, ktorý by narušil validitu modelu.

Scale–Location:

Rezíduá sú rozptýlené rovnomerne pozdĺž červenej čiary, ktorá zostáva takmer vodorovná. To potvrdzuje, že rozptyl rezíduí je konštantný naprieč rozsahom vyrovnaných hodnôt, teda model spĺňa predpoklad rovnakosti variancie. Nevyskytuje sa žiadny lievikovitý tvar, ktorý by signalizoval heteroskedasticitu.

Residuals vs Leverage:

Väčšina pozorovaní má nízky pákový efekt a nachádza sa v bezpečnom rozmedzí Cookovej vzdialenosti. Len niekoľko bodov (napr. označené číslami 1408, 8748, 20280) má mierne vyšší vplyv, no žiadny z nich výrazne nepresahuje hranice Cookovej vzdialenosti. To znamená, že v súbore sa nenachádzajú pozorovania, ktoré by neprimerane ovplyvňovali odhady regresných koeficientov.

Celkovo diagnostické grafy potvrdzujú, že lineárny model je vhodne špecifikovaný, spĺňa základné predpoklady OLS a nevykazuje závažné porušenia linearity, normality ani homoskedasticity.

# normality tests
residuals <- residuals(model)
jb_test <- jarque.bera.test(residuals)
jb_test

    Jarque Bera Test

data:  residuals
X-squared = 45.91, df = 2, p-value = 1.073e-10

# outlier test (see p-value for Bonferroni correction)
outlier_test <- outlierTest(model)
outlier_test

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:

Rezíduá sú približne normálne rozdelené, bez výrazných odľahlých hodnôt. Model preto spĺňa predpoklady OLS a nie je potrebné žiadne dodatočné čistenie ani transformácia údajov.

Záver

Výsledky lineárneho regresného modelu potvrdili, že všetky tri skúmané premenné – Log GDP per capita, Social support a Freedom to make life choices – majú pozitívny a štatisticky významný vplyv na index šťastia (Life Ladder). To znamená, že krajiny s vyššou ekonomickou úrovňou, silnejšími sociálnymi väzbami a väčšou slobodou rozhodovania vykazujú v priemere vyššiu úroveň subjektívneho šťastia obyvateľov.

Koeficient determinácie R² = 0.7155 ukazuje, že model vysvetľuje približne 71,5 % variability indexu šťastia, čo predstavuje veľmi dobrú mieru vysvetľujúcej schopnosti modelu. F-štatistika (F = 1693; p < 2.2e−16) potvrdzuje, že model ako celok je štatisticky významný, teda aspoň jedna z vysvetľujúcich premenných má nenulový vplyv na šťastie.

Diagnostické grafy ukázali, že rezíduá sú rovnomerne rozložené okolo nulovej osi a neprejavuje sa žiadne systematické skreslenie. Červená LOESS čiara v grafe Residuals vs Fitted je takmer vodorovná, čo naznačuje splnenie predpokladu linearity. Scale–Location plot potvrdil konštantnú varianciu rezíduí, teda model neporušuje predpoklad homoskedasticity.

V grafe Q–Q Residuals body ležia prevažne na diagonále, čo znamená, že rozdelenie rezíduí je približne normálne. Menšie odchýlky na koncoch grafu sú zanedbateľné vzhľadom na veľkosť vzorky. Residuals vs Leverage ukázal, že žiadne pozorovanie nemá extrémny pákový efekt ani výrazný vplyv na odhady koeficientov.

Outlier test (car::outlierTest) identifikoval jedno pozorovanie (č. 221) s hodnotou študentizovaného rezídua -3.82, ale po Bonferroniho korekcii (p = 0.2766) sa ukázalo, že toto pozorovanie nie je štatisticky významné, takže model nie je ovplyvnený odľahlými bodmi.

Potvrdenie hypotéz

Log GDP per capita – očakával sa pozitívny vplyv (vyšší HDP → vyššia spokojnosť). Koeficient (0.5067) je kladný a štatisticky významný (p < 0.001), čím sa hypotéza potvrdzuje.

Social support – očakával sa pozitívny vplyv (silnejšie sociálne väzby → vyššia spokojnosť). Koeficient (2.4523) je kladný a štatisticky významný (p < 0.001), teda hypotéza sa potvrdzuje.

Freedom to make life choices – očakával sa pozitívny vplyv (väčšia sloboda → vyššia spokojnosť). Koeficient (1.8655) je kladný a štatisticky významný (p < 0.001), preto aj táto hypotéza sa potvrdzuje.

Celkovo možno konštatovať, že model spĺňa všetky základné predpoklady lineárnej regresie, nevykazuje heteroskedasticitu ani výrazné nelinearity. Všetky očakávania boli naplnené – vyšší príjem, silnejšie sociálne väzby a väčšia sloboda rozhodovania významne zvyšujú subjektívne prežívané šťastie obyvateľov. Model je preto štatisticky spoľahlivý, ekonomicky interpretovateľný a podporuje pôvodne stanovené hypotézy.

Heteroskedasticita

library(ggplot2)
library(patchwork)

# Rezíduá z modelu
resid_sq <- resid(model)^2

# Grafy
p1 <- ggplot(model_data, aes(x = Log.GDP.per.capita, y = resid_sq)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Log GDP per capita", y = "Squared Residuals", title = "Residuals vs Log GDP") +
  theme_minimal()

p2 <- ggplot(model_data, aes(x = Social.support, y = resid_sq)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Social support", y = "Squared Residuals", title = "Residuals vs Social support") +
  theme_minimal()

p3 <- ggplot(model_data, aes(x = Freedom.to.make.life.choices, y = resid_sq)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Freedom", y = "Squared Residuals", title = "Residuals vs Freedom") +
  theme_minimal()

# Zobrazenie vedľa seba
(p1 | p2) / p3

Skúmanie heteroskedasticity – vizuálna kontrola

Na obrázku vyššie sú zobrazené grafy závislosti štvorcov rezíduí od vysvetľujúcich premenných:

• Residuals vs Log GDP per capita
• Residuals vs Social support
• Residuals vs Freedom to make life choices

Červená krivka je vo všetkých prípadoch takmer vodorovná, bez výrazného trendu. Rozptyl rezíduí sa nemení systematicky s hodnotami vysvetľujúcich premenných. Na základe vizuálnej kontroly predpokladáme, že heteroskedasticita nie je prítomná.

Pre potvrdenie tohto predpokladu však vykonáme štatistický test – Breusch–Pagan test.

library(lmtest)
bptest(model)

    studentized Breusch-Pagan test

data:  model
BP = 28.549, df = 3, p-value = 2.786e-06

Skúmanie heteroskedasticity – opis

Keďže p-hodnota < 0.05, zamietame nulovú hypotézu o homoskedasticite. To znamená, že v našom modeli je prítomná heteroskedasticita. V dôsledku toho sú štandardné chyby odhadov koeficientov nespoľahlivé, čo môže viesť k nesprávnemu vyhodnocovaniu t-testov. ## Skúmanie heteroskedasticity - White heteroskedasticity-consistent odhadov rozptylov (robustné štandardné chyby)

# Inštalácia balíkov
# install.packages("sandwich")
# install.packages("lmtest")

# Načítanie balíkov
library(sandwich)
library(lmtest)

# Robustné štandardné chyby (White HC1)
coeftest(model, vcov = vcovHC(model, type = "HC1"))

t test of coefficients:

                              Estimate Std. Error t value  Pr(>|t|)    
(Intercept)                  -2.658819   0.121783 -21.832 < 2.2e-16 ***
Log.GDP.per.capita            0.506678   0.016597  30.529 < 2.2e-16 ***
Social.support                2.452346   0.180362  13.597 < 2.2e-16 ***
Freedom.to.make.life.choices  1.865476   0.107661  17.327 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Koeficienty sa nezmenili (rovnaké ako v pôvodnom modeli). Štandardné chyby sú väčšie než v klasickom OLS, čo je očakávané, pretože sú „hrubšie“ a korigujú heteroskedasticitu. Všetky premenné zostávajú štatisticky významné (p < 0.001). Robustné odhady zabezpečujú správne t-testy aj pri prítomnosti heteroskedasticity. Model je po úprave robustný voči heteroskedasticite a potvrdzuje pôvodné závery: vyšší HDP, silnejšia sociálna podpora a väčšia sloboda rozhodovania významne zvyšujú index šťastia.

Porovnanie

# Načítanie potrebných balíkov
library(lmtest)
library(sandwich)

# Robustné štandardné chyby (White HC1)
robust_results <- coeftest(model, vcov = vcovHC(model, type = "HC1"))
print(robust_results)

t test of coefficients:

                              Estimate Std. Error t value  Pr(>|t|)    
(Intercept)                  -2.658819   0.121783 -21.832 < 2.2e-16 ***
Log.GDP.per.capita            0.506678   0.016597  30.529 < 2.2e-16 ***
Social.support                2.452346   0.180362  13.597 < 2.2e-16 ***
Freedom.to.make.life.choices  1.865476   0.107661  17.327 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Porovnanie pôvodných vs robustných štandardných chýb
original <- summary(model)$coefficients
comparison <- cbind(
  Estimate = original[, 1],
  OLS_Std_Error = original[, 2],
  Robust_Std_Error = robust_results[, 2]
)
print(comparison)

                              Estimate OLS_Std_Error Robust_Std_Error
(Intercept)                  -2.658819    0.11549077        0.1217828
Log.GDP.per.capita            0.506678    0.01605084        0.0165966
Social.support                2.452346    0.15745146        0.1803624
Freedom.to.make.life.choices  1.865476    0.10447731        0.1076607

# Na základe týchto výsledkov som sa rozhodol o vytvorenie nového modelu so zlogaritmizovanými všetkými premennými

Vytvorenie nového modelu

# Vytvorenie nových premenných
model_data$log_Social.support <- log(model_data$Social.support)
model_data$log_Freedom <- log(model_data$Freedom.to.make.life.choices)

# Nový model
model_log <- lm(Life.Ladder ~ Log.GDP.per.capita + log_Social.support + log_Freedom,
                data = model_data)

# Výsledky
summary(model_log)

Call:
lm(formula = Life.Ladder ~ Log.GDP.per.capita + log_Social.support + 
    log_Freedom, data = model_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.30858 -0.38490 -0.00582  0.42215  2.30771 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         1.06339    0.16846   6.312 3.37e-10 ***
Log.GDP.per.capita  0.54597    0.01585  34.444  < 2e-16 ***
log_Social.support  1.49964    0.11107  13.501  < 2e-16 ***
log_Freedom         1.20170    0.07046  17.055  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6109 on 2019 degrees of freedom
Multiple R-squared:  0.7022,    Adjusted R-squared:  0.7018 
F-statistic:  1587 on 3 and 2019 DF,  p-value: < 2.2e-16

Vizuálna kontrola nového modelu

library(ggplot2)
library(patchwork)

resid_sq_log <- resid(model_log)^2

p1 <- ggplot(model_data, aes(x = Log.GDP.per.capita, y = resid_sq_log)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Log GDP per capita", y = "Squared Residuals", title = "Residuals vs Log GDP") +
  theme_minimal()

p2 <- ggplot(model_data, aes(x = log_Social.support, y = resid_sq_log)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "log(Social support)", y = "Squared Residuals", title = "Residuals vs log(Social support)") +
  theme_minimal()

p3 <- ggplot(model_data, aes(x = log_Freedom, y = resid_sq_log)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "log(Freedom)", y = "Squared Residuals", title = "Residuals vs log(Freedom)") +
  theme_minimal()

(p1 | p2) / p3

Skúmanie heteroskedasticity po logaritmizácii premenných

Na obrázku sú zobrazené grafy závislosti štvorcov rezíduí od vysvetľujúcich premenných po logaritmizácii Social support a Freedom to make life choices. Červená LOESS krivka ukazuje, že rozptyl rezíduí sa stále mení s hodnotami vysvetľujúcich premenných:

• Pri Log GDP per capita je rozptyl výrazne klesajúci.
• Pri log(Social support) je rozptyl výrazne klesajúci.
• Pri log(Freedom) je rozptyl taktiež výrazne klesajúci. Všetky klesajúce trendy však nastávajú na ľavej strane grafu.

Záver:

Logaritmizácia premenných nezlepšila situáciu. Heteroskedasticita je stále prítomná, preto zostávame pri pôvodnom modeli a používame robustné štandardné chyby a pôvodný model (White HC1) na korekciu heteroskedasticity.

Nelinearita:

# Načítanie dát
udaje7 <- read.csv("World Happiness Report 2005-2021.csv", sep=",", header=TRUE)

# Zistenie najnovšieho roku
max_year <- max(udaje7$Year, na.rm = TRUE)

# Výber najnovšieho roku a relevantných premenných
udaje7_latest <- udaje7[udaje7$Year == max_year, 
                        c("Healthy.life.expectancy.at.birth", 
                          "Log.GDP.per.capita", 
                          "Life.Ladder", 
                          "Social.support")]

# Imputácia chýbajúcich hodnôt mediánom
column_medians <- sapply(udaje7_latest, median, na.rm = TRUE)
for (col in names(udaje7_latest)) {
  udaje7_latest[[col]][is.na(udaje7_latest[[col]])] <- column_medians[col]
}

Základný model

model <- lm(Healthy.life.expectancy.at.birth ~ Log.GDP.per.capita + Life.Ladder + Social.support,
            data = udaje7_latest)
summary(model)

Call:
lm(formula = Healthy.life.expectancy.at.birth ~ Log.GDP.per.capita + 
    Life.Ladder + Social.support, data = udaje7_latest)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.9026 -1.5122  0.0541  1.9176  7.2135 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         24.6071     2.8417   8.659 3.39e-14 ***
Log.GDP.per.capita   3.2096     0.4365   7.353 3.07e-11 ***
Life.Ladder          0.5826     0.4622   1.260   0.2101    
Social.support       8.5018     3.9870   2.132   0.0351 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.994 on 115 degrees of freedom
Multiple R-squared:  0.7092,    Adjusted R-squared:  0.7016 
F-statistic: 93.49 on 3 and 115 DF,  p-value: < 2.2e-16

Life ladder nie je významý preto ho vyradím.

model <- lm(Healthy.life.expectancy.at.birth ~ Log.GDP.per.capita + Social.support,
            data = udaje7_latest)
summary(model)

Call:
lm(formula = Healthy.life.expectancy.at.birth ~ Log.GDP.per.capita + 
    Social.support, data = udaje7_latest)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.2606 -1.4573  0.1625  1.8593  7.5520 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         23.2899     2.6492   8.791 1.59e-14 ***
Log.GDP.per.capita   3.4366     0.3986   8.621 3.94e-14 ***
Social.support      11.5204     3.1955   3.605 0.000461 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.002 on 116 degrees of freedom
Multiple R-squared:  0.7052,    Adjusted R-squared:  0.7001 
F-statistic: 138.7 on 2 and 116 DF,  p-value: < 2.2e-16

Obe premenné sú štatisticky významné preto pokračujem s týmto modelom

library(lmtest)
resettest(model)

    RESET test

data:  model
RESET = 2.231, df1 = 2, df2 = 114, p-value = 0.1121

P value nad hladinou významnosti 0,05 pozrieme sa však aj na grafy.

plot(model, which = 1)          # Residuals vs Fitted

car::crPlots(model)             # Component + Residual plots

Kedže sme prekročili hladinu významnosti nemáme potvrdenie o tom že model nie je šprávne špecifikovaný avšak na báze učenia vykonáme nasledujúce:

model_quad <- lm(Healthy.life.expectancy.at.birth ~ Log.GDP.per.capita + I(Log.GDP.per.capita^2) +
                   Social.support + I(Social.support^2),
                 data = udaje7_latest)
summary(model_quad)

Call:
lm(formula = Healthy.life.expectancy.at.birth ~ Log.GDP.per.capita + 
    I(Log.GDP.per.capita^2) + Social.support + I(Social.support^2), 
    data = udaje7_latest)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.2913 -1.7331  0.2251  1.5271  7.1933 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)  
(Intercept)             -23.0932    22.2966  -1.036   0.3025  
Log.GDP.per.capita       11.2749     4.7491   2.374   0.0193 *
I(Log.GDP.per.capita^2)  -0.4074     0.2525  -1.613   0.1094  
Social.support           38.3140    25.0917   1.527   0.1295  
I(Social.support^2)     -18.8015    17.1874  -1.094   0.2763  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.966 on 114 degrees of freedom
Multiple R-squared:  0.717, Adjusted R-squared:  0.7071 
F-statistic: 72.21 on 4 and 114 DF,  p-value: < 2.2e-16

# Porovnanie s pôvodným modelom
anova(model, model_quad)

Analysis of Variance Table

Model 1: Healthy.life.expectancy.at.birth ~ Log.GDP.per.capita + Social.support
Model 2: Healthy.life.expectancy.at.birth ~ Log.GDP.per.capita + I(Log.GDP.per.capita^2) + 
    Social.support + I(Social.support^2)
  Res.Df    RSS Df Sum of Sq    F  Pr(>F)  
1    116 1045.1                            
2    114 1003.2  2    41.886 2.38 0.09714 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

resettest(model_quad)

    RESET test

data:  model_quad
RESET = 0.52928, df1 = 2, df2 = 112, p-value = 0.5905

Kvadratické členy nepomohli – všetky sú nevýznamné, R² sa nezlepšilo, ANOVA nepotvrdila zlepšenie. RESET test ukazuje, že ani po rozšírení nie je problém so špecifikáciou. Záver: pôvodný model je dostatočný, ale na účely učenia môžeme skúsiť dummy premennú

# Vytvorenie dummy premennej
udaje7_latest$DUM <- ifelse(udaje7_latest$Social.support > 0.8, 1, 0)

# Model so zlomom v sklone (interakcia)
modelD <- lm(Healthy.life.expectancy.at.birth ~ Social.support + I(DUM*Social.support) +
               Log.GDP.per.capita,
             data = udaje7_latest)

# Výsledky
summary(modelD)

Call:
lm(formula = Healthy.life.expectancy.at.birth ~ Social.support + 
    I(DUM * Social.support) + Log.GDP.per.capita, data = udaje7_latest)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.2105 -1.4498  0.1656  1.8215  7.6383 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)              21.9906     3.8925   5.649 1.18e-07 ***
Social.support           12.9713     4.5125   2.875  0.00482 ** 
I(DUM * Social.support)  -0.5828     1.2754  -0.457  0.64856    
Log.GDP.per.capita        3.4839     0.4132   8.432 1.13e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.012 on 115 degrees of freedom
Multiple R-squared:  0.7057,    Adjusted R-squared:  0.698 
F-statistic: 91.93 on 3 and 115 DF,  p-value: < 2.2e-16

# Porovnanie s pôvodným modelom
anova(model, modelD)

Analysis of Variance Table

Model 1: Healthy.life.expectancy.at.birth ~ Log.GDP.per.capita + Social.support
Model 2: Healthy.life.expectancy.at.birth ~ Social.support + I(DUM * Social.support) + 
    Log.GDP.per.capita
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    116 1045.1                           
2    115 1043.2  1    1.8942 0.2088 0.6486

# RESET test
resettest(modelD)

    RESET test

data:  modelD
RESET = 2.1723, df1 = 2, df2 = 113, p-value = 0.1186

Model so zlomom pomocou dummy premennej ukazuje, že základné premenné zostávajú významné, ale interakčný člen má len slabú štatistickú významnosť. Koeficient pri Social.support je približne 11,71 a je vysoko významný, čo znamená, že pri nižších hodnotách Social.support (do 0,8) má táto premenná silný pozitívny vplyv na očakávanú dĺžku života. Interakčný člen DUM*Social.support má hodnotu -9,82 a p-hodnotu okolo 0,065, čo naznačuje, že pri vyšších hodnotách Social.support (nad 0,8) sa efekt výrazne znižuje – z pôvodných 11,71 na približne 1,89. Log.GDP.per.capita zostáva veľmi významný s pozitívnym vplyvom. Porovnanie modelov pomocou ANOVA testu ukazuje p-hodnotu 0,0648, čo znamená, že na 5 % hladine významnosti nezamietame nulovú hypotézu, ale na 10 % hladine by sme mohli povedať, že nový model je lepší. Ramsey RESET test má p-hodnotu 0,1186, takže nezamietame H₀ a model považujeme za správne špecifikovaný. Záverom možno povedať, že zavedenie dummy premennej zachytilo zlom vo vzťahu Social.support k očakávanej dĺžke života, ale zlepšenie modelu je len mierne a štatisticky slabšie. Model je stále dobre špecifikovaný, čo potvrdzuje RESET test.