library(ggplot2)
library(patchwork)
library(lmtest)
library(sandwich)
library(dplyr)
# Načítanie dát
df <- read.csv("C:/Users/Lenovo/Documents/School/Mgr/R - Ekonometria/Cvicenia/My dataset/ai_job_dataset.csv")

# Rýchly pohľad
glimpse(df)
## Rows: 15,000
## Columns: 19
## $ job_id                 <chr> "AI00001", "AI00002", "AI00003", "AI00004", "AI…
## $ job_title              <chr> "AI Research Scientist", "AI Software Engineer"…
## $ salary_usd             <int> 90376, 61895, 152626, 80215, 54624, 123574, 796…
## $ salary_currency        <chr> "USD", "USD", "USD", "USD", "EUR", "EUR", "GBP"…
## $ experience_level       <chr> "SE", "EN", "MI", "SE", "EN", "SE", "MI", "EN",…
## $ employment_type        <chr> "CT", "CT", "FL", "FL", "PT", "CT", "FL", "FL",…
## $ company_location       <chr> "China", "Canada", "Switzerland", "India", "Fra…
## $ company_size           <chr> "M", "M", "L", "M", "S", "M", "S", "L", "L", "M…
## $ employee_residence     <chr> "China", "Ireland", "South Korea", "India", "Si…
## $ remote_ratio           <int> 50, 100, 0, 50, 100, 50, 0, 0, 0, 0, 100, 0, 10…
## $ required_skills        <chr> "Tableau, PyTorch, Kubernetes, Linux, NLP", "De…
## $ education_required     <chr> "Bachelor", "Master", "Associate", "PhD", "Mast…
## $ years_experience       <int> 9, 1, 2, 7, 0, 7, 3, 0, 7, 5, 8, 15, 5, 0, 6, 0…
## $ industry               <chr> "Automotive", "Media", "Education", "Consulting…
## $ posting_date           <chr> "2024-10-18", "2024-11-20", "2025-03-18", "2024…
## $ application_deadline   <chr> "2024-11-07", "2025-01-11", "2025-04-07", "2025…
## $ job_description_length <int> 1076, 1268, 1974, 1345, 1989, 819, 1936, 1286, …
## $ benefits_score         <dbl> 5.9, 5.2, 9.4, 8.6, 6.6, 5.9, 6.3, 7.6, 9.3, 5.…
## $ company_name           <chr> "Smart Analytics", "TechCorp Inc", "Autonomous …
# Jednoduchý model - mzda podľa skúseností, remote podielu a dĺžky popisu práce
model <- lm(salary_usd ~ years_experience + remote_ratio + job_description_length, data = df)

summary(model)
## 
## Call:
## lm(formula = salary_usd ~ years_experience + remote_ratio + job_description_length, 
##     data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -128744  -24882   -5978   18816  254186 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            65617.9816  1078.0473  60.867   <2e-16 ***
## years_experience        8013.6687    59.9279 133.722   <2e-16 ***
## remote_ratio               3.6081     8.1431   0.443    0.658    
## job_description_length    -0.3716     0.5768  -0.644    0.519    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 40700 on 14996 degrees of freedom
## Multiple R-squared:  0.544,  Adjusted R-squared:  0.5439 
## F-statistic:  5963 on 3 and 14996 DF,  p-value: < 2.2e-16

Vizuálne testovanie heteroskedasticity

resids_sq <- resid(model)^2

p1 <- ggplot(df, aes(x = years_experience, y = resids_sq)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", color = "red", se = FALSE) +
  labs(x = "Years of Experience", y = "Squared Residuals",
       title = "Squared Residuals vs Experience") +
  theme_minimal()

p2 <- ggplot(df, aes(x = remote_ratio, y = resids_sq)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", color = "red", se = FALSE) +
  labs(x = "Remote Ratio", y = "Squared Residuals",
       title = "Squared Residuals vs Remote Ratio") +
  theme_minimal()

p1 + p2
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 50
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 50
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0

Interpretácia:

Graf 1 – Squared Residuals vs Experience Na osi X sú roky skúseností, na osi Y štvorce rezíduí (teda veľkosť chýb modelu). Vidno, že s rastúcimi skúsenosťami sa rozptyl rezíduí mierne zvyšuje – červená LOESS krivka stúpa. To znamená, že pri vyššom počte skúseností model robí väčšie chyby – variabilita miezd je väčšia u seniorných pracovníkov. → Záver: slabá až stredná heteroskedasticita spôsobená premennou years_experience.

Graf 2 – Squared Residuals vs Remote Ratio Tu je na osi X percento práce na diaľku. Červená krivka má výrazný vlnovitý (oscilujúci) tvar – pri 0 %, 50 % a 100 % sa rozptyl správa odlišne. To ukazuje, že model nie je konzistentný naprieč rôznymi typmi remote práce – niektoré skupiny (napr. plne remote vs. onsite) majú výrazne iný rozptyl chýb. → Záver: veľmi jasná heteroskedasticita – rozptyl rezíduí sa výrazne mení v závislosti od remote_ratio.

Breusch–Pagan test

bptest(model)
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 1958.7, df = 3, p-value < 2.2e-16

Interpretácia výsledku:

p-hodnota je extrémne malá (< 0.001) → nulová hypotéza o konštantnom rozptyle rezíduí sa jednoznačne zamieta. To znamená, že v modeli je výrazná heteroskedasticita – rozptyl chýb nie je konštantný. Štatistika BP = 1958,7 pri troch stupňoch voľnosti (tri vysvetľujúce premenné) je veľmi vysoká → rozdiely v rozptyle sú silné a systematické. → Záver: heteroskedasticita je jasne prítomná — reziduá sa menia s hodnotami vysvetľujúcich premenných

Náprava pomocou White robustných odhadov

Ak test indikuje heteroskedasticitu, upravíme štandardné chyby koeficientov:

coeftest(model, vcov = vcovHC(model, type = "HC1"))
## 
## t test of coefficients:
## 
##                           Estimate  Std. Error  t value Pr(>|t|)    
## (Intercept)            65617.98160  1003.62230  65.3812   <2e-16 ***
## years_experience        8013.66866    75.69810 105.8635   <2e-16 ***
## remote_ratio               3.60807     8.08247   0.4464   0.6553    
## job_description_length    -0.37160     0.57041  -0.6515   0.5148    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretácia:

Konštanta aj skúsenosti ostávajú vysoko štatisticky významné – aj po korekcii. → To znamená, že vzťah medzi skúsenosťami a mzdou je robustný a spoľahlivý.

Remote ratio a job_description_length nemajú štatisticky významný vplyv (p > 0.05). → Po korekcii sa ukazuje, že rozdiely v platoch medzi remote a on-site pozíciami nie sú štatisticky významné, keď zohľadníme heteroskedasticitu.

White korekcia neodstraňuje heteroskedasticitu, len opravuje štandardné chyby tak, aby t-testy boli korektné.

Záver

Ak Breusch–Pagan test ukáže nízke p-value a vizuálne grafy potvrdzujú trend vo variabilite rezíduí, znamená to, že v modeli existuje heteroskedasticita – teda rozptyl chýb závisí od vysvetľujúcich premenných.
Po aplikácii Whiteovej korekcie môžeme pokračovať v interpretácii výsledkov, pretože koeficienty zostávajú nestranné a testy už zohľadňujú nekonštantný rozptyl.