Práca s databázou - import údajov, grafy, štatistiky

Importovanie vlastnej databázy

head(world_population)                            # niekolko prvych riadkov

colnames(world_population)

 [1] "Rank"                        "CCA3"                       
 [3] "Country.Territory"           "Capital"                    
 [5] "Continent"                   "X2022.Population"           
 [7] "X2020.Population"            "X2015.Population"           
 [9] "X2010.Population"            "X2000.Population"           
[11] "X1990.Population"            "X1980.Population"           
[13] "X1970.Population"            "Area..km.."                 
[15] "Density..per.km.."           "Growth.Rate"                
[17] "World.Population.Percentage"

Grafy

ggplot2 - knižnica pre grafy

Výber a následné triedenie

library(ggplot2)
SVK_CZE_data <- world_population %>%
  filter(CCA3 %in% c("SVK", "CZE")) %>%  # Filtrovanie na Slovensko a Česko
  select(CCA3, `2020 Population`, `2010 Population`, `2000 Population`, `1990 Population`, `1980 Population`) # Vyberáme len relevantné stĺpce

Error in world_population %>% filter(CCA3 %in% c("SVK", "CZE")) %>% select(CCA3,  : 
  could not find function "%>%"

Scatter plot

library(dplyr)
library(tidyr)
library(ggplot2)

# Predpokladám, že dataset je uložený v premennej `SVK_CZE_data`

# 1. Vyberieme potrebné stĺpce a dáta len pre Slovensko a Česko
SVK_CZE_sel <- SVK_CZE_data %>%
  filter(CCA3 %in% c("SVK", "CZE")) %>%  # Filtrovanie na Slovensko a Česko
  select(CCA3, `2020 Population`, `2010 Population`, `2000 Population`, `1990 Population`, `1980 Population`)  # Vyberáme len relevantné stĺpce

# 2. Premeníme dáta do dlhého formátu (pre graf)
SVK_CZE_long <- SVK_CZE_sel %>%
  pivot_longer(cols = c(`2020 Population`, `2010 Population`, `2000 Population`, `1990 Population`, `1980 Population`),
               names_to = "year",     # Stĺpec s rokmi
               values_to = "population")  # Stĺpec s populáciou

# 3. Vytvoríme stĺpcový graf - Zobrazíme populáciu Slovenska a Česka v rokoch 1980, 1990, 2000, 2010 a 2020
ggplot(SVK_CZE_long, aes(x = year,y = population, fill = CCA3, group = CCA3)) + geom_bar(stat = "identity", position = "dodge", width = 0.8) + # Stĺpcový graf s oddelenými stĺpcami pre každú krajinu 
  geom_text(aes(label = population), position = position_dodge(width = 0.3), vjust = -0.1) + # Pridanie označení hodnôt 
  theme_minimal() + # Minimalistický dizajn grafu 
  labs(title = "Vývoj populácie Slovenska a Česka", # Upravený nadpis 
       x = "Rok", 
       y = "Populácia", 
       fill = "Krajina") + # Popis farieb podľa krajiny 
  theme(axis.text.x = element_text(angle = 30, hjust = 1)) # Rotácia textu na osi X

Boxplot

# Bar plot with grouping
library(ggplot2)

ggplot(SVK_CZE_long, aes(x = year, y = population, color = CCA3, group = CCA3)) +
  geom_line(linewidth = 1.2) +                 # Čiary pre každú krajinu
  geom_point(size = 3) +                       # Bodky na čiarach
  geom_text(aes(label = population),           # Hodnoty nad bodkami
            vjust = -0.5, size = 2,            # Malé písmo, posunuté mierne nahor
            show.legend = FALSE) +             # Bez legendy pre text
  theme_minimal() +
  labs(
    title = "Vývoj populácie Slovenska a Česka",
    x = "Rok",
    y = "Populácia",
    color = "Krajina"
  ) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

Základné štatistiky

knitr - tabuľka

library(dplyr)
library(knitr)

# 1. Vyberieme potrebné stĺpce a dáta len pre Slovensko a Česko
SVK_CZE_sel <- SVK_CZE_data %>%
  filter(CCA3 %in% c("SVK", "CZE")) %>%  # Filtrovanie na Slovensko a Česko
  select(CCA3, `2020 Population`, `2010 Population`, `2000 Population`, `1990 Population`, `1980 Population`)  # Relevantné stĺpce

# 2. Premeníme dáta do dlhého formátu
SVK_CZE_long <- SVK_CZE_sel %>%
  pivot_longer(cols = c(`2020 Population`, `2010 Population`, `2000 Population`, `1990 Population`, `1980 Population`),
               names_to = "year",
               values_to = "population")
# 3. Vypočítame základné štatistiky populácie podľa krajiny a roku
SVK_CZE_stats <- SVK_CZE_long %>%
  group_by(CCA3, year) %>%
  summarise(
    n      = n(),
    mean   = mean(population, na.rm = TRUE),
    sd     = sd(population, na.rm = TRUE),
    min    = min(population, na.rm = TRUE),
    q25    = quantile(population, 0.25, na.rm = TRUE),
    median = median(population, na.rm = TRUE),
    q75    = quantile(population, 0.75, na.rm = TRUE),
    max    = max(population, na.rm = TRUE),
    .groups = "drop"
  )

# 4. Vytvoríme tabuľku pomocou knitr
kable(SVK_CZE_stats, digits = 0, caption = "Základné štatistiky populácie Slovenska a Česka (1980–2020)")

Základné štatistiky populácie Slovenska a Česka (1980–2020)

CCA3	year	n	mean	sd	min	q25	median	q75	max
CZE	1980 Population	1	10270060	NA	10270060	10270060	10270060	10270060	10270060
CZE	1990 Population	1	10301192	NA	10301192	10301192	10301192	10301192	10301192
CZE	2000 Population	1	10234710	NA	10234710	10234710	10234710	10234710	10234710
CZE	2010 Population	1	10464749	NA	10464749	10464749	10464749	10464749	10464749
CZE	2020 Population	1	10530953	NA	10530953	10530953	10530953	10530953	10530953
SVK	1980 Population	1	4973883	NA	4973883	4973883	4973883	4973883	4973883
SVK	1990 Population	1	5261305	NA	5261305	5261305	5261305	5261305	5261305
SVK	2000 Population	1	5376690	NA	5376690	5376690	5376690	5376690	5376690
SVK	2010 Population	1	5396424	NA	5396424	5396424	5396424	5396424	5396424
SVK	2020 Population	1	5456681	NA	5456681	5456681	5456681	5456681	5456681

Testovanie hypotéz

t-test: Porovnanie populácie Slovenska a Česka v rokoch 2020

# 1. Vyberieme populáciu Slovenska a Česka v roku 2020
SVK_CZE_2020 <- SVK_CZE_data %>%
  filter(CCA3 %in% c("SVK", "CZE")) %>%
  select(CCA3, `2020 Population`)

# 2. Vykonáme t-test na porovnanie populácie medzi SVK a CZE v roku 2020
SVK_pop <- SVK_CZE_long %>% filter(CCA3 == "SVK") %>% pull(population)
CZE_pop <- SVK_CZE_long %>% filter(CCA3 == "CZE") %>% pull(population)

# 3. Výpis výsledku testu
t_test_result <- t.test(SVK_pop, CZE_pop)
print(t_test_result)

    Welch Two Sample t-test

data:  SVK_pop and CZE_pop
t = -48.899, df = 7.0276, p-value = 3.655e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -5312182 -4822491
sample estimates:
mean of x mean of y 
  5292997  10360333 

ANOVA: Comparing Reading Scores Across Programs

anova_result <- aov(population ~ year, data = SVK_CZE_long)
summary(anova_result)

            Df    Sum Sq   Mean Sq F value Pr(>F)
year         4 1.663e+11 4.157e+10   0.003      1
Residuals    5 6.424e+13 1.285e+13               

Linear Regression: Predicting Math Scores

# 3. Vytvoríme lineárny model populácie podľa roku a krajiny
model <- lm(population ~ year + CCA3, data = SVK_CZE_long)
summary(model)

Call:
lm(formula = population ~ year + CCA3, data = SVK_CZE_long)

Residuals:
        1         2         3         4         5         6         7 
   3467.9     494.4 -104658.1  -13724.6  114420.4   -3467.9    -494.4 
        8         9        10 
 104658.1   13724.6 -114420.4 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         10155640      85286 119.077 2.98e-08 ***
year1990 Population   159277     110104   1.447   0.2216    
year2000 Population   183728     110104   1.669   0.1705    
year2010 Population   308615     110104   2.803   0.0487 *  
year2020 Population   371846     110104   3.377   0.0279 *  
CCA3SVK             -5067336      69636 -72.769 2.14e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 110100 on 4 degrees of freedom
Multiple R-squared:  0.9992,    Adjusted R-squared:  0.9983 
F-statistic:  1062 on 5 and 4 DF,  p-value: 2.478e-06

# install.packages(c("broom", "kableExtra", "dplyr", "stringr"))
library(broom)
library(dplyr)
library(kableExtra)
library(stringr)

# Your model (already fitted)
# model <- lm(ESG.INDEX ~ RETURN.ON.ASSETS + FIRM.SIZE + DEBT.TO.ASSET, data = udaje.2013)

coef.tbl <- tidy(model, conf.int = TRUE) %>%
  mutate(
    term = recode(term,
      "(Intercept)" = "Intercept",
      "year1990 Population" = "Rok 1990",
      "year2000 Population" = "Rok 2000",
      "year2010 Population" = "Rok 2010",
      "year2020 Population" = "Rok 2020",
      "CCA3SVK" = "Slovensko"
    ),
    stars = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      p.value < 0.1   ~ "·",
      TRUE            ~ ""
    )
  ) %>%
  transmute(
    Term = term,
    Estimate = estimate,
    `Std. Error` = std.error,
    `t value` = statistic,
    `p value` = p.value,
    `95% CI` = str_c("[", round(conf.low, 0), ", ", round(conf.high, 0), "]"),
    Sig = stars
  )

coef.tbl %>%
  kable(
    digits = 0,
    caption = "OLS regresné koeficienty (population ~ year + CCA3)"
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed")) %>%
  column_spec(1, bold = TRUE) %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2") %>%
  footnote(
    general = "Signif. codes: *** p<0.001, ** p<0.01, * p<0.05, · p<0.1.",
    threeparttable = TRUE)

OLS regresné koeficienty (population ~ year + CCA3)

Term	Estimate	Std. Error	t value	p value	95% CI	Sig
Intercept	10155640	85287	119	0	[9918846, 10392433]	***
Rok 1990	159277	110104	1	0	[-146422, 464976]
Rok 2000	183728	110104	2	0	[-121970, 489427]
Rok 2010	308615	110104	3	0	[2916, 614314]	*
Rok 2020	371845	110104	3	0	[66147, 677544]	*
Slovensko	-5067336	69636	-73	0	[-5260677, -4873995]	***
Note:
Signif. codes: *** p<0.001, ** p<0.01, * p<0.05, · p<0.1.

NA

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

# Načítanie dát
udaje <- read.csv("world_population.csv", dec=".", sep=",", header = TRUE)

# Kontrola názvov stĺpcov
names(udaje)

 [1] "Rank"                        "CCA3"                       
 [3] "Country.Territory"           "Capital"                    
 [5] "Continent"                   "X2022.Population"           
 [7] "X2020.Population"            "X2015.Population"           
 [9] "X2010.Population"            "X2000.Population"           
[11] "X1990.Population"            "X1980.Population"           
[13] "X1970.Population"            "Area..km.."                 
[15] "Density..per.km.."           "Growth.Rate"                
[17] "World.Population.Percentage"

# Výber len číselných údajov (uprav podľa presných názvov stĺpcov)
udaje.2020 <- udaje[, c("X2020.Population", "Area..km..", "Density..per.km..")]

# Konverzia na numeric
udaje.2020 <- as.data.frame(lapply(udaje.2020, as.numeric))

# Nahradenie NA mediánmi
column_medians <- sapply(udaje.2020, median, na.rm = TRUE)
for (col in names(udaje.2020)) {
  udaje.2020[[col]][is.na(udaje.2020[[col]])] <- column_medians[col]
}

# Vykreslenie boxplotov
par(mfrow = c(2, 2))
par(mar = c(4, 4, 2, 1))
for (col in names(udaje.2020)) {
  boxplot(udaje.2020[[col]], main = col, xlab = "Hodnoty", col = "lightpink", border = "orange")
}

Boxploty – Populácia, Rozloha a Hustota obyvateľstva

Interpretácia grafov

1. X2020.Population

Väčšina hodnôt sa koncentruje pri veľmi nízkych hodnotách populácie.

Viditeľných je niekoľko výrazných odľahlých hodnôt (outliers) – niektoré krajiny/územia majú populáciu výrazne vyššiu než ostatné.

Rozdelenie je silne pravostranné (pozitívne zošikmené) – väčšina pozorovaní má malé hodnoty, zatiaľ čo niekoľko extrémnych pozorovaní ťahá medián smerom nadol.

2. Area..km..

Väčšina krajín má pomerne malú rozlohu, pričom zopár veľmi veľkých území spôsobuje silnú pravostrannú šikmosť.

Viaceré body ďaleko od hlavného poľa dát naznačujú odľahlé pozorovania (veľké krajiny).

Podobne ako pri populácii, dáta sú veľmi nesymetrické – odporúča sa zvážiť transformáciu (napr. log(area)).

3. Density..per.km..

Hustota obyvateľstva má väčšinu hodnôt sústredených pri nízkych číslach, no niekoľko oblastí má extrémne vysokú hustotu.

To sa prejavuje odľahlými bodmi nad horným fúzom boxplotu.

Rozdelenie je pravostranné, teda existujú niektoré extrémne husto obývané regióny.

Transformácia (napr. log) by mohla pomôcť zlepšiť symetriu rozdelenia.

Lineárna regresia

#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 = population ~ year + CCA3, data = SVK_CZE_long)

Residuals:
        1         2         3         4         5         6         7         8         9        10 
   3467.9     494.4 -104658.1  -13724.6  114420.4   -3467.9    -494.4  104658.1   13724.6 -114420.4 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         10155640      85286 119.077 2.98e-08 ***
year1990 Population   159277     110104   1.447   0.2216    
year2000 Population   183728     110104   1.669   0.1705    
year2010 Population   308615     110104   2.803   0.0487 *  
year2020 Population   371846     110104   3.377   0.0279 *  
CCA3SVK             -5067336      69636 -72.769 2.14e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 110100 on 4 degrees of freedom
Multiple R-squared:  0.9992,    Adjusted R-squared:  0.9983 
F-statistic:  1062 on 5 and 4 DF,  p-value: 2.478e-06

# Nastaviť rozloženie 2 x 2
par(mfrow = c(2, 2))

# Vykresliť všetky 4 diagnostické grafy modelu
plot(model)

# (Voliteľné) pridať spoločný nadpis
mtext("Diagnostické grafy regresného modelu", outer = TRUE, cex = 1.2, font = 2)

# Resetovať layout
par(mfrow = c(1, 1))

Interpretácia grafov

Residuals vs Fitted

Reziduály sú rozložené okolo nuly bez výrazného trendu, čo naznačuje dobré prispôsobenie modelu. Mierne zakrivenie čiary poukazuje na možnú nelinearitu alebo vplyv niekoľkých odľahlých pozorovaní.

Normal Q–Q plot

Väčšina bodov leží blízko diagonály, takže rozdelenie rezíduí je približne normálne. Menšie odchýlky na okrajoch naznačujú len mierne porušenie normality.

Scale–Location

Rozptyl rezíduí sa javí ako približne konštantný, čo podporuje predpoklad homoskedasticity. Mierne výkyvy na začiatku môžu byť spôsobené niekoľkými extrémnymi bodmi.

Residuals vs Factor Levels

Reziduály sa v jednotlivých faktorových úrovniach pohybujú okolo nuly, bez systematického trendu. Body 5 a 10 pôsobia ako možné vplyvné pozorovania, ktoré môžu ovplyvňovať model.

install.packages("tseries")  # len raz
library(tseries)   

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

    Jarque Bera Test

data:  residuals
X-squared = 0.1133, df = 2, p-value = 0.9449

# outlier test (see p-value for Bonferroni correction)
install.packages("car")  

Error in install.packages : Updating loaded packages

library(car)  
outlierTest(model)

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

model2 <- lm(Density..per.km.. ~ X2020.Population + Area..km.., data = udaje.2020)
summary(model2)

Call:
lm(formula = Density..per.km.. ~ X2020.Population + Area..km.., 
    data = udaje.2020)

Residuals:
    Min      1Q  Median      3Q     Max 
 -493.6  -417.8  -361.8  -217.8 22677.4 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       4.948e+02  1.434e+02   3.452 0.000662 ***
X2020.Population  2.622e-08  1.124e-06   0.023 0.981404    
Area..km..       -7.495e-05  8.647e-05  -0.867 0.386995    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2071 on 231 degrees of freedom
Multiple R-squared:  0.003987,  Adjusted R-squared:  -0.004636 
F-statistic: 0.4624 on 2 and 231 DF,  p-value: 0.6304

# Diagnostické grafy
par(mfrow = c(2, 2))
par(mar = c(4, 4, 2, 1))
plot(model2)
par(mfrow = c(1, 1))

Interpretácia grafov

Residuals vs Fitted

Reziduály sa väčšinou pohybujú okolo nuly, čo naznačuje, že model nemá výrazné systematické chyby v predikciách. Niekoľko bodov s veľkými reziduálmi (napr. 120, 135, 188) predstavuje možné odľahlé pozorovania.

Normal Q–Q plot

Väčšina bodov sleduje diagonálu, no koniec grafu sa výrazne odchyľuje smerom nahor. To znamená, že reziduály nie sú dokonale normálne rozdelené a existujú extrémne hodnoty v pravej časti rozdelenia.

Scale–Location

Rozptyl rezíduí sa javí ako relatívne konštantný, bez výrazného rozširovania alebo zužovania. Len niekoľko bodov s vysokými štandardizovanými reziduálmi môže naznačovať miernu heteroskedasticitu.

Residuals vs Leverage

Väčšina bodov má nízke hodnoty leverage, čo je dobré, no niekoľko pozorovaní (napr. 120, 135, 172) má väčší vplyv podľa Cookovej vzdialenosti. Tieto pozorovania môžu výraznejšie ovplyvňovať odhad parametrov a je vhodné ich skontrolovať.

residuals <- residuals(model)
jarque.bera.test(residuals)
outlierTest(model)

Zadanie 6 - Heteroskedasticita

library(dplyr)
library(tidyr)
library(ggplot2)
library(patchwork)

# 1️⃣ Preformátovanie dát do dlhého formátu
world_population_long <- world_population %>%
  select(Continent, Country.Territory, 
         X2022.Population, X2020.Population, X2015.Population, 
         X2010.Population, X2000.Population, X1990.Population, X1980.Population, X1970.Population) %>%
  pivot_longer(cols = starts_with("X"), 
               names_to = "Year", 
               values_to = "Population") %>%
  mutate(Year = as.numeric(gsub("X|\\.Population", "", Year))) %>%  # odstráni "X" a ".Population", zostáva rok
  group_by(Country.Territory) %>%
  arrange(Year) %>%
  mutate(GrowthRate = (Population / lag(Population) - 1) * 100) %>%
  ungroup() %>%
  filter(!is.na(GrowthRate))

# 2️⃣ Lineárny model
model <- lm(GrowthRate ~ Year, data = world_population_long)

# 3️⃣ Predikcie a reziduá
world_population_long$pred_model <- predict(model)
world_population_long$resid_model <- resid(model)

# 4️⃣ Grafy heteroskedasticity
p1 <- ggplot(world_population_long, aes(x = Year, y = resid_model^2)) +
  geom_point(alpha = 1, color = "steelblue") +
  geom_smooth(method = "loess", se = FALSE, color = "brown") +
  labs(x = "Rok", y = "Štvorce rezíduí", title = "Residuals² vs Rok") +
  theme_minimal()

p2 <- ggplot(world_population_long, aes(x = GrowthRate, y = resid_model^2)) +
  geom_point(alpha = 0.6, color = "lightgreen") +
  geom_smooth(method = "loess", se = FALSE, color = "pink") +
  labs(x = "Tempo rastu populácie (%)", y = "Štvorce rezíduí", title = "Residuals² vs Growth Rate") +
  theme_minimal()

# 5️⃣ Zobrazenie oboch grafov vedľa seba
p1 + p2

library(ggplot2)
library(patchwork)

model <- lm(GrowthRate ~ Continent, data = world_population_long)
world_population_long$resid_model <- resid(model)

p1 <- ggplot(world_population_long, aes(x = Year, y = resid_model^2)) +
  geom_point(alpha = 0.6, color = "grey") +
  geom_smooth(method = "loess", se = FALSE, color = "purple") +
  labs(x = "Rok", y = "Štvorce rezíduí", title = "Residuals² vs Rok") +
  theme_minimal() +
  theme(
  plot.background = element_rect(fill = "#f5f5dc", color = NA),  # béžová
  panel.background = element_rect(fill = "#f5f5dc", color = NA)
) 

p2 <- ggplot(world_population_long, aes(x = GrowthRate, y = resid_model^2)) +
  geom_point(alpha = 0.6, color = "darkgreen") +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Tempo rastu populácie (%)", y = "Štvorce rezíduí", title = "Residuals² vs Growth Rate") +
  theme_minimal() +
  theme(
  plot.background = element_rect(fill = "#f5f5dc", color = NA),  # béžová
  panel.background = element_rect(fill = "#f5f5dc", color = NA)
) 
  

p1 + p2

Graf: Štvorce rezíduí vs Rok

V tomto grafe sledujeme, či sa veľkosť chýb modelu mení v čase. Fialová vyhladená krivka je takmer rovná, čo naznačuje, že rozptyl rezíduí je stabilný naprieč rokmi. To znamená, že premenná „rok“ pravdepodobne nespôsobuje heteroskedasticitu. Model teda v tomto smere neukazuje žiadny systematický problém.

Graf: Štvorce rezíduí vs Tempo rastu populácie

Tu skúmame, či veľkosť chýb závisí od tempa rastu populácie. Červená krivka má výrazný oblúkovitý tvar, čo naznačuje, že model robí väčšie chyby pri extrémnych hodnotách rastu. To je vizuálny náznak heteroskedasticity – teda toho, že rozptyl rezíduí nie je konštantný. Tento problém potvrdzuje aj výsledok Breusch–Paganovho testu.

# 📦 Načítanie knižníc
library(dplyr)
library(tidyr)
library(ggplot2)
library(lmtest)
library(patchwork)

# 1️⃣ Preformátovanie dát do dlhého formátu
world_population_long <- world_population %>%
  select(Continent, Country.Territory, 
         X2022.Population, X2020.Population, X2015.Population, 
         X2010.Population, X2000.Population, X1990.Population, X1980.Population, X1970.Population) %>%
  pivot_longer(cols = starts_with("X"), 
               names_to = "Year", 
               values_to = "Population") %>%
  mutate(Year = as.numeric(gsub("X|\\.Population", "", Year))) %>%  # z názvu stĺpca vyber rok
  group_by(Country.Territory) %>%
  arrange(Year) %>%
  mutate(GrowthRate = (Population / lag(Population) - 1) * 100) %>%
  ungroup() %>%
  filter(!is.na(GrowthRate))

# 2️⃣ Vytvorenie lineárneho modelu
model <- lm(GrowthRate ~ Year, data = world_population_long)
summary(model)

Call:
lm(formula = GrowthRate ~ Year, data = world_population_long)

Residuals:
    Min      1Q  Median      3Q     Max 
-69.121  -7.491  -2.093   5.668 239.690 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1087.20768   59.39472   18.30   <2e-16 ***
Year          -0.53527    0.02962  -18.07   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 17.63 on 1636 degrees of freedom
Multiple R-squared:  0.1664,    Adjusted R-squared:  0.1659 
F-statistic: 326.6 on 1 and 1636 DF,  p-value: < 2.2e-16

# 3️⃣ Breusch-Pagan test na heteroskedasticitu
bptest(model)

    studentized Breusch-Pagan test

data:  model
BP = 19.018, df = 1, p-value = 1.295e-05

Breusch–Paganov test overuje, či sa veľkosť chýb modelu mení v závislosti od vysvetľujúcich premenných – teda či je prítomná heteroskedasticita. Výsledok testu (BP = 19.018, p-hodnota < 0.001) naznačuje, že rozptyl rezíduí nie je konštantný. To znamená, že model robí systematicky väčšie chyby pri určitých hodnotách vstupných premenných. Takýto problém môže ovplyvniť spoľahlivosť odhadov a vyžaduje úpravu modelu.

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

# Load the package
library(lmtest)

model2 <- lm(GrowthRate ~ +1 + I(log(Year)),data=world_population_long)
summary(model2)

Call:
lm(formula = GrowthRate ~ +1 + I(log(Year)), data = world_population_long)

Residuals:
    Min      1Q  Median      3Q     Max 
-69.089  -7.493  -2.125   5.670 239.708 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8153.92     450.83   18.09   <2e-16 ***
I(log(Year)) -1070.57      59.29  -18.06   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 17.63 on 1636 degrees of freedom
Multiple R-squared:  0.1662,    Adjusted R-squared:  0.1657 
F-statistic:   326 on 1 and 1636 DF,  p-value: < 2.2e-16

# Run the Breusch–Pagan test
bptest(model2)

    studentized Breusch-Pagan test

data:  model2
BP = 18.994, df = 1, p-value = 1.311e-05

Tento model skúma vzťah medzi rastom populácie a logaritmicky transformovaným rokom. Výsledky ukazujú, že log(Year) má silne negatívny a štatisticky významný vplyv na rast populácie. Napriek transformácii však Breusch–Paganov test (BP = 18.994, p < 0.001) stále naznačuje prítomnosť heteroskedasticity. To znamená, že rozptyl chýb sa mení a model by mohol byť ďalej upravený, napríklad použitím robustných metód.

#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) 1087.207679   66.659480  16.310 < 2.2e-16 ***
Year          -0.535267    0.033097 -16.173 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1