S využitím databázy World Population Data.

library(zoo)

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(tseries)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(lmtest)
library(sandwich)
library(car)

## Loading required package: carData

rm(list=ls())

Úvod do problému

Rozhodla som sa modelovať mieru rastu populácie (growth.rate) v závislosti od troch vysvetľujúcich premenných, a to celkového počtu obyvateľov v roku 2023 (X2023.population), rozlohy krajiny v km² (area..km..) a hustoty obyvateľstva (density..km..). Pracovná hypotéza hovorí o štatisticky významnom vplyve všetkých troch vysvetľujúcich premenných. Očakáva sa, že:

väčšia populácia môže viesť k nižšej miere rastu (negatívny vplyv – efekt nasýtenia),
väčšia rozloha by mohla byť spojená s mierne vyšším rastom (pozitívny vplyv – viac priestoru na rozvoj),
vyššia hustota obyvateľstva môže viesť k nižšiemu rastu (negatívny vplyv – preplnenie, urbanizácia).

Cieľom modelu je teda zistiť, ktoré z týchto faktorov štatisticky významne ovplyvňujú rast populácie krajín sveta.

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

udaje <- read.csv("world_population_data.csv", dec = ".", sep = ",", header = TRUE)
names(udaje)

##  [1] "rank"             "cca3"             "country"          "continent"       
##  [5] "X2023.population" "X2022.population" "X2020.population" "X2015.population"
##  [9] "X2010.population" "X2000.population" "X1990.population" "X1980.population"
## [13] "X1970.population" "area..km.."       "density..km.."    "growth.rate"     
## [17] "world.percentage"

udaje_2023 <- udaje[, c("X2023.population", "area..km..", "density..km..", "growth.rate")]

udaje_2023$`growth.rate` <- as.numeric(gsub("%", "", as.character(udaje_2023$`growth.rate`)))
udaje_2023$`growth.rate` <- udaje_2023$`growth.rate` / 100
head(udaje_2023$`growth.rate`)

## [1]  0.0081 -0.0002  0.0050  0.0074  0.0198  0.0241

model <- lm(growth.rate ~ X2023.population + area..km.. + density..km.., data = udaje_2023)

summary(model)

## 
## Call:
## lm(formula = growth.rate ~ X2023.population + area..km.. + density..km.., 
##     data = udaje_2023)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.084415 -0.007330 -0.001634  0.007371  0.039881 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.002e-02  8.819e-04  11.360   <2e-16 ***
## X2023.population -1.453e-12  6.628e-12  -0.219    0.827    
## area..km..       -3.446e-11  5.177e-10  -0.067    0.947    
## density..km..    -4.690e-07  4.111e-07  -1.141    0.255    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01239 on 230 degrees of freedom
## Multiple R-squared:  0.005895,   Adjusted R-squared:  -0.007071 
## F-statistic: 0.4546 on 3 and 230 DF,  p-value: 0.7143

Boxploty

library(ggplot2)

# Predpríprava dát
udaje_plot <- udaje_2023[, c("X2023.population", "area..km..", "density..km..", "growth.rate")]

# Odstránenie NA a percent z growth rate
udaje_plot$growth.rate <- as.numeric(gsub("%", "", as.character(udaje_plot$growth.rate)))
udaje_plot$X2023.population <- as.numeric(as.character(udaje_plot$X2023.population))
udaje_plot$area..km.. <- as.numeric(as.character(udaje_plot$area..km..))
udaje_plot$density..km.. <- as.numeric(as.character(udaje_plot$density..km..))

# Ošetrenie NA hodnot mediánom
for(col in names(udaje_plot)){
  udaje_plot[[col]][is.na(udaje_plot[[col]])] <- median(udaje_plot[[col]], na.rm = TRUE)
}

# Population (log scale, horizontálny)
ggplot(udaje_plot, aes(y = X2023.population)) +
  geom_boxplot(fill = "skyblue", outlier.colour = "red", outlier.shape = 16, outlier.size = 2) +
  scale_y_log10() +
  labs(title = "Population 2023",
       y = "Population (log10)",
       x = "") +
  theme_minimal(base_size = 14)

# Area (log scale)
ggplot(udaje_plot, aes(y = area..km..)) +
  geom_boxplot(fill = "lightgreen", outlier.colour = "red", outlier.shape = 16, outlier.size = 2) +
  scale_y_log10() +
  labs(title = "Area (km²)",
       y = "Area (log10 km²)",
       x = "") +
  theme_minimal(base_size = 14)

# Density (pôvodná škála)
ggplot(udaje_plot, aes(y = density..km..)) +
  geom_boxplot(fill = "orange", outlier.colour = "red", outlier.shape = 16, outlier.size = 2) +
  labs(title = "Population Density (km²)",
       y = "Density",
       x = "") +
  theme_minimal(base_size = 14)

# Growth rate (pôvodná škála, %)
ggplot(udaje_plot, aes(y = growth.rate)) +
  geom_boxplot(fill = "pink", outlier.colour = "red", outlier.shape = 16, outlier.size = 2) +
  labs(title = "Growth Rate (%)",
       y = "Growth rate (%)",
       x = "") +
  theme_minimal(base_size = 14)

Lineárna regresia

Model odhadujeme príkazom lm()

udaje_model <- udaje_2023[, c("growth.rate", "X2023.population", "area..km..", "density..km..")]
udaje_model$growth.rate <- as.numeric(gsub("%","", as.character(udaje_model$growth.rate)))
model <- lm(`growth.rate` ~ `X2023.population` + `area..km..` + `density..km..`, data = udaje_model)

summary(model)

## 
## Call:
## lm(formula = growth.rate ~ X2023.population + area..km.. + density..km.., 
##     data = udaje_model)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.084415 -0.007330 -0.001634  0.007371  0.039881 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.002e-02  8.819e-04  11.360   <2e-16 ***
## X2023.population -1.453e-12  6.628e-12  -0.219    0.827    
## area..km..       -3.446e-11  5.177e-10  -0.067    0.947    
## density..km..    -4.690e-07  4.111e-07  -1.141    0.255    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01239 on 230 degrees of freedom
## Multiple R-squared:  0.005895,   Adjusted R-squared:  -0.007071 
## F-statistic: 0.4546 on 3 and 230 DF,  p-value: 0.7143

Na základe odhadnutého lineárneho modelu možno konštatovať, že žiadna z vysvetľujúcich premenných (rozloha, hustota obyvateľstva) nemá štatisticky významný vplyv na mieru rastu populácie. Model ako celok nie je štatisticky významný (p = 0.7143) a vysvetľuje iba 0,6 % variability rastu populácie (R² = 0.0059). Intercept (≈ 0.01) naznačuje, že priemerný ročný rast populácie krajín je približne 1 %. Na zlepšenie modelu by bolo vhodné zahrnúť ďalšie premenné, napríklad demografické alebo ekonomické ukazovatele.

Inými slovami – miera rastu populácie závisí pravdepodobne od iných faktorov (napr. ekonomická úroveň, veková štruktúra, migrácia, pôrodnosť a úmrtnosť), nie od samotnej veľkosti krajiny či jej hustoty.

Diagnostické grafy regresného modelu

par(mfrow = c(2, 2))  # horný okraj pre nadpis
plot(model)

Diagnostické grafy regresného modelu

par(mfrow = c(1, 1))  # reset layout

Výsledky diagnostiky

Residuals vs. fitted

Reziduá sú veľmi blízko nuly a väčšina bodov sa zhromažďuje pri vyšších hodnotách predikcie (okolo 0.008–0.01). Nie je tam žiadny výrazný vzor alebo oblúk, takže model nemá veľké systematické chyby. Malá variabilita reziduí → dáta majú nízky rozptyl rastu populácie.

Q-Q plot

Body ležia približne na diagonále, hoci niektoré sú mierne odchýlené na krajoch. Znamená to, že reziduá sú približne normálne rozložené, ale môže byť mierna špicatosť. Pre tento dataset je to v poriadku.

Scale-Location

Body sú veľmi nízko na osi Y a tesne pri sebe. Rozptyl reziduí je veľmi malý a rovnomerný, teda homoskedasticita je zachovaná. Graf nevykazuje žiadny lievikovitý tvar.

Residuals vs Leverage

Všetky body majú veľmi nízku pákovú hodnotu (leverage), žiadny bod nevystupuje ako výrazný outlier.Cookova vzdialenosť je veľmi malá → žiadne pozorovanie nemá neprimeraný vplyv na model.

Záver

Reziduá sú blízko nuly a homogénne → model spĺňa predpoklad lineárnosti a homoskedasticity.
Q-Q graf naznačuje len miernu odchýlku od normality, nie kritickú.
Žiadne extrémne alebo vplyvné pozorovania → model nie je ťahaný outlierom.
Model je teda štatisticky bezpečný na interpretáciu, aj keď samotný R² a významnosť koeficientov ukazujú, že vysvetľuje len veľmi malú časť variability rastu populácie.

model <- lm(growth.rate ~ I(log(X2023.population)) + area..km.. + density..km.., data = udaje_2023)
summary(model)

## 
## Call:
## lm(formula = growth.rate ~ I(log(X2023.population)) + area..km.. + 
##     density..km.., data = udaje_2023)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.087225 -0.006642 -0.000998  0.007031  0.039392 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              -6.828e-03  4.555e-03  -1.499 0.135213    
## I(log(X2023.population))  1.149e-03  3.055e-04   3.761 0.000215 ***
## area..km..               -7.333e-10  4.802e-10  -1.527 0.128134    
## density..km..            -3.084e-07  4.013e-07  -0.768 0.443058    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01203 on 230 degrees of freedom
## Multiple R-squared:  0.06328,    Adjusted R-squared:  0.05107 
## F-statistic:  5.18 on 3 and 230 DF,  p-value: 0.001759

par(mfrow = c(2, 2))
plot(model)

Diagnostické grafy regresného modelu

par(mfrow = c(1, 1))

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

## 
##  Jarque Bera Test
## 
## data:  residuals
## X-squared = 1383.8, df = 2, p-value < 2.2e-16

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

##    rstudent unadjusted p-value Bonferroni p
## 41 -8.28151         1.0166e-14   2.3788e-12

V modeli sa zistilo, že rast populácie je pozitívne spojený s veľkosťou populácie – čím väčšia krajina podľa počtu obyvateľov, tým vyšší je jej rast. Naopak, rozloha krajiny a hustota obyvateľstva nemajú významný vplyv, ich účinok je v rámci náhodnej variability.

Rezíduá modelu nevykazujú presne normálne rozdelenie, ale vzhľadom na veľký počet pozorovaní je model stále spoľahlivý a vhodný na analýzu. Nezaznamenali sme žiadne závažné odchýlky od linearity, takže predpoklady regresie sú v podstate splnené.

library(ggplot2)

ggplot(udaje_2023, aes(x = log(X2023.population), y = growth.rate, 
                       size = density..km.., color = area..km..)) +
  geom_point(alpha = 0.7) +
  scale_color_gradient(low = "lightblue", high = "darkblue") +
  labs(x = "Logaritmus populácie", y = "Rast populácie",
       size = "Hustota", color = "Rozloha") +
  theme_minimal()

Interpretácia:

Graf umožňuje vizuálne sledovať, ako rast populácie súvisí s veľkosťou populácie, a zároveň zobrazuje ďalšie dve premenné naraz (rozlohu a hustotu). Väčšie bubliny (vyššia hustota) sa dajú ľahko identifikovať a farebná škála ukazuje, ktoré krajiny sú väčšie podľa rozlohy. Je vidieť, že krajiny s vyššou populáciou majú mierne vyšší rast, zatiaľ čo hustota a rozloha sa v rámci tohto grafu vizuálne prejavujú rôznorodo.

library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

plot_ly(udaje_2023, x = ~log(X2023.population), y = ~area..km.., z = ~growth.rate,
        type = "scatter3d", mode = "markers",
        marker = list(size = 5, color = ~density..km.., colorscale = 'Viridis')) %>%
  layout(scene = list(xaxis = list(title = 'Log populácie'),
                      yaxis = list(title = 'Rozloha'),
                      zaxis = list(title = 'Rast')))

Interpretácia:

Graf umožňuje vizuálne sledovať trojrozmerný vzťah medzi populáciou, rozlohou a rastom. Krajiny s vysokou populáciou a veľkou rozlohou možno identifikovať ako extrémne body v priestore grafu. Farebné kódovanie hustoty obyvateľstva pridáva ďalšiu dimenziu informácie – napríklad krajiny s vysokou hustotou možno ľahko vizuálne odlíšiť od tých s nízkou hustotou. Interaktivita (otočenie a zoom) umožňuje lepšie skúmať súvislosti a odľahlé hodnoty, ktoré by v 2D grafe mohli zostať skryté.

Heteroskedasticita

Prítomnosť heteroskedasticity, čiže nekonštantného rozptylu náhodnej zložky, môže viesť k nesprávnemu odhadu štandardných chýb regresného koeficientu a tým aj k nespoľahlivému t-testu významnosti. Preto je vhodné heteroskedasticitu najprv detegovať (vizuálne a testami) a v prípade potreby ju eliminovať.

V mojom prípade som skúmala rezíduá modelu, ktorý predikoval growth.rate na základe demografických premenných: X2023.population, area..km.. a density..km…

library(ggplot2)
library(patchwork)

# Rezíduá modelu
udaje_2023$resid2 <- resid(model)^2   

# Nahradenie nulových/NA hodnôt pre log škálu
udaje_2023$X2023.population[udaje_2023$X2023.population <= 0 | is.na(udaje_2023$X2023.population)] <- 1
udaje_2023$density..km..[udaje_2023$density..km.. <= 0 | is.na(udaje_2023$density..km..)] <- 1
udaje_2023$area..km..[udaje_2023$area..km.. <= 0 | is.na(udaje_2023$area..km..)] <- 1

# Grafy
p1 <- ggplot(udaje_2023, aes(x = X2023.population, y = resid2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  scale_x_log10(labels = scales::comma) +
  labs(x = "Population (2023)",
       y = "Squared Residuals",
       title = "Squared Residuals \nvs Population (2023)") +
  theme_minimal()

p2 <- ggplot(udaje_2023, aes(x = density..km.., y = resid2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  scale_x_log10() +
  labs(x = "Population Density \n(per km²)",
       y = "Squared Residuals",
       title = "Squared Residuals \nvs Population Density") +
  theme_minimal()

p3 <- ggplot(udaje_2023, aes(x = area..km.., y = resid2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  scale_x_log10() +
  labs(x = "Area (km²)",
       y = "Squared Residuals",
       title = "Squared Residuals \nvs Area") +
  theme_minimal()

# Zobrazenie vedľa seba
p1 + p2 + p3

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Zlogaritmované grafy

library(ggplot2)
library(patchwork)

# Vytvoríme stĺpec so štvorcami rezíduí
udaje_2023$resid2 <- resid(model)^2

# Ošetrenie nulových alebo NA hodnôt pred logaritmom
udaje_2023$X2023.population[udaje_2023$X2023.population <= 0 | is.na(udaje_2023$X2023.population)] <- 1
udaje_2023$density..km..[udaje_2023$density..km.. <= 0 | is.na(udaje_2023$density..km..)] <- 1
udaje_2023$area..km..[udaje_2023$area..km.. <= 0 | is.na(udaje_2023$area..km..)] <- 1

# Spoločný štýl pre všetky grafy
common_theme <- theme(
  plot.title = element_text(size = 12, margin = margin(b = 8)),
  plot.margin = margin(10, 10, 10, 10))

# Graf 1: Rezíduá vs log(populácia)
p1 <- ggplot(udaje_2023, aes(x = X2023.population, y = resid2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  scale_x_log10(labels = scales::comma) +
  labs(x = "Population (2023, \n(log))",
    y = "Squared Residuals",
    title = "Squared Residuals \nvs Population (2023)") + theme_minimal() + common_theme

# Graf 2: Rezíduá vs hustota obyvateľstva 
p2 <- ggplot(udaje_2023, aes(x = density..km.., y = resid2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  scale_x_log10(labels = scales::comma) +
  labs(x = "Population Density (per km², \n(log))",
    y = "Squared Residuals",
    title = "Squared Residuals \nvs Population Density") +
  theme_minimal() + common_theme

# Graf 3: Rezíduá vs rozloha krajiny 
p3 <- ggplot(udaje_2023, aes(x = area..km.., y = resid2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  scale_x_log10(labels = scales::comma) +
  labs(x = "Area (km², \n(log))",
    y = "Squared Residuals",
    title = "Squared Residuals \nvs Area") + theme_minimal() + common_theme

# Zlúčenie všetkých troch grafov vedľa seba 
(p1 + p2 + p3) +
  plot_layout(ncol = 3) +
  plot_annotation(
    title = "Squared Residuals vs Demographic Variables (2023)",
    theme = theme(plot.title = element_text(hjust = 0.5, size = 14)))

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Zlogaritmované grafy bez odľahlých hodnôt

Pre presnejšiu vizualizáciu som tiež odstránila odľahlé hodnoty rezíduí (pomocou interkvartilového rozpätia, IQR), čím sa grafy stali ešte prehľadnejšími a indikovali stabilný rozptyl.

library(ggplot2)
library(patchwork)

udaje_2023$resid2 <- resid(model)^2 

udaje_2023$X2023.population[udaje_2023$X2023.population <= 0 | is.na(udaje_2023$X2023.population)] <- 1
udaje_2023$density..km..[udaje_2023$density..km.. <= 0 | is.na(udaje_2023$density..km..)] <- 1
udaje_2023$area..km..[udaje_2023$area..km.. <= 0 | is.na(udaje_2023$area..km..)] <- 1

Q1 <- quantile(udaje_2023$resid2, 0.25)
Q3 <- quantile(udaje_2023$resid2, 0.75)
IQR_val <- Q3 - Q1
lower <- Q1 - 1.5 * IQR_val
upper <- Q3 + 1.5 * IQR_val

udaje_filtered <- udaje_2023[udaje_2023$resid2 >= lower & udaje_2023$resid2 <= upper, ]

common_theme <- theme(
  plot.title = element_text(size = 12, margin = margin(b = 8)),
  plot.margin = margin(10, 10, 10, 10))

p1 <- ggplot(udaje_filtered, aes(x = X2023.population, y = resid2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  scale_x_log10(labels = scales::comma) +
  labs(x = "Population 2023 \n(log)",
       y = "Squared Residuals",
       title = "Squared Residuals \nvs Population") +
  theme_minimal() + common_theme

p2 <- ggplot(udaje_filtered, aes(x = density..km.., y = resid2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  scale_x_log10(labels = scales::comma) +
  labs(x = "Population Density \n(log)",
       y = "Squared Residuals",
       title = "Squared Residuals \nvs Population Density") +
  theme_minimal() + common_theme

p3 <- ggplot(udaje_filtered, aes(x = area..km.., y = resid2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  scale_x_log10(labels = scales::comma) +
  labs(x = "Area (km², \n(log))",
       y = "Squared Residuals",
       title = "Squared Residuals \nvs Area") +
  theme_minimal() + common_theme

(p1 + p2 + p3) +
  plot_layout(ncol = 3) +
  plot_annotation(
    title = "Squared Residuals vs Demographic Variables (2023)") &
  theme(plot.title = element_text(hjust = 0.5, size = 14))

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

# Install (if not yet installed)
# install.packages("lmtest")

# Load the package
library(lmtest)

# Run the Breusch–Pagan test
bptest(model)

## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 4.0687, df = 3, p-value = 0.2541

summary(udaje_model[, c("X2023.population", "area..km..", "density..km..")])

##  X2023.population      area..km..        density..km..    
##  Min.   :5.180e+02   Min.   :4.000e-01   Min.   :    0.0  
##  1st Qu.:4.226e+05   1st Qu.:2.650e+03   1st Qu.:   39.5  
##  Median :5.644e+06   Median :8.120e+04   Median :   97.5  
##  Mean   :3.437e+07   Mean   :5.814e+05   Mean   :  451.3  
##  3rd Qu.:2.325e+07   3rd Qu.:4.304e+05   3rd Qu.:  242.8  
##  Max.   :1.429e+09   Max.   :1.710e+07   Max.   :21403.0

model_log <- lm(growth.rate ~ log(X2023.population) + log(area..km..) + log(density..km.. + 1), data = udaje_model)

# Install (if not yet installed)
# install.packages("lmtest")

# Load the package
library(lmtest)

# Run the Breusch–Pagan test
bptest(model_log)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_log
## BP = 3.1807, df = 3, p-value = 0.3646

#install.packages("sandwich")
#install.packages("lmtest")
library(sandwich)
library(lmtest)
coeftest(model, vcov = vcovHC(model))

## 
## t test of coefficients:
## 
##                             Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)              -6.8280e-03  3.5698e-03 -1.9127   0.05703 .  
## I(log(X2023.population))  1.1490e-03  2.6711e-04  4.3016 2.508e-05 ***
## area..km..               -7.3330e-10  3.0318e-10 -2.4187   0.01635 *  
## density..km..            -3.0835e-07  5.1667e-07 -0.5968   0.55122    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretácia skúmania heteroskedasticity

Po spustení Breusch–Pagan testu pre logaritmizovaný model som získala:

– BP test pre pôvodný model: p-hodnota = 0.2541

– BP test pre log-transformovaný model: p-hodnota = 0.3646

Pre oba modely sú p-hodnoty väčšie ako 0,05, čo znamená, že neexistuje štatisticky významný dôkaz heteroskedasticity rezíduí.

Logaritmická transformácia premenných navyše ešte mierne zlepšila stabilitu rozptylu rezíduí.

Grafická kontrola štvorcov rezíduí potvrdzuje, že rozptyl rezíduí je približne konštantný naprieč hodnotami vysvetľujúcich premenných.

Nelineárne špecifikácie

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)

rm(list = ls())

udaje <- read.csv("world_population_data.csv", dec = ".", sep = ",", header = TRUE)

udaje_2023 <- udaje[, c("X2023.population", "area..km..", "density..km..", "growth.rate")]

# growth.rate máš v %, odstránime znak % a premeníme na čislo (v PERCENTECH, NIE ako podiel)
udaje_2023$growth.rate <- as.numeric(gsub("%", "", as.character(udaje_2023$growth.rate)))

column_medians <- sapply(udaje_2023, median, na.rm = TRUE)

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

udaje_2023 <- udaje_imputed
udaje <- udaje_2023

model_lin <- lm(growth.rate ~ +1 + X2023.population + area..km.. + density..km.., data = udaje)
summary(model_lin)

## 
## Call:
## lm(formula = growth.rate ~ +1 + X2023.population + area..km.. + 
##     density..km.., data = udaje)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.4415 -0.7330 -0.1634  0.7371  3.9881 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.002e+00  8.819e-02  11.360   <2e-16 ***
## X2023.population -1.453e-10  6.628e-10  -0.219    0.827    
## area..km..       -3.446e-09  5.177e-08  -0.067    0.947    
## density..km..    -4.690e-05  4.111e-05  -1.141    0.255    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.239 on 230 degrees of freedom
## Multiple R-squared:  0.005895,   Adjusted R-squared:  -0.007071 
## F-statistic: 0.4546 on 3 and 230 DF,  p-value: 0.7143

RESET test pre základný model

library(lmtest)

# Ramsey RESET pre základný lineárny model
resettest(model_lin)                           # default (power 2:3, fitted)

## 
##  RESET test
## 
## data:  model_lin
## RESET = 1.2307, df1 = 2, df2 = 228, p-value = 0.294

resettest(model_lin, power = 2, type = "fitted")

## 
##  RESET test
## 
## data:  model_lin
## RESET = 2.4608, df1 = 1, df2 = 229, p-value = 0.1181

resettest(model_lin, power = 2:3, type = "fitted")

## 
##  RESET test
## 
## data:  model_lin
## RESET = 1.2307, df1 = 2, df2 = 228, p-value = 0.294

RESET test pre logaritmovaný model

model_log <- lm(growth.rate ~ log(X2023.population) + log(area..km..) + log(density..km.. + 1), data = udaje)

resettest(model_log)

## 
##  RESET test
## 
## data:  model_log
## RESET = 1.5932, df1 = 2, df2 = 228, p-value = 0.2055

resettest(model_log, power = 2, type = "fitted")

## 
##  RESET test
## 
## data:  model_log
## RESET = 0.50262, df1 = 1, df2 = 229, p-value = 0.4791

resettest(model_log, power = 2:3, type = "fitted")

## 
##  RESET test
## 
## data:  model_log
## RESET = 1.5932, df1 = 2, df2 = 228, p-value = 0.2055

Interpretácia RESET testu

Na overenie, či sú naše lineárne a logaritmické modely správne špecifikované, sme vykonali Ramseyho RESET test. Tento test skúma, či by pridanie nelineárnych členov predikovaných hodnôt významne zlepšilo model.

Lineárny model — RESET test: p = 0.294 a p = 0.1181

Interpretácia: p-hodnota > 0.05 → neexistuje štatisticky významný dôkaz nesprávnej špecifikácie. Model je teda vhodne špecifikovaný a netreba pridávať nelineárne členy ani ďalšie premenné.

Logaritmický model — RESET test: p = 0.2055 a p = 0.4791

Interpretácia: rovnako ako pri lineárnom modeli, neexistuje dôkaz nesprávnej špecifikácie.

Záver: Oba modely, lineárny aj logaritmický, sú podľa RESET testu štatisticky správne špecifikované. To znamená, že funkčná forma modelu je vhodná a nie sú zjavné chýbajúce premenné alebo potreba nelineárnych transformácií.

1. Residuals vs Fitted

# Lineárny model
plot(model_lin, which = 1, main = "Residuals vs Fitted (Lineárny model)")

# Logaritmický model
plot(model_log, which = 1, main = "Residuals vs Fitted (Logaritmický model)")

2. Component + Residual (C+R grafy)

library(car)

# Lineárny model
crPlots(model_lin, main = "C+R plots (Lineárny model)")

# Logaritmický model
crPlots(model_log, main = "C+R plots (Logaritmický model)")

3. Nelineárna špecifikácia - porovnanie základného a modifikovaného modelu

# Pôvodný model (lineárny)
model_lin <- lm(growth.rate ~ X2023.population + area..km.. + density..km.., data = udaje_2023)
summary(model_lin)

## 
## Call:
## lm(formula = growth.rate ~ X2023.population + area..km.. + density..km.., 
##     data = udaje_2023)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.4415 -0.7330 -0.1634  0.7371  3.9881 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.002e+00  8.819e-02  11.360   <2e-16 ***
## X2023.population -1.453e-10  6.628e-10  -0.219    0.827    
## area..km..       -3.446e-09  5.177e-08  -0.067    0.947    
## density..km..    -4.690e-05  4.111e-05  -1.141    0.255    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.239 on 230 degrees of freedom
## Multiple R-squared:  0.005895,   Adjusted R-squared:  -0.007071 
## F-statistic: 0.4546 on 3 and 230 DF,  p-value: 0.7143

# Modifikovaný model s kvadratickými členmi pre X2023.population a density
model_quad <- lm(growth.rate ~ X2023.population + area..km.. + density..km.. + 
                   I(X2023.population^2) + I(density..km..^2), data = udaje_2023)
summary(model_quad)

## 
## Call:
## lm(formula = growth.rate ~ X2023.population + area..km.. + density..km.. + 
##     I(X2023.population^2) + I(density..km..^2), data = udaje_2023)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.5325 -0.7013 -0.1495  0.7482  3.9937 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.005e+00  9.956e-02  10.091   <2e-16 ***
## X2023.population       3.428e-09  2.388e-09   1.435   0.1525    
## area..km..            -4.533e-08  5.597e-08  -0.810   0.4188    
## density..km..         -2.762e-04  1.436e-04  -1.924   0.0556 .  
## I(X2023.population^2) -2.534e-18  1.653e-18  -1.533   0.1266    
## I(density..km..^2)     1.285e-08  7.628e-09   1.685   0.0933 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.231 on 228 degrees of freedom
## Multiple R-squared:  0.02799,    Adjusted R-squared:  0.00667 
## F-statistic: 1.313 on 5 and 228 DF,  p-value: 0.2593

# Porovnanie modelov pomocou ANOVA
anova(model_lin, model_quad)

# RESET test pre nelineárny model
library(lmtest)
resettest(model_quad)

## 
##  RESET test
## 
## data:  model_quad
## RESET = 1.133, df1 = 2, df2 = 226, p-value = 0.3239

4. Rozšírený RESET test pre nelineárny model

# Vytvorenie rozšíreného modelu s kvadratickými členmi
model_rozsireny <- lm(growth.rate ~ X2023.population + area..km.. + density..km.. + I(X2023.population^2) + I(area..km..^2) + I(density..km..^2), data = udaje_2023)

# Zobrazenie výsledkov modelu
summary(model_rozsireny)

## 
## Call:
## lm(formula = growth.rate ~ X2023.population + area..km.. + density..km.. + 
##     I(X2023.population^2) + I(area..km..^2) + I(density..km..^2), 
##     data = udaje_2023)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.5790 -0.6944 -0.1406  0.7645  4.0497 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            9.401e-01  1.026e-01   9.167   <2e-16 ***
## X2023.population       1.567e-09  2.501e-09   0.627   0.5316    
## area..km..             2.591e-07  1.435e-07   1.805   0.0723 .  
## density..km..         -2.317e-04  1.436e-04  -1.614   0.1080    
## I(X2023.population^2) -1.582e-18  1.689e-18  -0.937   0.3499    
## I(area..km..^2)       -2.302e-14  1.001e-14  -2.300   0.0224 *  
## I(density..km..^2)     1.075e-08  7.613e-09   1.412   0.1593    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.219 on 227 degrees of freedom
## Multiple R-squared:  0.05012,    Adjusted R-squared:  0.02501 
## F-statistic: 1.996 on 6 and 227 DF,  p-value: 0.0672

# Porovnanie pôvodného a rozšíreného modelu ANOVA testom
anova(model_lin, model_rozsireny)

# Rozšírený RESET test
library(lmtest)
resettest(model_rozsireny)

## 
##  RESET test
## 
## data:  model_rozsireny
## RESET = 5.7269, df1 = 2, df2 = 225, p-value = 0.00375

5. Transformácia pomocou dummy premennej a lineárnej lomenej funkcie

# Zavedenie dummy premennej DUM podľa zvoleného prahu
udaje_2023$DUM <- ifelse(udaje_2023$X2023.population < 1e7, 0, 1)

# Model so zlomom v autonómnom člene
modelD_auto <- lm(growth.rate ~ +1 + DUM + X2023.population + area..km.. + density..km..,
                  data = udaje_2023)
summary(modelD_auto)

## 
## Call:
## lm(formula = growth.rate ~ +1 + DUM + X2023.population + area..km.. + 
##     density..km.., data = udaje_2023)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.8198 -0.6671 -0.0756  0.7208  4.2177 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       7.689e-01  1.058e-01   7.270 5.63e-12 ***
## DUM               6.513e-01  1.730e-01   3.764 0.000212 ***
## X2023.population -5.454e-10  6.533e-10  -0.835 0.404648    
## area..km..       -4.722e-08  5.168e-08  -0.914 0.361814    
## density..km..    -2.980e-05  4.024e-05  -0.741 0.459714    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.205 on 229 degrees of freedom
## Multiple R-squared:  0.06382,    Adjusted R-squared:  0.04747 
## F-statistic: 3.903 on 4 and 229 DF,  p-value: 0.004366

# Model so zlomom v sklone nadroviny (interakcia X2023.population a DUM)
modelD_sklon <- lm(growth.rate ~ +1 + X2023.population + I(DUM*X2023.population) +
                      area..km.. + density..km..,
                    data = udaje_2023)
summary(modelD_sklon)

## 
## Call:
## lm(formula = growth.rate ~ +1 + X2023.population + I(DUM * X2023.population) + 
##     area..km.. + density..km.., data = udaje_2023)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.4614 -0.7358 -0.1598  0.7289  4.0130 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                1.024e+00  1.016e-01  10.071   <2e-16 ***
## X2023.population          -1.498e-08  3.435e-08  -0.436    0.663    
## I(DUM * X2023.population)  1.481e-08  3.430e-08   0.432    0.666    
## area..km..                -5.222e-09  5.203e-08  -0.100    0.920    
## density..km..             -4.672e-05  4.118e-05  -1.134    0.258    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.242 on 229 degrees of freedom
## Multiple R-squared:  0.006704,   Adjusted R-squared:  -0.01065 
## F-statistic: 0.3864 on 4 and 229 DF,  p-value: 0.8183

# Porovnanie pôvodného lineárneho modelu s modelom so zlomom
anova(model_lin, modelD_sklon)

# Rozšírený RESET test pre model so zlomom
library(lmtest)
resettest(modelD_sklon)

## 
##  RESET test
## 
## data:  modelD_sklon
## RESET = 1.5274, df1 = 2, df2 = 227, p-value = 0.2193

Záver

Nelineárne špecifikácie

Najprv som odhadla lineárny model, kde som rast populácie (growth.rate) vysvetľovala pomocou premenných X2023.population, area..km.. a density..km…

Výsledky ukazujú, že žiadna z vysvetľujúcich premenných nebola štatisticky významná (p > 0,05 okrem interceptu) a model vysvetľuje len veľmi malú časť variability rastu populácie (R^2 ≈ 0,006).

RESET test

Na overenie správnej špecifikácie som použila Ramseyho RESET test:

Lineárny model: p = 0,294 a 0,1181

Logaritmický model: p = 0,2055 a 0,4791

Interpretácia: p-hodnota > 0,05 znamená, že neexistuje štatisticky významný dôkaz nesprávnej špecifikácie. Funkčná forma modelu je vhodná, netreba pridávať ďalšie nelineárne členy ani nové premenné.

Kvadratická rozšírená špecifikácia

Zaviedla som kvadratické členy pre X2023.population a density..km… Tento modifikovaný model mierne zvýšil R² a znížil reziduálnu sumu štvorcov. Porovnanie s pôvodným lineárnym modelom pomocou ANOVA testu: p ≈ 0,077 → efekt nie je silne štatisticky významný. RESET test: p = 0,324 → stále neexistuje dôkaz zásadnej nesprávnej špecifikácie.

Rozšírený model s viacerými kvadratickými členmi

Pridala som kvadratické členy pre všetky tri premenné (X2023.population^2, area..km..^2, density..km..^2).

Upravený koeficient determinácie sa mierne zvýšil (R^2_adj ≈ 0,025).

ANOVA test vs základný model: p = 0,0158 → zavedenie kvadratických členov významne zlepšilo model.

Transformácia pomocou dummy premennej

Zaviedla som dummy premennú DUM podľa prahu X2023.population < 1e7.

Model so zlomom v interceptoch (modelD_auto)

– Zlom nebol štatisticky významný pri vysvetľovaní rastu populácie.

– R² sa mierne zvýšil.

Model so zlomom v sklone (modelD_sklon)

– Interakcia DUM * X2023.population nebola štatisticky významná (p ≈ 0,666).

– ANOVA porovnanie s pôvodným lineárnym modelom: p = 0,666 → pridanie zlomu nezlepšilo model.

– RESET test: p = 0,219 → funkčná forma modelu zostáva vhodná.

Záverečné zhodnotenie

Lineárny aj logaritmický model sú podľa RESET testu štatisticky vhodné a funkčná forma modelu je primeraná.

Kvadratické členy mierne zlepšili model, najmä pri premenných area a density.

Transformácia pomocou dummy premennej nebola potrebná, nepriniesla významné zlepšenie modelu.

Celkovo je jednoduchý lineárny alebo logaritmický model dostatočný na opis rastu populácie podľa týchto premenných.

Ekonometria - cv. 4 + 5 + 6 + 7

Bartáková Terézia

Úvod do problému

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

Boxploty

Lineárna regresia

Diagnostické grafy regresného modelu

Výsledky diagnostiky

Residuals vs. fitted

Q-Q plot

Scale-Location

Residuals vs Leverage

Záver

Heteroskedasticita

Zlogaritmované grafy

Zlogaritmované grafy bez odľahlých hodnôt

Interpretácia skúmania heteroskedasticity

Nelineárne špecifikácie

RESET test pre základný model

RESET test pre logaritmovaný model

Interpretácia RESET testu

1. Residuals vs Fitted

2. Component + Residual (C+R grafy)

3. Nelineárna špecifikácia - porovnanie základného a modifikovaného modelu

4. Rozšírený RESET test pre nelineárny model

5. Transformácia pomocou dummy premennej a lineárnej lomenej funkcie

Záver