Cvičenie 5 – Základná diagnostika regresného modelu

Moja vysnívaná kariéra: Dátový vedec

Autor

Martin Bakič

Publikované

9. novembra 2026

knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE
)

1 Úvod

V tomto cvičení pokračujem v analýze svojej databázy Latest Data Science Salaries.

Moja hlavná téma je Moja vysnívaná kariéra: Dátový vedec.

V predchádzajúcom cvičení som vytvoril lineárny regresný model, ktorým som skúmal výšku platu pracovnej pozície Data Scientist.

V tomto cvičení sa zameriam najmä na:

uloženie a obnovenie workspace,
interpretáciu regresných koeficientov a t-testov,
F-test celého regresného modelu,
koeficient determinácie (R^2),
grafickú analýzu rezíduí,
test normality rezíduí,
identifikáciu odľahlých a vplyvných pozorovaní,
modifikáciu pôvodného regresného modelu.

2 Príprava databázy

Databázu mám uloženú v priečinku údaje.

cesta <- "údaje/Latest_Data_Science_Salaries.csv"

file.exists(cesta)

[1] TRUE

Výsledkom by mala byť hodnota TRUE.

2.1 Načítanie údajov

udaje <- read.csv(
  cesta,
  sep = ",",
  dec = ".",
  na.strings = c("", "NA"),
  check.names = FALSE,
  stringsAsFactors = FALSE
)

head(udaje)

       Job Title Employment Type Experience Level Expertise Level Salary
1  Data Engineer       Full-Time           Senior          Expert 210000
2  Data Engineer       Full-Time           Senior          Expert 165000
3  Data Engineer       Full-Time           Senior          Expert 185900
4  Data Engineer       Full-Time           Senior          Expert 129300
5 Data Scientist       Full-Time           Senior          Expert 140000
6 Data Scientist       Full-Time           Senior          Expert 126000
       Salary Currency Company Location Salary in USD Employee Residence
1 United States Dollar    United States        210000      United States
2 United States Dollar    United States        165000      United States
3 United States Dollar    United States        185900      United States
4 United States Dollar    United States        129300      United States
5 United States Dollar    United States        140000      United States
6 United States Dollar    United States        126000      United States
  Company Size Year
1       Medium 2023
2       Medium 2023
3       Medium 2023
4       Medium 2023
5       Medium 2023
6       Medium 2023

Názvy premenných upravím do formátu vhodného pre R.

names(udaje) <- make.names(
  trimws(
    names(udaje)
  ),
  unique = TRUE
)

names(udaje)

 [1] "Job.Title"          "Employment.Type"    "Experience.Level"  
 [4] "Expertise.Level"    "Salary"             "Salary.Currency"   
 [7] "Company.Location"   "Salary.in.USD"      "Employee.Residence"
[10] "Company.Size"       "Year"

2.2 Výber pracovnej pozície Data Scientist

Keďže mojou vysnívanou kariérou je pozícia Data Scientist, vyberiem iba tieto záznamy.

data_scientist <- udaje[
  udaje$Job.Title == "Data Scientist",
]

Vyberiem premenné potrebné pre regresnú analýzu.

data_scientist <- data_scientist[, c(
  "Salary.in.USD",
  "Year",
  "Experience.Level",
  "Company.Size"
)]

Odstránim prípadné neúplné riadky.

data_scientist <- data_scientist[
  complete.cases(data_scientist),
]

Počet použitých pozorovaní:

nrow(data_scientist)

[1] 635

2.3 Úprava premenných

Úroveň skúseností nastavím ako faktor.

Ako základnú kategóriu použijem Entry.

data_scientist$Experience.Level <- factor(
  data_scientist$Experience.Level,
  levels = c(
    "Entry",
    "Mid",
    "Senior",
    "Executive"
  )
)

Veľkosť firmy nastavím ako faktor.

Ako základnú kategóriu použijem Small.

data_scientist$Company.Size <- factor(
  data_scientist$Company.Size,
  levels = c(
    "Small",
    "Medium",
    "Large"
  )
)

Pre jednoduchšiu interpretáciu roku vytvorím premennú, ktorá vyjadruje počet rokov od roku 2020.

data_scientist$Rok.od.2020 <-
  data_scientist$Year - 2020

Skontrolujem výslednú databázu.

head(data_scientist)

   Salary.in.USD Year Experience.Level Company.Size Rok.od.2020
5         140000 2023           Senior       Medium           3
6         126000 2023           Senior       Medium           3
7         170000 2023           Senior       Medium           3
8         130000 2023           Senior       Medium           3
17        185900 2023           Senior       Medium           3
18        129300 2023           Senior       Medium           3

str(data_scientist)

'data.frame':   635 obs. of  5 variables:
 $ Salary.in.USD   : int  140000 126000 170000 130000 185900 129300 225000 156400 61520 49216 ...
 $ Year            : int  2023 2023 2023 2023 2023 2023 2023 2023 2023 2023 ...
 $ Experience.Level: Factor w/ 4 levels "Entry","Mid",..: 3 3 3 3 3 3 3 3 2 2 ...
 $ Company.Size    : Factor w/ 3 levels "Small","Medium",..: 2 2 2 2 2 2 2 2 2 2 ...
 $ Rok.od.2020     : num  3 3 3 3 3 3 3 3 3 3 ...

3 Ukladanie a obnovovanie workspace

Workspace predstavuje objekty, ktoré máme aktuálne uložené v pracovnej pamäti R.

Ich zoznam zobrazím príkazom:

ls()

[1] "cesta"          "data_scientist" "udaje"

3.1 Uloženie workspace

Celý aktuálny workspace uložím do súboru:

save.image(
  "pamat_Cvicenie5.RData"
)

Takto môžem uložiť databázu aj ostatné vytvorené objekty a v ďalšom cvičení ich nemusím vytvárať úplne od začiatku.

3.2 Mazanie vybraných objektov

Vytvorím dva skúšobné objekty.

aaa <- 6
bbb <- "text"

aaa

[1] 6

bbb

[1] "text"

Vymažem iba tieto dva objekty.

rm(
  list = c(
    "aaa",
    "bbb"
  )
)

ls()

[1] "cesta"          "data_scientist" "udaje"

3.3 Vyčistenie celého workspace

Pred úplným vyčistením ešte raz uložím aktuálnu pamäť.

save.image(
  "pamat_Cvicenie5.RData"
)

Teraz môžem odstrániť všetky objekty.

rm(
  list = ls()
)

ls()

character(0)

Workspace je teraz prázdny.

3.4 Opätovné načítanie workspace

Načítam späť uložené objekty.

load(
  "pamat_Cvicenie5.RData"
)

ls()

[1] "cesta"          "data_scientist" "udaje"

Takýmto spôsobom môžem prenášať svoje pracovné prostredie medzi jednotlivými cvičeniami.

4 Lineárna regresia

Pokračujem s regresným modelom, v ktorom vysvetľujem plat Data Scientist pomocou:

roku,
úrovne skúseností,
veľkosti firmy.

Model má tvar:

\[ Salary_i = \beta_0 + \beta_1 Rok_i + \beta_2 Experience_i + \beta_3 CompanySize_i + u_i \]

Model odhadnem pomocou funkcie lm().

model <- lm(
  Salary.in.USD ~
    Rok.od.2020 +
    Experience.Level +
    Company.Size,
  data = data_scientist
)

summary(model)


Call:
lm(formula = Salary.in.USD ~ Rok.od.2020 + Experience.Level + 
    Company.Size, data = data_scientist)

Residuals:
    Min      1Q  Median      3Q     Max 
-135659  -35648   -3527   30791  290271 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  32197      11879   2.710 0.006904 ** 
Rok.od.2020                  11301       3344   3.379 0.000771 ***
Experience.LevelMid          16348       8152   2.005 0.045359 *  
Experience.LevelSenior       77673       7751  10.022  < 2e-16 ***
Experience.LevelExecutive    94384      14391   6.559 1.13e-10 ***
Company.SizeMedium           30753      11129   2.763 0.005889 ** 
Company.SizeLarge            11859      11557   1.026 0.305211    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 53430 on 628 degrees of freedom
Multiple R-squared:  0.3627,    Adjusted R-squared:  0.3566 
F-statistic: 59.57 on 6 and 628 DF,  p-value: < 2.2e-16

5 Regresné koeficienty a t-testy

Výsledky jednotlivých regresných koeficientov získam zo summary(model).

coef_mat <- summary(model)$coefficients

coef_mat

                          Estimate Std. Error   t value     Pr(>|t|)
(Intercept)               32196.96  11879.013  2.710407 6.903661e-03
Rok.od.2020               11301.33   3344.190  3.379392 7.714495e-04
Experience.LevelMid       16348.02   8152.397  2.005302 4.535937e-02
Experience.LevelSenior    77672.94   7750.556 10.021596 4.988176e-22
Experience.LevelExecutive 94383.80  14390.562  6.558729 1.134515e-10
Company.SizeMedium        30753.02  11128.988  2.763326 5.889235e-03
Company.SizeLarge         11859.18  11556.846  1.026160 3.052112e-01

Vytvorím prehľadnú tabuľku.

coef_tab <- data.frame(
  Premenna = rownames(coef_mat),
  Estimate = coef_mat[, 1],
  Std.Error = coef_mat[, 2],
  t.value = coef_mat[, 3],
  p.value = coef_mat[, 4],
  row.names = NULL
)

Pridám označenie štatistickej významnosti.

coef_tab$Significance <- ifelse(
  coef_tab$p.value < 0.001,
  "***",
  ifelse(
    coef_tab$p.value < 0.01,
    "**",
    ifelse(
      coef_tab$p.value < 0.05,
      "*",
      ifelse(
        coef_tab$p.value < 0.10,
        ".",
        ""
      )
    )
  )
)

knitr::kable(
  coef_tab,
  digits = 4,
  row.names = FALSE,
  col.names = c(
    "Premenná",
    "Odhad",
    "Štandardná chyba",
    "t hodnota",
    "p-hodnota",
    "Významnosť"
  ),
  caption = "Regresné koeficienty a výsledky t-testov"
)

Regresné koeficienty a výsledky t-testov
Premenná	Odhad	Štandardná chyba	t hodnota	p-hodnota	Významnosť
(Intercept)	32196.96	11879.013	2.7104	0.0069	**
Rok.od.2020	11301.33	3344.190	3.3794	0.0008	***
Experience.LevelMid	16348.02	8152.397	2.0053	0.0454	*
Experience.LevelSenior	77672.94	7750.556	10.0216	0.0000	***
Experience.LevelExecutive	94383.80	14390.562	6.5587	0.0000	***
Company.SizeMedium	30753.02	11128.988	2.7633	0.0059	**
Company.SizeLarge	11859.18	11556.846	1.0262	0.3052

6 Interpretácia t-testu

Pri každom regresnom koeficiente testujeme:

\[ H_0:\beta_j=0 \]

oproti:

\[ H_1:\beta_j\neq0. \]

Testovacia štatistika je:

\[ t_j = \frac{\hat{\beta}_j} {SE(\hat{\beta}_j)}. \]

Ak je p-hodnota menšia ako 0,05, nulovú hypotézu na 5 % hladine významnosti zamietame.

To znamená, že daný regresný koeficient považujeme za štatisticky významný.

6.1 Interpretácia roku

Koeficient pri premennej Rok.od.2020 je:

coef(model)["Rok.od.2020"]

Rok.od.2020 
   11301.33

Jeho p-hodnota je:

summary(model)$coefficients[
  "Rok.od.2020",
  4
]

[1] 0.0007714495

Pri rovnakých skúsenostiach a rovnakej veľkosti firmy sa zvýšenie roku o jednu jednotku spája so zmenou predikovaného platu približne o:

1.130133^{4} USD.

Koeficient pri roku je štatisticky významný na 5 % hladine významnosti.

6.2 Interpretácia Senior úrovne

coef(model)["Experience.LevelSenior"]

Experience.LevelSenior 
              77672.94

Senior Data Scientist má pri rovnakom roku a rovnakej veľkosti firmy predikovaný plat približne o

7.767294^{4} USD

vyšší oproti základnej kategórii Entry.

P-hodnota:

summary(model)$coefficients[
  "Experience.LevelSenior",
  4
]

[1] 4.988176e-22

Tento rozdiel je štatisticky významný.

6.3 Interpretácia Mid úrovne

Mid Data Scientist má oproti Entry Data Scientist pri ostatných nezmenených podmienkach rozdiel v predikovanom plate:

coef(model)["Experience.LevelMid"]

Experience.LevelMid 
           16348.02

P-hodnota:

summary(model)$coefficients[
  "Experience.LevelMid",
  4
]

[1] 0.04535937

Koeficient je štatisticky významný na 5 % hladine významnosti.

6.4 Interpretácia Executive úrovne

coef(model)["Experience.LevelExecutive"]

Experience.LevelExecutive 
                  94383.8

Executive Data Scientist má oproti Entry úrovni pri ostatných nezmenených premenných predikovaný plat vyšší približne o

9.43838^{4} USD.

Koeficient je štatisticky významný.

7 Označenie štatistickej významnosti

Pri výsledkoch používam tradičné označenie:

Symbol	p-hodnota	Interpretácia
`***`	p < 0.001	veľmi silná štatistická významnosť
`**`	p < 0.01	významné na 1 % hladine
`*`	p < 0.05	významné na 5 % hladine
`.`	p < 0.10	významné na 10 % hladine
bez symbolu	p >= 0.10	štatisticky nevýznamné

8 Kvalita vyrovnania modelu

Kvalitu regresného modelu posúdim pomocou:

koeficientu determinácie (R^2),
upraveného (R^2),
reziduálnej štandardnej chyby,
F-testu,
AIC a BIC.

Najskôr získam potrebné hodnoty.

model_summary <- summary(model)

f_stat <- model_summary$fstatistic

F_hodnota <- unname(
  f_stat[1]
)

df1 <- unname(
  f_stat[2]
)

df2 <- unname(
  f_stat[3]
)

F_p_value <- pf(
  F_hodnota,
  df1,
  df2,
  lower.tail = FALSE
)

Vytvorím tabuľku.

suhrn_modelu <- data.frame(

  R.squared =
    model_summary$r.squared,

  Adj.R.squared =
    model_summary$adj.r.squared,

  Sigma =
    model_summary$sigma,

  F.stat =
    F_hodnota,

  F.p.value =
    F_p_value,

  df =
    df1,

  df.residual =
    df.residual(model),

  N =
    nobs(model),

  AIC =
    AIC(model),

  BIC =
    BIC(model),

  Deviance =
    deviance(model)
)

knitr::kable(
  suhrn_modelu,
  digits = 4,
  row.names = FALSE,
  caption = "Súhrnné charakteristiky regresného modelu"
)

Súhrnné charakteristiky regresného modelu
R.squared	Adj.R.squared	Sigma	F.stat	F.p.value	df	df.residual	N	AIC	BIC	Deviance
0.3627	0.3566	53427.47	59.5651	0	6	628	635	15636.34	15671.96	1.792623e+12

9 Koeficient determinácie R²

Koeficient determinácie má tvar:

\[ R^2 = 1 - \frac{RSS}{TSS}. \]

Kde:

\[ RSS = \sum_{i=1}^{n} (y_i-\hat y_i)^2 \]

\[ TSS = \sum_{i=1}^{n} (y_i-\bar y)^2. \]

V mojom modeli je:

model_summary$r.squared

[1] 0.3626892

Teda model vysvetľuje približne:

round(
  model_summary$r.squared * 100,
  2
)

[1] 36.27

percent variability platov Data Scientist v použitej vzorke.

Hodnota (R^2) je približne 0.363, teda model zachytáva približne 36.3 % variability platu.

To neznamená, že rovnaké percento predikcií je správnych. (R^2) meria vysvetlenú variabilitu vysvetľovanej premennej.

10 Manuálne overenie R²

Najskôr vypočítam reziduálny súčet štvorcov.

RSS <- sum(
  residuals(model)^2
)

RSS

[1] 1.792623e+12

Celkový súčet štvorcov:

TSS <- sum(
  (
    data_scientist$Salary.in.USD -
      mean(data_scientist$Salary.in.USD)
  )^2
)

TSS

[1] 2.812792e+12

Manuálne vypočítané (R^2):

R2_manual <- 1 -
  RSS / TSS

R2_manual

[1] 0.3626892

Porovnanie:

c(
  R2_z_modelu =
    model_summary$r.squared,

  R2_manual =
    R2_manual
)

R2_z_modelu   R2_manual 
  0.3626892   0.3626892

Výsledky by mali byť rovnaké.

11 Upravený koeficient determinácie

Upravený koeficient determinácie zohľadňuje aj počet vysvetľujúcich premenných.

model_summary$adj.r.squared

[1] 0.3566002

Na rozdiel od obyčajného (R^2) nemusí po pridaní ďalšej premennej vždy rásť.

Je preto vhodnejší pri porovnávaní modelov s rôznym počtom vysvetľujúcich premenných.

12 Reziduálna štandardná chyba

sigma(model)

[1] 53427.47

Reziduálna štandardná chyba vyjadruje typickú veľkosť odchýlky medzi skutočným a modelom predikovaným platom.

Má rovnakú jednotku ako vysvetľovaná premenná, teda v našom prípade USD.

Overím jej výpočet.

sqrt(
  deviance(model) /
    df.residual(model)
)

[1] 53427.47

Táto hodnota by mala byť rovnaká ako výsledok sigma(model).

13 F-test celého modelu

F-test skúma spoločnú štatistickú významnosť vysvetľujúcich premenných.

Testujeme:

\[ H_0: \beta_1 = \beta_2 = \dots = \beta_k = 0 \]

oproti:

\[ H_1: \text{aspoň jeden smernicový koeficient nie je nulový}. \]

F-štatistika:

F_hodnota

[1] 59.56508

P-hodnota:

F_p_value

[1] 2.438836e-58

Keďže p-hodnota je 2.44e-58, na 5 % hladine významnosti

zamietame nulovú hypotézu.

To znamená, že

regresný model je ako celok štatisticky významný.

Je dôležité rozlišovať:

t-test hodnotí jednotlivý regresný koeficient,
F-test hodnotí model ako celok.

14 Grafická analýza kvality modelu

Pri diagnostike lineárneho regresného modelu môžem použiť príkaz:

plot(model)

Ten vytvorí štyri základné diagnostické grafy.

par(
  mfrow = c(
    2,
    2
  )
)

plot(model)

Štyri základné diagnostické grafy regresného modelu

par(
  mfrow = c(
    1,
    1
  )
)

15 Interpretácia diagnostických grafov

15.1 Residuals vs Fitted

Prvý graf porovnáva rezíduá s predikovanými hodnotami.

V ideálnom prípade by mali byť rezíduá náhodne rozptýlené okolo nuly a nemali by vytvárať výrazný systematický tvar.

Červená vyhladená čiara by mala byť približne horizontálna.

Výrazné zakrivenie by mohlo naznačovať, že lineárny model nezachytáva celý tvar vzťahu medzi premennými.

Ak by sa rozptyl rezíduí smerom doprava výrazne zväčšoval alebo zmenšoval, mohlo by ísť o heteroskedasticitu.

15.2 Normal Q-Q

Druhý graf slúži na vizuálne posúdenie normality rezíduí.

Ak sú rezíduá približne normálne rozdelené, body by mali ležať približne na diagonálnej čiare.

Výrazné odchýlky najmä na koncoch grafu naznačujú problém s normalitou alebo existenciu extrémnych hodnôt.

Normalitu následne overím aj formálnym Jarqueovým-Berovým testom.

15.3 Scale-Location

Tretí graf skúma stabilitu rozptylu rezíduí.

Ideálne by mali byť body rozptýlené približne rovnomerne a červená čiara by mala byť približne horizontálna.

Výrazný rast alebo pokles môže naznačovať heteroskedasticitu.

15.4 Residuals vs Leverage

Posledný graf spája veľkosť rezíduí a leverage jednotlivých pozorovaní.

Leverage ukazuje, ako neobvyklé je pozorovanie vzhľadom na hodnoty vysvetľujúcich premenných.

Cookova vzdialenosť zároveň pomáha identifikovať pozorovania, ktoré môžu výraznejšie ovplyvňovať výsledné regresné koeficienty.

Vysoká leverage sama osebe neznamená, že je dané pozorovanie chybné.

16 Test normality rezíduí

Na formálne posúdenie normality použijem Jarqueov-Berov test.

Testujeme:

\[ H_0: \text{rezíduá majú normálne rozdelenie} \]

oproti:

\[ H_1: \text{rezíduá nemajú normálne rozdelenie}. \]

Jarqueov-Berov test využíva šikmosť a špicatosť rezíduí.

Testovacia štatistika má tvar:

\[ JB = \frac{n}{6} \left( S^2 + \frac{(K-3)^2}{4} \right). \]

16.1 Výpočet Jarqueovho-Berovho testu

Najskôr získam rezíduá.

res <- residuals(model)

n_res <- length(res)

Vypočítam centrované rezíduá.

res_c <- res -
  mean(res)

Druhý, tretí a štvrtý moment:

m2 <- mean(
  res_c^2
)

m3 <- mean(
  res_c^3
)

m4 <- mean(
  res_c^4
)

Šikmosť:

S <- m3 /
  (
    m2^(3 / 2)
  )

S

[1] 0.6184228

Špicatosť:

K <- m4 /
  (
    m2^2
  )

K

[1] 4.504448

Jarqueova-Berova štatistika:

JB <- n_res / 6 *
  (
    S^2 +
      (
        (K - 3)^2 / 4
      )
  )

JB

[1] 100.3605

P-hodnota:

JB_p <- pchisq(
  JB,
  df = 2,
  lower.tail = FALSE
)

JB_p

[1] 1.61066e-22

Výsledok zobrazím v tabuľke.

JB_tabulka <- data.frame(
  JB_statistika = JB,
  p_hodnota = JB_p,
  Skewness = S,
  Kurtosis = K
)

knitr::kable(
  JB_tabulka,
  digits = 4,
  row.names = FALSE,
  caption = "Jarqueov-Berov test normality rezíduí"
)

Jarqueov-Berov test normality rezíduí
JB_statistika	p_hodnota	Skewness	Kurtosis
100.3605	0	0.6184	4.5044

17 Interpretácia Jarqueovho-Berovho testu

Rozhodovacie pravidlo je:

ak (p < 0.05), zamietame normalitu rezíduí,
ak (p ), normalitu nezamietame.

V mojom modeli je p-hodnota:

1.61e-22.

Preto na 5 % hladine významnosti

zamietam nulovú hypotézu o normálnom rozdelení rezíduí.

Formálny test je vhodné posudzovať spolu s Q-Q grafom.

18 Test odľahlých pozorovaní

Odľahlé pozorovania môžem skúmať pomocou externe studentizovaných rezíduí.

Použijem:

rstud <- rstudent(model)

Externe studentizované rezíduum má tvar:

\[ t_i = \frac{e_i} {\hat{\sigma}_{(i)} \sqrt{1-h_{ii}}}. \]

Pre každé pozorovanie vypočítam obojstrannú p-hodnotu.

df_outlier <-
  df.residual(model) - 1

p_unadjusted <- 2 *
  pt(
    abs(rstud),
    df = df_outlier,
    lower.tail = FALSE
  )

Keďže naraz testujem veľa pozorovaní, použijem Bonferroniho korekciu.

p_bonferroni <- pmin(
  1,
  nobs(model) *
    p_unadjusted
)

Vytvorím tabuľku.

model_data <- model.frame(model)

outlier_table <- data.frame(

  Pozorovanie =
    rownames(model_data),

  Plat_USD =
    model_data$Salary.in.USD,

  Studentizovane_reziduum =
    rstud,

  p_hodnota =
    p_unadjusted,

  Bonferroni_p =
    p_bonferroni,

  row.names = NULL
)

Zoradím ju podľa Bonferroniho p-hodnoty.

outlier_table <- outlier_table[
  order(
    outlier_table$Bonferroni_p
  ),
]

knitr::kable(
  head(
    outlier_table,
    10
  ),
  digits = 5,
  row.names = FALSE,
  caption = "Najvýraznejšie pozorovania podľa testu odľahlosti"
)

Najvýraznejšie pozorovania podľa testu odľahlosti
Pozorovanie	Plat_USD	Studentizovane_reziduum	p_hodnota	Bonferroni_p
3296	412000	5.64894	0.00000	0.00002
1752	370000	3.70111	0.00023	0.14830
2393	350000	3.53566	0.00044	0.27733
458	359170	3.49194	0.00051	0.32590
5	140000	-0.64693	0.51792	1.00000
6	126000	-0.90954	0.36342	1.00000
7	170000	-0.08479	0.93245	1.00000
8	130000	-0.83448	0.40433	1.00000
17	185900	0.21303	0.83137	1.00000
18	129300	-0.84761	0.39698	1.00000

19 Identifikácia štatisticky významných odľahlých pozorovaní

Za významne odľahlé budem považovať pozorovania, pri ktorých:

\[ p_{Bonferroni} < 0.05. \]

odlahle <- outlier_table[
  outlier_table$Bonferroni_p < 0.05,
]

Počet takýchto pozorovaní:

nrow(
  odlahle
)

[1] 1

Ak sa takéto pozorovania nachádzajú, zobrazím ich.

if (
  nrow(odlahle) == 0
) {

  cat(
    "Model neobsahuje štatisticky významné odľahlé pozorovania podľa Bonferroniho korekcie."
  )

} else {

  knitr::kable(
    odlahle,
    digits = 5,
    row.names = FALSE,
    caption = "Štatisticky významné odľahlé pozorovania"
  )
}

Štatisticky významné odľahlé pozorovania
Pozorovanie	Plat_USD	Studentizovane_reziduum	p_hodnota	Bonferroni_p
3296	412000	5.64894	0	2e-05

Odľahlé pozorovanie však nemusí byť automaticky chybou v databáze.

Môže predstavovať skutočného Data Scientist s nezvyčajne vysokým alebo nízkym platom.

Preto ho nemožno automaticky odstrániť iba preto, že ho test označil za odľahlé.

20 Leverage a Cookova vzdialenosť

Pri diagnostike vplyvných pozorovaní je dôležitá aj leverage.

Leverage získam pomocou:

leverage <- hatvalues(model)

Cookovu vzdialenosť získam:

cook <- cooks.distance(model)

Počet parametrov modelu:

p_model <- length(
  coef(model)
)

p_model

[1] 7

Priemerná leverage je:

p_model /
  nobs(model)

[1] 0.01102362

Ako orientačnú hranicu použijem:

\[ \frac{2p}{n}. \]

hranica_leverage <-
  2 *
  p_model /
  nobs(model)

hranica_leverage

[1] 0.02204724

Pre Cookovu vzdialenosť môžem ako orientačnú hodnotu použiť:

\[ \frac{4}{n}. \]

hranica_cook <-
  4 /
  nobs(model)

hranica_cook

[1] 0.006299213

20.1 Tabuľka diagnostických hodnôt

diagnostika <- data.frame(

  Pozorovanie =
    rownames(
      model.frame(model)
    ),

  Leverage =
    leverage,

  Cookova_vzdialenost =
    cook,

  Standardizovane_reziduum =
    rstandard(model),

  row.names = NULL
)

diagnostika <- diagnostika[
  order(
    diagnostika$Cookova_vzdialenost,
    decreasing = TRUE
  ),
]

Zobrazím desať najvýraznejších pozorovaní.

knitr::kable(
  head(
    diagnostika,
    10
  ),
  digits = 4,
  row.names = FALSE,
  caption = "Pozorovania s najvyššou Cookovou vzdialenosťou"
)

Pozorovania s najvyššou Cookovou vzdialenosťou
Pozorovanie	Leverage	Cookova_vzdialenost	Standardizovane_reziduum
3296	0.0295	0.1319	5.5148
383	0.0527	0.0377	-2.1777
1941	0.0527	0.0348	2.0916
190	0.0527	0.0345	2.0820
1510	0.0527	0.0245	-1.7546
800	0.0473	0.0242	1.8452
828	0.0925	0.0203	-1.1802
3298	0.0494	0.0145	1.3976
939	0.0527	0.0134	1.2983
2516	0.0399	0.0130	-1.4798

21 Graf leverage a Cookovej vzdialenosti

par(
  mfrow = c(
    1,
    2
  )
)

plot(
  leverage,
  type = "h",
  lwd = 2,
  xlab = "Pozorovanie",
  ylab = "Leverage",
  main = "Leverage"
)

abline(
  h = hranica_leverage,
  lty = 2
)

plot(
  cook,
  type = "h",
  lwd = 2,
  xlab = "Pozorovanie",
  ylab = "Cookova vzdialenosť",
  main = "Cookova vzdialenosť"
)

abline(
  h = hranica_cook,
  lty = 2
)

par(
  mfrow = c(
    1,
    1
  )
)

Leverage ukazuje nezvyčajnú polohu pozorovania vzhľadom na vysvetľujúce premenné.

Cookova vzdialenosť kombinuje informáciu o leverage a veľkosti rezídua.

Vysoká leverage preto nie je automaticky to isté ako vplyvné pozorovanie.

22 Úprava pôvodného modelu

Diagnostika môže ukázať, že jednoduchý lineárny vzťah medzi časom a platom nemusí úplne zachytiť vývoj platov.

V mojej databáze mám iba štyri roky – 2020 až 2023.

Preto vyskúšam model, ktorý okrem lineárneho trendu obsahuje aj kvadratický člen roku.

Nový model má tvar:

\[ Salary_i = \beta_0 + \beta_1 Rok_i + \beta_2 Rok_i^2 + \beta_3 Experience_i + \beta_4 CompanySize_i + u_i. \]

model2 <- lm(
  Salary.in.USD ~
    Rok.od.2020 +
    I(Rok.od.2020^2) +
    Experience.Level +
    Company.Size,
  data = data_scientist
)

summary(model2)


Call:
lm(formula = Salary.in.USD ~ Rok.od.2020 + I(Rok.od.2020^2) + 
    Experience.Level + Company.Size, data = data_scientist)

Residuals:
    Min      1Q  Median      3Q     Max 
-138013  -35873   -3706   30553  269280 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  46814      13399   3.494  0.00051 ***
Rok.od.2020                 -15314      11907  -1.286  0.19886    
I(Rok.od.2020^2)              6766       2906   2.328  0.02021 *  
Experience.LevelMid          15966       8126   1.965  0.04986 *  
Experience.LevelSenior       77409       7724  10.022  < 2e-16 ***
Experience.LevelExecutive    92891      14354   6.471 1.96e-10 ***
Company.SizeMedium           37706      11485   3.283  0.00108 ** 
Company.SizeLarge            18497      11864   1.559  0.11947    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 53240 on 627 degrees of freedom
Multiple R-squared:  0.3682,    Adjusted R-squared:  0.3611 
F-statistic: 52.19 on 7 and 627 DF,  p-value: < 2.2e-16

23 Koeficienty upraveného modelu

coef_mat2 <-
  summary(model2)$coefficients

coef_tab2 <- data.frame(
  Premenna =
    rownames(coef_mat2),

  Estimate =
    coef_mat2[, 1],

  Std.Error =
    coef_mat2[, 2],

  t.value =
    coef_mat2[, 3],

  p.value =
    coef_mat2[, 4],

  row.names = NULL
)

knitr::kable(
  coef_tab2,
  digits = 4,
  row.names = FALSE,
  col.names = c(
    "Premenná",
    "Odhad",
    "Štandardná chyba",
    "t hodnota",
    "p-hodnota"
  ),
  caption = "Koeficienty upraveného regresného modelu"
)

Koeficienty upraveného regresného modelu
Premenná	Odhad	Štandardná chyba	t hodnota	p-hodnota
(Intercept)	46814.35	13399.132	3.4938	0.0005
Rok.od.2020	-15314.05	11906.726	-1.2862	0.1989
I(Rok.od.2020^2)	6765.95	2905.864	2.3284	0.0202
Experience.LevelMid	15966.27	8125.504	1.9650	0.0499
Experience.LevelSenior	77408.69	7724.249	10.0215	0.0000
Experience.LevelExecutive	92891.22	14354.491	6.4712	0.0000
Company.SizeMedium	37706.18	11485.045	3.2831	0.0011
Company.SizeLarge	18497.36	11864.025	1.5591	0.1195

24 Porovnanie modelov pomocou AIC

Pri porovnávaní modelov s rovnakou vysvetľovanou premennou uprednostňujeme nižšiu hodnotu AIC.

porovnanie_modelov <- data.frame(

  Model = c(
    "Pôvodný model",
    "Model s kvadratickým trendom"
  ),

  AIC = c(
    AIC(model),
    AIC(model2)
  ),

  BIC = c(
    BIC(model),
    BIC(model2)
  ),

  R2 = c(
    summary(model)$r.squared,
    summary(model2)$r.squared
  ),

  Upraveny_R2 = c(
    summary(model)$adj.r.squared,
    summary(model2)$adj.r.squared
  )
)

knitr::kable(
  porovnanie_modelov,
  digits = 4,
  row.names = FALSE,
  caption = "Porovnanie pôvodného a upraveného modelu"
)

Porovnanie pôvodného a upraveného modelu
Model	AIC	BIC	R2	Upraveny_R2
Pôvodný model	15636.34	15671.96	0.3627	0.3566
Model s kvadratickým trendom	15632.87	15672.95	0.3682	0.3611

Na základe AIC je vhodnejší:

model s kvadratickým trendom, pretože má nižšiu hodnotu AIC.

Rozdiel AIC je:

abs(
  AIC(model) -
    AIC(model2)
)

[1] 3.466917

25 Grafická diagnostika upraveného modelu

par(
  mfrow = c(
    2,
    2
  )
)

plot(model2)

Diagnostické grafy upraveného regresného modelu

par(
  mfrow = c(
    1,
    1
  )
)

26 Normalita rezíduí upraveného modelu

Vypočítam Jarqueov-Berov test aj pre druhý model.

res2 <- residuals(model2)

n_res2 <- length(res2)

res2_c <- res2 -
  mean(res2)

m2_2 <- mean(
  res2_c^2
)

m3_2 <- mean(
  res2_c^3
)

m4_2 <- mean(
  res2_c^4
)

S2 <- m3_2 /
  (
    m2_2^(3 / 2)
  )

K2 <- m4_2 /
  (
    m2_2^2
  )

JB2 <- n_res2 / 6 *
  (
    S2^2 +
      (
        (K2 - 3)^2 / 4
      )
  )

JB2_p <- pchisq(
  JB2,
  df = 2,
  lower.tail = FALSE
)

JB2_tabulka <- data.frame(
  JB_statistika = JB2,
  p_hodnota = JB2_p,
  Skewness = S2,
  Kurtosis = K2
)

knitr::kable(
  JB2_tabulka,
  digits = 4,
  row.names = FALSE,
  caption = "Jarqueov-Berov test upraveného modelu"
)

Jarqueov-Berov test upraveného modelu
JB_statistika	p_hodnota	Skewness	Kurtosis
71.191	0	0.566	4.1872

Pri upravenom modeli je p-hodnota:

3.48e-16.

Preto:

normalitu rezíduí ani pri upravenom modeli na 5 % hladine významnosti neprijímame.

27 Odľahlé pozorovania upraveného modelu

rstud2 <- rstudent(model2)

df_outlier2 <-
  df.residual(model2) - 1

p_unadjusted2 <- 2 *
  pt(
    abs(rstud2),
    df = df_outlier2,
    lower.tail = FALSE
  )

p_bonferroni2 <- pmin(
  1,
  nobs(model2) *
    p_unadjusted2
)

model_data2 <- model.frame(model2)

outlier_table2 <- data.frame(

  Pozorovanie =
    rownames(model_data2),

  Plat_USD =
    model_data2$Salary.in.USD,

  Studentizovane_reziduum =
    rstud2,

  p_hodnota =
    p_unadjusted2,

  Bonferroni_p =
    p_bonferroni2,

  row.names = NULL
)

outlier_table2 <- outlier_table2[
  order(
    outlier_table2$Bonferroni_p
  ),
]

knitr::kable(
  head(
    outlier_table2,
    10
  ),
  digits = 5,
  row.names = FALSE,
  caption = "Odľahlé pozorovania upraveného modelu"
)

Odľahlé pozorovania upraveného modelu
Pozorovanie	Plat_USD	Studentizovane_reziduum	p_hodnota	Bonferroni_p
3296	412000	5.32401	0.00000	0.00009
1752	370000	3.66944	0.00026	0.16751
2393	350000	3.64554	0.00029	0.18352
458	359170	3.45961	0.00058	0.36685
5	140000	-0.69361	0.48818	1.00000
6	126000	-0.95724	0.33881	1.00000
7	170000	-0.12936	0.89712	1.00000
8	130000	-0.88189	0.37817	1.00000
17	185900	0.16957	0.86541	1.00000
18	129300	-0.89508	0.37109	1.00000

28 Porovnanie normality oboch modelov

normalita_porovnanie <- data.frame(

  Model = c(
    "Pôvodný model",
    "Kvadratický model"
  ),

  JB = c(
    JB,
    JB2
  ),

  p_hodnota = c(
    JB_p,
    JB2_p
  )
)

knitr::kable(
  normalita_porovnanie,
  digits = 5,
  row.names = FALSE,
  caption = "Porovnanie normality rezíduí oboch modelov"
)

Porovnanie normality rezíduí oboch modelov
Model	JB	p_hodnota
Pôvodný model	100.36046	0
Kvadratický model	71.19103	0

Nižšia Jarqueova-Berova štatistika znamená, že rozdelenie rezíduí je z pohľadu šikmosti a špicatosti bližšie k normálnemu rozdeleniu.

Ani nižšia hodnota JB však automaticky neznamená splnenie normality. Rozhodnutie robím podľa p-hodnoty a zároveň pomocou Q-Q grafu.

29 Záverečné porovnanie modelov

Pôvodný model:

\[ Salary_i = \beta_0 + \beta_1 Rok_i + \beta_2 Experience_i + \beta_3 CompanySize_i + u_i \]

Upravený model pridáva aj kvadratický časový trend:

\[ Salary_i = \beta_0 + \beta_1 Rok_i + \beta_2 Rok_i^2 + \beta_3 Experience_i + \beta_4 CompanySize_i + u_i. \]

Pri rozhodovaní medzi modelmi sledujem:

AIC,
upravené (R^2),
diagnostické grafy,
normalitu rezíduí,
odľahlé a vplyvné pozorovania.

Na základe AIC vychádza ako lepší

upravený model s kvadratickým trendom.

Samotné AIC však nestačí. Model musí byť zároveň ekonomicky a vecne interpretovateľný a jeho diagnostické vlastnosti by mali byť prijateľné.

30 Záver

V tomto cvičení som pokračoval v analýze databázy Latest Data Science Salaries a pracovnej pozície Data Scientist.

Najskôr som si ukázal, ako uložiť celý workspace pomocou save.image() a ako ho neskôr obnoviť pomocou load().

Následne som odhadol viacnásobný lineárny regresný model platu Data Scientist.

Pri jednotlivých regresných koeficientoch som interpretoval Studentove t-testy a p-hodnoty.

Kvalitu modelu ako celku som posudzoval pomocou koeficientu determinácie (R^2), upraveného (R^2), reziduálnej štandardnej chyby a F-testu.

Ďalej som pomocou plot(model) vykonal grafickú diagnostiku rezíduí.

Normalitu rezíduí som formálne overil Jarqueovým-Berovým testom.

Pomocou externe studentizovaných rezíduí a Bonferroniho korekcie som identifikoval potenciálne odľahlé pozorovania.

S využitím leverage a Cookovej vzdialenosti som posúdil aj možné vplyvné pozorovania.

Nakoniec som vytvoril druhý regresný model s kvadratickým časovým trendom a oba modely som porovnal pomocou AIC, upraveného (R^2) a diagnostiky rezíduí.

Toto cvičenie ukazuje, že pri dátovej analýze nestačí regresný model iba odhadnúť. Dôležité je tiež overiť jeho štatistickú významnosť, kvalitu vyrovnania a základné diagnostické vlastnosti.

31 Uloženie workspace na ďalšie cvičenie

Na konci si uložím celý aktuálny workspace.

save.image(
  "pamat_Cvicenie5.RData"
)

Tento súbor môžem v ďalšom cvičení načítať pomocou:

# load("pamat_Cvicenie5.RData")