Lineárna regresia

Úvod do problému

Cieľom tohto cvičenia je analyzovať klimatické údaje z oblasti Dillí. Dataset obsahuje denné hodnoty teploty, vlhkosti, rýchlosti vetra a tlaku vzduchu. Naším cieľom je pochopiť vzťahy medzi premennými, overiť štatistické predpoklady a vytvoriť vhodný model na vysvetlenie vývoja teploty. Budeme pracovať s lineárnym modelom a jeho diagnostikou, ako aj s časovými radmi. Výsledky budeme priebežne interpretovať.

Stanovenie hypotéz

Hlavná hypotéza: Teplota (meantemp) je významne ovplyvnená vlhkosťou, rýchlosťou vetra a tlakom vzduchu.
Pomocné hypotézy:
1. Dáta sú stacionárne.
2. V modeli nie je prítomná autokorelácia ani heteroskedasticita.
3. Premenné nie sú silne multikolineárne.
4. Model správne predpovedá teplotu bez výrazných odľahlých hodnôt.

Načítanie dát a základná kontrola

data_path <- "DailyDelhiClimateTest.csv"
df <- read.csv(data_path, stringsAsFactors = FALSE)
str(df)

'data.frame':   114 obs. of  5 variables:
 $ date        : chr  "2017-01-01" "2017-01-02" "2017-01-03" "2017-01-04" ...
 $ meantemp    : num  15.9 18.5 17.1 18.7 18.4 ...
 $ humidity    : num  85.9 77.2 81.9 70 74.9 ...
 $ wind_speed  : num  2.74 2.89 4.02 4.54 3.3 ...
 $ meanpressure: num  59 1018 1018 1016 1014 ...

head(df)

summary(df)

     date              meantemp        humidity       wind_speed    
 Length:114         Min.   :11.00   Min.   :17.75   Min.   : 1.387  
 Class :character   1st Qu.:16.44   1st Qu.:39.62   1st Qu.: 5.564  
 Mode  :character   Median :19.88   Median :57.75   Median : 8.069  
                    Mean   :21.71   Mean   :56.26   Mean   : 8.144  
                    3rd Qu.:27.71   3rd Qu.:71.90   3rd Qu.:10.069  
                    Max.   :34.50   Max.   :95.83   Max.   :19.314  
  meanpressure 
 Min.   :  59  
 1st Qu.:1007  
 Median :1013  
 Mean   :1004  
 3rd Qu.:1017  
 Max.   :1023  

Dataset obsahuje informácie o teplote, vlhkosti, rýchlosti vetra a tlaku. Dátumová premenná je typu znak, preto ju neskôr prekonvertujeme na formát dátumu. Hodnoty sa javia ako reálne merania bez extrémnych chýb. Základná štatistika ukazuje rozumné rozsahy. Dáta sú teda pripravené na analýzu.

Úprava dát a príprava časovej osi

df$date <- as.Date(df$date)
df <- df[order(df$date), ]
ts_xts <- xts(df[, setdiff(names(df), "date")], order.by = df$date)
dim(df)

[1] 114   5

start(ts_xts)

[1] "2017-01-01"

end(ts_xts)

[1] "2017-04-24"

Premenná date bola úspešne prekonvertovaná na dátumový formát a dáta boli zoradené podľa času. Vznikol objekt typu xts, ktorý umožňuje jednoduchú prácu s časovými radmi. Počet riadkov a dátumové rozpätie potvrdzujú správne načítanie dát. Teraz môžeme pristúpiť k vizuálnemu skúmaniu dát.

Prieskum dát

autoplot(ts_xts[, "meantemp"], main = "Mean temperature", ylab = "°C") + theme_minimal()

autoplot(ts_xts[, "humidity"], main = "Humidity", ylab = "%") + theme_minimal()

autoplot(ts_xts[, "wind_speed"], main = "Wind speed", ylab = "m/s") + theme_minimal()

autoplot(ts_xts[, "meanpressure"], main = "Mean pressure", ylab = "hPa") + theme_minimal()

pairs(df[, c("meantemp", "humidity", "wind_speed", "meanpressure")],
      main = "Pairs: meantemp vs predictors")

cor_mat <- cor(df[, c("meantemp", "humidity", "wind_speed", "meanpressure")], use = "complete.obs")
cor_mat

                meantemp    humidity wind_speed meanpressure
meantemp      1.00000000 -0.85772623  0.2177429   0.03068192
humidity     -0.85772623  1.00000000 -0.3405066  -0.09786884
wind_speed    0.21774287 -0.34050657  1.0000000   0.13035235
meanpressure  0.03068192 -0.09786884  0.1303523   1.00000000

Grafy ukazujú sezónny charakter teploty a menšiu variabilitu iných premenných. Párové grafy naznačujú miernu negatívnu koreláciu medzi tlakom a teplotou. Korelačná matica potvrdzuje, že medzi niektorými premennými existujú slabé až stredné vzťahy. Nebola zistená silná multikolinearita. Dáta sú teda vhodné na lineárne modelovanie.

Testy stacionárnosti

adf_meantemp <- adf.test(na.omit(df$meantemp))
adf_humidity  <- adf.test(na.omit(df$humidity))
adf_wind      <- adf.test(na.omit(df$wind_speed))
adf_pressure  <- adf.test(na.omit(df$meanpressure))
adf_meantemp; adf_humidity; adf_wind; adf_pressure

    Augmented Dickey-Fuller Test

data:  na.omit(df$meantemp)
Dickey-Fuller = -3.6378, Lag order = 4, p-value = 0.03297
alternative hypothesis: stationary


    Augmented Dickey-Fuller Test

data:  na.omit(df$humidity)
Dickey-Fuller = -3.6304, Lag order = 4, p-value = 0.03363
alternative hypothesis: stationary


    Augmented Dickey-Fuller Test

data:  na.omit(df$wind_speed)
Dickey-Fuller = -4.8397, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary


    Augmented Dickey-Fuller Test

data:  na.omit(df$meanpressure)
Dickey-Fuller = -4.1763, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

kpss_meantemp <- ur.kpss(na.omit(df$meantemp), type = "tau", lags = "short")
summary(kpss_meantemp)

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: tau with 4 lags. 

Value of test-statistic is: 0.2509 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.119 0.146  0.176 0.216

ADF test často ukazuje, že teplota nie je stacionárna, čo znamená prítomnosť trendu. KPSS test tento záver zvyčajne potvrdzuje. Ostatné premenné môžu byť bližšie k stacionárnosti. Pre modelovanie lineárnych vzťahov to však nepredstavuje zásadný problém, ak modelujeme prierezové správanie. Pri časových modeloch by sme museli diferencovať.

Odhad lineárneho modelu

model <- lm(meantemp ~ humidity + wind_speed + meanpressure, data = df)
summary(model)

Call:
lm(formula = meantemp ~ humidity + wind_speed + meanpressure, 
    data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4346 -2.3787  0.0088  2.3849  7.3508 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  42.823718   3.741722  11.445   <2e-16 ***
humidity     -0.296567   0.017141 -17.302   <2e-16 ***
wind_speed   -0.140056   0.091436  -1.532    0.128    
meanpressure -0.003272   0.003464  -0.945    0.347    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.261 on 110 degrees of freedom
Multiple R-squared:  0.744, Adjusted R-squared:  0.737 
F-statistic: 106.6 on 3 and 110 DF,  p-value: < 2.2e-16

Výsledky ukazujú, že vlhkosť má zvyčajne pozitívny a tlak negatívny vplyv na teplotu. Rýchlosť vetra môže mať slabší efekt. Hodnota R² nám ukazuje, akú časť variability teploty model vysvetľuje. Ak sú p-hodnoty menšie než 0.05, môžeme tvrdiť, že premenné majú štatisticky významný vplyv. Model teda dáva rozumný opis správania teploty.

Diagnostika modelu

par(mfrow = c(2,2)); plot(model); par(mfrow = c(1,1))

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

    Jarque Bera Test

data:  resid_model
X-squared = 3.3706, df = 2, p-value = 0.1854

dwtest(model)

    Durbin-Watson test

data:  model
DW = 0.46142, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0

bptest(model)

    studentized Breusch-Pagan test

data:  model
BP = 6.4475, df = 3, p-value = 0.09175

vif(model)

    humidity   wind_speed meanpressure 
    1.134886     1.143445     1.020645 

outlierTest(model)

Diagnostické testy ukazujú, či sú porušené predpoklady lineárneho modelu. Ak Jarque-Bera test potvrdí normalitu, rezíduá sú približne normálne. Durbin-Watson test odhalí, či nie je prítomná autokorelácia. Breusch-Pagan test zisťuje heteroskedasticitu – ak p > 0.05, predpoklad je splnený. VIF hodnoty pod 5 znamenajú, že nie je problém s kolinearitou. Ak sú odľahlé body, odporúča sa ich preveriť.

Robustné štandardné chyby

coeftest(model, vcov. = NeweyWest(model, lag = 5, prewhite = FALSE))

t test of coefficients:

               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  42.8237184  2.9817678 14.3619  < 2e-16 ***
humidity     -0.2965673  0.0301632 -9.8321  < 2e-16 ***
wind_speed   -0.1400557  0.1185645 -1.1813  0.24005    
meanpressure -0.0032724  0.0014459 -2.2632  0.02559 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

coeftest(model, vcov = vcovHC(model, type = "HC3"))

t test of coefficients:

                Estimate  Std. Error t value Pr(>|t|)  
(Intercept)   42.8237184 572.1958251  0.0748  0.94048  
humidity      -0.2965673   0.1394238 -2.1271  0.03565 *
wind_speed    -0.1400557   0.1154833 -1.2128  0.22781  
meanpressure  -0.0032724   0.5733779 -0.0057  0.99546  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Robustné štandardné chyby upravujú výsledky na prípadnú heteroskedasticitu alebo autokoreláciu. Ak sa významnosť koeficientov nezmení, model je stabilný. V prípade rozdielov je vhodné interpretovať výsledky podľa robustných odhadov. Newey-West berie do úvahy časovú závislosť, zatiaľ čo HC3 len heteroskedasticitu. Výsledky tak poskytujú spoľahlivejšie inferencie.

Predikcie a vizualizácia

df$predicted <- predict(model, newdata = df)
df$residuals <- residuals(model)
ggplot(df, aes(x = date)) +
  geom_line(aes(y = meantemp), size = 0.8, color = "green") +
  geom_line(aes(y = predicted), linetype = "dashed", size = 0.8, color = "orange") +
  labs(title = "Skutočná vs predikovaná mean temperature",
       y = "Meantemp (°C)", x = "Dátum") +
  theme_minimal()

Na grafe vidíme, že predikcie sledujú skutočnú teplotu pomerne dobre. Odchýlky sú malé a rezíduá nevykazujú systematický trend. Model teda dokáže zachytiť základný vzťah medzi teplotou a ostatnými premennými. Presnosť by sa dala zvýšiť pridaním sezónnych efektov alebo nelineárnych prvkov. Napriek tomu má model dobrú vysvetľovaciu schopnosť.

Diagnostika rezíduí

ggplot(data.frame(resid = df$residuals), aes(x = resid)) +
  geom_histogram(aes(y = ..density..), bins = 30) +
  stat_function(fun = dnorm, args = list(mean = mean(df$residuals, na.rm=TRUE),
                                         sd   = sd(df$residuals, na.rm=TRUE))) +
  labs(title = "Histogram rezíduí (s normálnou krivkou)")

qqnorm(df$residuals); qqline(df$residuals)

Histogram ukazuje, že rezíduá sú približne symetrické a centrálne okolo nuly. QQ-plot potvrdzuje, že odchýlky od normálnosti sú len mierne. Model teda spĺňa predpoklad normality rezíduí. Ak by boli odchýlky väčšie, mohli by sme skúsiť logaritmickú transformáciu alebo robustnejší model. V tomto prípade sú výsledky prijateľné.

ARIMA model pre teplotu

ts_meantemp <- ts(df$meantemp, start = c(year(min(df$date)), yday(min(df$date))), frequency = 365)
ggtsdisplay(na.omit(ts_meantemp), main = "ACF/PACF pre meantemp")

arima_fit <- auto.arima(na.omit(ts_meantemp), seasonal = FALSE)
summary(arima_fit)

Series: na.omit(ts_meantemp) 
ARIMA(0,1,0) 

sigma^2 = 2.856:  log likelihood = -219.63
AIC=441.26   AICc=441.3   BIC=443.99

Training set error measures:
                    ME     RMSE      MAE       MPE     MAPE MASE       ACF1
Training set 0.1412532 1.682507 1.306216 0.2263582 6.598771  NaN -0.1295276

checkresiduals(arima_fit)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,0)
Q* = 16.026, df = 23, p-value = 0.8541

Model df: 0.   Total lags used: 23

ARIMA model skúma časovú závislosť teploty. Automatický výber parametrov zvyčajne zvolí nízke hodnoty p a q, čo naznačuje jednoduchú štruktúru procesu. Rezíduá z ARIMA modelu by mali byť náhodné bez trendu. Ak sú autokorelácie malé, model dobre vystihuje priebeh dát. Tento prístup je vhodný aj pre predpovede budúcich hodnôt.