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

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

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.