Špecifikácia modelu

Cieľom analýzy je modelovať výsledné skúškové skóre študentov (exam_score) na základe:

hours_studied – počet hodín štúdia, sleep_hours – priemerný počet hodín spánku, attendance_percent – percentuálna účasť na vyučovaní.

#######################################################################
# PRIPRAVA UDAJOV
#######################################################################
df <- read_csv("student_exam_scores.csv") %>%
select(exam_score, hours_studied, sleep_hours, attendance_percent)

# Imputácia mediánom

column_medians <- sapply(df, median, na.rm = TRUE)
for(col in names(df)){
df[[col]][is.na(df[[col]])] <- column_medians[col]
}

summary(df)

   exam_score    hours_studied     sleep_hours    attendance_percent
 Min.   :17.10   Min.   : 1.000   Min.   :4.000   Min.   : 50.30    
 1st Qu.:29.50   1st Qu.: 3.500   1st Qu.:5.300   1st Qu.: 62.20    
 Median :34.05   Median : 6.150   Median :6.700   Median : 75.25    
 Mean   :33.95   Mean   : 6.325   Mean   :6.622   Mean   : 74.83    
 3rd Qu.:38.75   3rd Qu.: 9.000   3rd Qu.:8.025   3rd Qu.: 87.42    
 Max.   :51.30   Max.   :12.000   Max.   :9.000   Max.   :100.00    

par(mfrow=c(2,2))
cols <- c("#d62828","#006d2c","#f4d35e","#ffffff") # červená, zelená, zlatá, biela
i <- 1
for(col in names(df)){
boxplot(df[[col]], main=col, col=cols[i %% length(cols) + 1], border=cols[(i+1) %% length(cols) + 1])
i <- i+1
}
par(mfrow=c(1,1))

Diagnostické grafy základného modelu nám poskytujú cenné informácie o vhodnosti špecifikácie regresie a o vlastnostiach rezíduí.

Graf rezíduá versus vyrovnané hodnoty ukazuje, že rezíduá sú rozptýlené náhodne okolo nulovej línie, bez výrazného zakrivenia alebo systematického vzoru, čo naznačuje, že lineárna špecifikácia modelu je primeraná a nie je potrebné transformovať premenné.

Q–Q graf, ktorý porovnáva kvantily rezíduí s kvantilmi normálneho rozdelenia, ukazuje, že body sa pohybujú veľmi blízko diagonálnej línie, len s menšími odchýlkami na okrajoch, čo potvrdzuje, že rezíduá sú približne normálne rozložené a predpoklad normality nie je porušený.

Scale–Location graf, ktorý znázorňuje štandardizované rezíduá voči vyrovnaným hodnotám, ukazuje rovnomerné rozptýlenie bodov bez zjavného rastúceho alebo klesajúceho vzoru, čo naznačuje, že predpoklad homoskedasticity je splnený.

Nakoniec, diagnostika pomocou Cookovej vzdialenosti a identifikácie odľahlých hodnôt neodhalila žiadne extrémne pozorovania, ktoré by významne deformovali model, takže všetky pozorovania sú spoľahlivé a model je stabilný. Celkovo tieto grafy podporujú tvrdenie, že lineárny model pre vysvetľovanie skúškového skóre študentov je adekvátny, rezíduá majú požadované vlastnosti a výsledky regresie sú interpretovateľné.

1. Test RESET (test chyby špecifikácie Ramseyho regresnej rovnice - Ramsey Reset Test)

Testovanie, či lineárny model pre exam_score je správne špecifikovaný, vykonáme pridaním mocnín vyrovnaných hodnôt (fitted values) do modelu. Pôvodný model:

exam_scoret​=β0​+β1​⋅hours_studiedt​+β2​⋅sleep_hourst​+β3​⋅attendance_percentt​+ut​

a rozšírený model pre test:

exam_scoret​=β0​+β1​⋅hours_studiedt​+β2​⋅sleep_hourst​+β3​⋅attendance_percentt​+γ2​y^{​t2​+γ3​y}​t3​+ut​

kde y^​t sú vyrovnané hodnoty z pôvodného modelu. Testuje sa hypotéza:

H0: model je správne špecifikovaný (y2 je rovné y3 a to je rovné 0) H1: model je nesprávne špecifikovaný (y2 nie je rovné 0 alebo y3 nie je rovné 0)

# základný model
model <- lm(exam_score ~ hours_studied + sleep_hours + attendance_percent, data=df)

# RESET test
library(lmtest)
resettest(model)

    RESET test

data:  model
RESET = 0.56353, df1 = 2, df2 = 194, p-value = 0.5701

Interpretácia:

Ak je p-hodnota testu menšia ako 0.05, zamietame nulovú hypotézu a model je zrejme nesprávne špecifikovaný. To by znamenalo, že lineárna forma nemusí byť úplne vhodná, môže chýbať niektorá premenná alebo je vhodná nelineárna transformácia (napr. logaritmus hodín štúdia alebo kvadrát premenných). Ak je p-hodnota väčšia ako 0.05, nulovú hypotézu neprijímame a lineárny model je adekvátny pre naše dáta.

2. Grafická analýza

Graf Residuals vs. Fitted

Najprv vykreslíme klasický graf rezíduá vs. vyrovnané hodnoty pre základný model. Tento graf nám ukáže, či rezíduá vykazujú náhodný vzor. Ak je rezíduí príliš zakrivený alebo systematický, môže to naznačovať potrebu nelineárnej transformácie premenných.

# Residuals vs Fitted

plot(model, which = 1, col="#d62828", pch=19, cex=0.6)

Ak rezíduá vykazujú náhodný rozptyl okolo nuly a nie sú žiadne jasné zakrivenia, model je lineárny a vhodne špecifikovaný. Ak by sa objavila zakrivená štruktúra, zvážime transformáciu niektorých premenných, napríklad logaritmus alebo druhú mocninu.

Grafy C+R **

# C+R Plots pre všetky vysvetľujúce premenné

car::crPlots(model, col="#006d2c", pch=19)

Ak krivka pre danú premennú sleduje približne priamku, lineárny predpoklad je vhodný.

Ak krivka vykazuje zakrivenie, odporúča sa transformácia danej premennej, napríklad logaritmus hours_studied alebo kvadratický člen I(attendance_percent^2).

V prípade dát môžeme očakávať, že hours_studied môže vykazovať miernu nelinearitu (napr. úbytok prínosu pri vysokom počte hodín štúdia), zatiaľ čo sleep_hours a attendance_percent môžu byť lineárne.

3. Porovnanie základného a modifikovaného modelu

Predpokladajme, že C+R grafy naznačujú, že vzťah medzi výsledným skóre študentov a premennými hours_studied a attendance_percent môže byť nelineárny. Zavedieme preto do modelu kvadratické členy pre tieto premenné. Pôvodný lineárny model:

model <- lm(exam_score ~ hours_studied + sleep_hours + attendance_percent, data = df)
model_poly <- lm(exam_score ~ hours_studied + I(hours_studied^2) +
                               sleep_hours +
                               attendance_percent + I(attendance_percent^2),
                 data = df)
summary(model_poly)

Call:
lm(formula = exam_score ~ hours_studied + I(hours_studied^2) + 
    sleep_hours + attendance_percent + I(attendance_percent^2), 
    data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.8885 -2.5420 -0.1928  2.9303  7.9036 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)   
(Intercept)              2.095799   8.359165   0.251  0.80230   
hours_studied            1.196918   0.399812   2.994  0.00311 **
I(hours_studied^2)       0.033364   0.029985   1.113  0.26721   
sleep_hours              0.576782   0.182831   3.155  0.00186 **
attendance_percent       0.395018   0.225507   1.752  0.08141 . 
I(attendance_percent^2) -0.001857   0.001512  -1.228  0.22086   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.848 on 194 degrees of freedom
Multiple R-squared:  0.6868,    Adjusted R-squared:  0.6787 
F-statistic: 85.08 on 5 and 194 DF,  p-value: < 2.2e-16

anova(model, model_poly)

Analysis of Variance Table

Model 1: exam_score ~ hours_studied + sleep_hours + attendance_percent
Model 2: exam_score ~ hours_studied + I(hours_studied^2) + sleep_hours + 
    attendance_percent + I(attendance_percent^2)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    196 2915.8                           
2    194 2873.1  2    42.717 1.4422 0.2389

library(lmtest)
resettest(model_poly)

    RESET test

data:  model_poly
RESET = 0.065584, df1 = 2, df2 = 192, p-value = 0.9365

V modifikovanom modeli môžu byť kvadratické členy štatisticky významné, čo naznačuje, že vzťah nie je priamo lineárny. Ak upravený koeficient determinácie R2adj stúpol, znamená to, že model lepšie vysvetľuje variabilitu skúškových výsledkov. Ak niektorý kvadratický člen (napr. pre attendance_percent) nie je štatisticky významný, môžeme ho vylúčiť a zachovať len významné kvadratické členy, čím získame jednoduchejší a interpretovateľný model.

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

Predpokladajme, že chceme analyzovať, či sa vzťah medzi počtom hodín štúdia (hours_studied) a skúškovým skóre mení nad určitým prahom. Vytvoríme dummy premennú DUM, ktorá nad určitým počtom hodín nadobúda hodnotu 1 a inak 0. Tým môžeme testovať:

1. Zlom v autonómnom člene – posun regresnej priamky pre študentov, ktorí študujú nad daný prah:

df$DUM <- ifelse(df$hours_studied > 10, 1, 0)  # prah 10 hodín, upravi podľa potreby
modelD_auto <- lm(exam_score ~ DUM + hours_studied + sleep_hours + attendance_percent, data=df)
summary(modelD_auto)

Call:
lm(formula = exam_score ~ DUM + hours_studied + sleep_hours + 
    attendance_percent, data = df)

Residuals:
   Min     1Q Median     3Q    Max 
-8.594 -2.584 -0.150  3.167  7.444 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        11.06387    1.98162   5.583 7.84e-08 ***
DUM                 0.57099    0.97955   0.583  0.56063    
hours_studied       1.58173    0.11836  13.364  < 2e-16 ***
sleep_hours         0.58107    0.18352   3.166  0.00179 ** 
attendance_percent  0.11937    0.01924   6.205 3.21e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.864 on 195 degrees of freedom
Multiple R-squared:  0.6827,    Adjusted R-squared:  0.6762 
F-statistic: 104.9 on 4 and 195 DF,  p-value: < 2.2e-16

2. Zlom v sklone – efekt hours_studied môže byť iný pre študentov nad a pod prahom:

modelD_sklon <- lm(exam_score ~ hours_studied + I(DUM*hours_studied) +
                                     sleep_hours + attendance_percent,
                   data=df)
summary(modelD_sklon)

Call:
lm(formula = exam_score ~ hours_studied + I(DUM * hours_studied) + 
    sleep_hours + attendance_percent, data = df)

Residuals:
   Min     1Q Median     3Q    Max 
-8.602 -2.557 -0.121  3.149  7.353 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)            11.09197    1.98099   5.599 7.24e-08 ***
hours_studied           1.57297    0.11888  13.232  < 2e-16 ***
I(DUM * hours_studied)  0.06085    0.08909   0.683  0.49540    
sleep_hours             0.58164    0.18346   3.170  0.00177 ** 
attendance_percent      0.11944    0.01923   6.211 3.12e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.862 on 195 degrees of freedom
Multiple R-squared:  0.6829,    Adjusted R-squared:  0.6764 
F-statistic:   105 on 4 and 195 DF,  p-value: < 2.2e-16

Týmto modelom umožníme, aby sa sklon vzťahu medzi hours_studied a exam_score líšil pre študentov, ktorí študujú menej alebo viac ako 10 hodín.

Porovnanie pôvodného lineárneho modelu a modelu so zlomom vykonáme ANOVA testom a overíme pomocou RESET testu:

anova(model, modelD_sklon)

Analysis of Variance Table

Model 1: exam_score ~ hours_studied + sleep_hours + attendance_percent
Model 2: exam_score ~ hours_studied + I(DUM * hours_studied) + sleep_hours + 
    attendance_percent
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    196 2915.8                           
2    195 2908.9  1    6.9593 0.4665 0.4954

library(lmtest)
resettest(modelD_sklon)

    RESET test

data:  modelD_sklon
RESET = 0.2949, df1 = 2, df2 = 193, p-value = 0.7449

Ak sa upravený koeficient determinácie R2adj zvýši, model lepšie vysvetľuje variabilitu výsledkov. Koeficient pri interakcii DUM*hours_studied nám ukazuje, ako sa efekt hodín štúdia mení pre študentov nad prahom. Ak je interakcia významná, vzťah medzi štúdiom a skóre nie je konštantný a zlom je opodstatnený. Ak nie, lineárny efekt hours_studied je postačujúci.