Έκδοση: 23 / 04 / 2024

Προετοιμασία περιβάλλοντος R
Οι παρακάτω εντολές, ελέγχουν αν έχουν εγκατασταθεί τα απαιτούμενα πακέτα για την εκτέλεση του συνόλου του κώδικα της ενότητας.

list.of.packages <- c("ggplot2", "plyr", "htmlTable", "ppcor", "psych", "car", "rms", "lmtest", "broom", "see", "performance", "forecast","tseries", "sarima", 'TSstudio')
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

Όλες οι συναρτήσεις που χρησιμοποιούνται στην ενότητα αυτή βρίσκονται στο αρχείο της γλώσσας R: MyRFunctions.R.

Για να γίνουν διαθέσιμες για χρήση αρκεί να φορτωθούν στο περιβάλλον της R, εκτελώντας την εντολή

source("https://utopia.duth.gr/epdiaman/files/kedivim/MyRFunctions.R")

(Στην περίπτωση όπου γίνει λήψη τοπικά το αρχείο MyRFunctions.R, πρέπει να αλλαξει ανάλογα και η διεύθυνση του αρχείου)

Η γλώσσα R περιέχει μία πληθώρα μεθόδων που μπορούν να χρησιμοποιηθούν για την ανάλυση χρονοσειρών. Μία καταγραφή τους γίνεται εδώ: https://cran.r-project.org/web/views/TimeSeries.html.

1 Ευθεία

Η ευθεία είναι το γεωμετρικό αντικείμενο που προκύπτει αν ενώσουμε δύο σημεία και προεκτείνουμε απεριόριστα και ως προς τις δύο κατευθύνσεις.

Σε περίπτωση όπου το επίπεδο εφοδιαστεί με ένα σύστημα αξόνων, η ευθεία αποκτά εξίσωση της μορφής

y = α x + β.

Το α ονομάζεται κλίση της ευθείας και δηλώνει τη μεταβολή στην τιμή του y όταν το x αυξηθεί κατά 1 μονάδα.

Το β ονομάζεται σταθερός όρος και αντιστοιχεί στην τιμή που παίρνει το y όταν το x είναι 0.

2 Απλή Γραμμική Παλινδρόμηση (Simple Linear Regression - SLR).

Η γραμμική παλινδρόμηση χρησιμοποιείται για την πρόβλεψη της τιμής μίας εξαρτημένης συνεχούς μεταβλητής από μία ή περισσότερες ανεξάρτητες μεταβλητές που ενδεχομένως το επηρεάζουν. Οι ανεξάρτητες μεταβλητές μπορούν να είναι συνεχείς είτε δίτιμες ποιοτικές μεταβλητές (0 - 1) ή ακόμα και ποιοτικές με περισσότερες τιμές αφού πρώτα αυτές αναπαρασταθούν από κατάλληλο πλήθος από δίτιμες ψευδομεταβλητές (dummy variables).

Στο dataframe my.df παρουσιάζονται οι βαθμοί στην πρόοδο (test) και στην τελική εξέταση (exam) από 14 φοιτητές.

test = c(5.5, 6, 7, 5, 8, 2, 9, 10, 2, 3, 4, 6.5, 8.5, 1)
exam = c(6.5, 6, 8, 7.5, 7, 4, 8, 10, 1, 5, 5, 6, 9, 5)
my.df = data.frame(test, exam)
Hmisc::label(my.df$test) = 'Πρόοδος'
Hmisc::label(my.df$exam) = 'Τελικές εξετάσεις'

Εκτελούμε την εντολή

sjPlot::plot_scatter(my.df, test, exam)

Είναι εμφανές πως τα σημεία βρίσκονται «κοντά» σε κάποια ευθεία, γεγονός που ερμηνεύεται ως υψηλή αναλογικότητα μεταξύ της test και της exam. Συνεπώς, έχει νόημα η αναζήτηση της με τη μέθοδο της απλής γραμμικής παλινδρόμησης.

Για τον υπολογισμό του γραμμικού μοντέλου αρκεί να εφαρμόσουμε τη συνάρτηση my.lm.

test = c(5.5, 6, 7, 5, 8, 2, 9, 10, 2, 3, 4, 6.5, 8.5, 1)
exam = c(6.5, 6, 8, 7.5, 7, 4, 8, 10, 1, 5, 5, 6, 9, 5)
my.df = data.frame(test, exam)
Hmisc::label(my.df$test) = 'Πρόοδος'
Hmisc::label(my.df$exam) = 'Τελικές εξετάσεις'


my.model = my_lm(my.df, exam ~ test)
my.model$htmlreport

					95% C.I.
Linear regression results for model exam ~ const + test
	B	SE	t	p	Lower	Upper
Constant	2.501	0.747	3.348	0.006	0.873	4.129
Πρόοδος	0.684	0.121	5.651	<0.001	0.42	0.947

Coefficient of determination
R² = 0.727

Homoscedasticity Breusch–Pagan test
x²(1) = 2.905, p = 0.088. Homoscedasticity assumption is not rejected. Model seems to be valid.

	Statistic	Description
Normality tests of models’ residuals
Skewness & Kurtosis
Skewness	-0.704	Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
Kurtosis	3.619	Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η₀: This sample is from a normal distribution vs Η₁: Not the Η₀.
Shapiro–Wilk	W = 0.943, p = 0.455	This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
Lilliefors	D = 0.136, p = 0.69	The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.

Οι συντελεστές της εξίσωσης είναι στην στήλη B. Παρατηρούμε ότι η ευθεία γραμμικής παλινδρόμησης έχει εξίσωση:

Εξετάσεις = 2,501 + 0,684 * Πρόοδος

Για την γραφική αναπαράσταση της ευθείας γραμμικής παλινδρόμησης μπορούμε να εφαρμόσουμε ξανά τη συνάρτηση plot_scatter:

sjPlot::plot_scatter(my.df, test, exam, fit.line = 'lm')

Αξιοποίηση της γραμμικής εξίσωσης πρόβλεψης
Η εξίσωση γραμμικής παλινδρόμησης αξιοποιείται με τρεις τρόπους.

Α) Ο πρώτος είναι για την πρόβλεψη των τιμών της εξαρτημένης από την/τις ανεξάρτητες. Για παράδειγμα, η πρόβλεψη για την τελική επίδοση στις εξετάσεις για έναν φοιτητή που είχε επίδοση 7 στην Πρόοδο είναι:

  Εξετάσεις = 2,501 + 0,684 * 7 = 7,289 ≈ 7,3

Β) Ο δεύτερος τρόπος είναι με την ανάγνωση των συντελεστών της εξίσωσης οι οποίοι αντιστοιχούν στην μεταβολή των τιμών της εξαρτημένης που αντιστοιχεί σε αύξηση μίας μονάδας στην ανεξάρτητη.

Στο παράδειγμά μας, προκύπτει πως για κάθε μία μονάδα μεγαλύτερη επίδοση στην Πρόοδο περιμένουμε 0,684 μονάδες μεγαλύτερη επίδοση στην τελική εξέταση.

Γ) Ο τρίτος τρόπος είναι η ερμηνεία του σταθερού όρου ο οποίος αντιστοιχεί στην τιμή της εξαρτημένης μεταβλητής όταν η τιμή της ανεξάρτητης είναι ίση με μηδέν.

Στο παράδειγμά μας, προκύπτει πως ο φοιτητής που θα γράψει 0 στην Πρόοδο περιμένουμε να γράψει 2,501 στις εξετάσεις.

Συντελεστής Προσδιορισμού
Η συνολική αποτελεσματικότητα του γραμμικού μοντέλου όπως αυτή εκφράζεται με το συντελεστή προσδιορισμού R² (Multiple R square) είναι ίση με 0,727 = 72,7% δηλαδή το 72,7% της συνολικής μεταβλητότητας της τελικής επίδοσης ενός φοιτητή ερμηνεύεται από την επίδοσή του στην Πρόοδο. Συμπεραίνουμε ανάλογα πως το υπόλοιπο 0,273 = 27,3% οφείλεται σε άλλους παράγοντες που δεν καταγράφονται στην έρευνά μας.

Στήλες t, p
Το στατιστικό t = B / Std.Error χρησιμοποιείται για να ελεγχθεί για κάθε μία από τις ανεξάρτητες μεταβλητές που συμμετέχουν στο γραμμικό μοντέλο, η στατιστική υπόθεση:

Η₀: Β = 0,

έναντι της:

Η₁: Β ≠ 0.

Εύλογο είναι πως προαπαιτείται η απόρριψη της στατιστικής υπόθεσης Η₀ για να αποδεχθούμε τη σημαντικότητα της αντίστοιχης ανεξάρτητης μεταβλητής ως σημαντικού παράγοντα επιρροής της εξαρτημένης μεταβλητής όπως επίσης και να προχωρήσουμε σε ερμηνεία του συντελεστή της.

Η στατιστική σημαντικότητα που αντιστοιχεί σε αυτόν τον (δίπλευρο) έλεγχο βρίσκεται στη στήλη p. Οι βαθμοί ελευθερίας για τον έλεγχο αυτό είναι df = n - 2 = 14 - 2 = 12.

Για τα δεδομένα του παραδείγματός μας παρατηρούμε πως t(12) =5, 651, p < 0,001, άρα η Πρόοδος είναι σημαντική μεταβλητή στο γραμμικό μοντέλο. Επιπλέον, παρατηρούμε πως η σταθερά (Constant) του μοντέλου είναι στατιστικά μη μηδενική ποσότητα στο γραμμικό μοντέλο που μελετούμε (t(12) = 3,348, p = 0,005).

Σημειώνεται, πως απόρριψη της σημαντικότητας της σταθεράς ενδέχεται να οδηγήσει σε επανάληψη της διαδικασίας επιλέγοντας υπολογισμό γραμμικού μοντέλου χωρίς σταθερά.

95% Διάστημα εμπιστοσύνης για τους συντελεστές
Το 95% διάστημα εμπιστοσύνης ερμηνεύεται ως μέτρο της ασφάλειας με την οποία μπορεί να καταγραφεί η επιρροή κάθε μίας από τις ανεξάρτητες μεταβλητές στην εξαρτημένη. Φανερά, αρμόζει να γίνει αναφορά μόνο στους συντελεστές που αντιστοιχούν σε μεταβλητές που υπολογίζονται να είναι στατιστικά σημαντικοί παράγοντες μεταβολής των τιμών της εξαρτημένης μεταβλητής.

Για τα δεδομένα του παραδείγματός μας καταγράφεται ότι αύξηση μίας μονάδας στην Πρόοδο αντιστοχεί σε μέση αύξηση από 0,42 έως 0,947 μονάδες την επίδοση στις τελικές εξετάσεις, ενώ επίδοση στην Πρόοδο ίση με 0 αντιστοχεί σε μέση επίδοση στις τελικές εξετάσεις από 0,873 έως 4,129 μονάδες.

Υπάρχουν δύο προϋποθέσεις εγκυρότητας ενός γραμμικού μοντέλου:
α) η κανονικότητα των υπολοίπων (residuals) του γραμμικού μοντέλου που θα προκύψει μετά το τέλος της διαδικασίας.
β) η ομοσκεδαστικότητα των παρατηρήσεων δηλαδή η ομοιογένεια της διακύμανσης των τιμών της εξαρτημένης μεταβλητής ανάμεσα στις τιμές των ανεξάρτητων μεταβλητών.

Οι παραπάνω προϋποθέσεις της κανονικότητας και της ομοσκεδαστικότητας ελέγχονται με παρατήρηση της κατανομής των υπολοίπων του μοντέλου, δηλαδή από τις διαφορές y_predicted - y_observed. Η συνάρτηση my.lm εφαρμόζει δύο καθιερωμένους στατιστικούς ελέγχους για την επαλήθευση αυτών των προϋποθέσεων.

Δοκιμασία ελέγχου κανονικότητας των υπολοίπων
Η υπόθεση πως τα υπόλοιπα (y_predicted - y_observed) έχουν κανονική κατανομή ελέγχεται με τη δοκιμασία Lilliefors η οποία είναι μία περισσότερο έγκυρη εκδοχή της δοκιμασίας Kolmogorov-Smirnov. Από το Output παρατηρούμε πως η υπόθεση αυτή δεν απορρίπτεται.

Δοκιμασία ελέγχου ομοσκεδαστικότητας
Η υπόθεση της ομοσκεδαστικότητας ελέγχεται με τη δοκιμασία Breusch-Pagan. Από το Output παρατηρούμε πως η υπόθεση της ομοσκεδαστικότητας δεν απορρίπτεται.

Αν ο ερευνητής επιθυμεί περαιτέρω ανάλυση, τότε μπορεί να ανακτήσει τα υπόλοιπα ως εξής:

fit = my.model$fit
round(fit$residuals, 3)

## Τελικές εξετάσεις 
##      1      2      3      4      5      6      7      8      9     10     11 
##  0.239 -0.603  0.713  1.581 -0.971  0.132 -0.654  0.662 -2.868  0.448 -0.236 
##     12     13     14 
## -0.945  0.688  1.815

Ένα επιπλέον χρήσιμο διάγραμμα αναπαράστασης είναι το εξής:

sjPlot::plot_residuals(fit)

Η αναπαράσταση της κατανομής τους μπορεί να γίνει με ένα ιστόγραμμα. Αν αυτό δεν μοιάζει με την κανονική κατανομή τότε υποδεικνύεται πως η γραμμική εξίσωση δεν είναι κατάλληλη για την επεξηγηση της σχέσης των δύο μεταβλητών.

my_histFrequency_with_normal_curve(fit$residuals)

Μία καλή αναπαράσταση των υπολοίπων και της επιρροής τους στην εγκυρότητα του γραμμικού μοντέλου είναι διαθέσιμη εδώ και εδώ.

2.1 Παρατηρήσεις

2.1.1 Τυποποιημένοι Συντελεστές (Standardized Coefficients)

Οι τυποποιημένοι συντελεστές για τις ανεξάρτητες μεταβλητές του μοντέλου θεωρούνται περισσότερο αξιόπιστες μαρτυρίες για την ακριβή σχέση της κάθε ανεξάρτητης μεταβλητή με την εξαρτημένη καθώς αναφέρουν την μεταβολή της εξαρτημένης μεταβλητής σε μονάδες τυπικής απόκλισης που επιφέρει αύξηση μίας τυπικής απόκλισης στην ανεξάρτητη μεταβλητή. Είναι εξαιρετικά χρήσιμες σε γραμμικά μοντέλα με πολλές ανεξάρτητες μεταβλητές όπου με απλή παρατήρηση των τιμών αυτών εύκολα ταξινομούνται οι ανεξάρτητες μεταβλητών ως προς την επιρροή τους στις τιμές της εξαρτημένης. Στην R μπορούν να υπολογιστούν με την εντολή

lm(scale(exam) ~ scale(test))

## 
## Call:
## lm(formula = scale(exam) ~ scale(test))
## 
## Coefficients:
## (Intercept)  scale(test)  
##   2.409e-17    8.526e-01

2.1.2 Γραμμικό μοντέλο χωρίς σταθερό όρο

Αν η εξαρτημένη μεταβλητή μπορεί από τη φύση της να είναι ακριβώς 0 όταν όλες οι ανεξάρτητες είναι 0 τότε ο ερευνητής έχει τη δυνατότητα να υπολογίσει γραμμικό μοντέλο που δεν θα περιέχει σταθερό όρο. Στην R αυτό μπορεί να γίνει με την εντολή:

lm(exam ~ test - 1)

## 
## Call:
## lm(formula = exam ~ test - 1)
## 
## Coefficients:
##  test  
## 1.047

2.1.3 Η προϋπόθεση της κανονικότητας

Η προϋπόθεση της κανονικότητας των υπολοίπων πολλές φορές συγχέεται με την κανονικότητα των κατανομών της εξαρτημένης ή των ανεξάρτητων μεταβλητών που εμπλέκονται σε ένα γραμμικό μοντέλο, μία συσχέτιση που είναι ανακριβής.
Τόσο η εξαρτημένη όσο και οι ανεξάρτητες μεταβλητές μπορούν να ακολουθούν μη κανονικές κατανομές αλλά τα υπόλοιπα του μοντέλου να είναι κανονικά καθιστώντας έγκυρη την αξιολόγηση του αντίστοιχου γραμμικού μοντέλου.

Παράδειγμα
Εξαρτημένη μεταβλητή Υ: δίτιμη (0 - 1) ποιοτική μεταβλητή, ανεξάρτητη μεταβλητή Χ = Υ + έναν τυχαίο αριθμό από τη Ν(0,1) κατανομή.

Φανερά, τόσο η Υ όσο και η Χ δεν ακολουθούν την κανονική κατανομή αλλά τα υπόλοιπα του γραμμικού μοντέλου Υ = α*Χ + β την ακολουθούν εκ κατασκευής. Η αναπαράσταση του παραπάνω παραδείγματος μπορεί να γίνει με τον εξής κώδικα R:

library(psych)
x = rbinom(100, 1, 0.3)
y = rnorm(length(x), 5 + x * 5, 1)
scatter.hist(x, y, correl=TRUE, density=TRUE, ellipse=FALSE, xlab="x", ylab="y")

fit = lm(y~x)
summary(fit)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.00664 -0.75043 -0.05146  0.58207  2.32124 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.7493     0.1128   42.09   <2e-16 ***
## x             5.5126     0.2027   27.20   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9374 on 98 degrees of freedom
## Multiple R-squared:  0.883,  Adjusted R-squared:  0.8818 
## F-statistic: 739.8 on 1 and 98 DF,  p-value: < 2.2e-16

2.1.4 Standardized και Studentized Residual Value

Ως υπόλοιπο (residual) μίας παρατήρησης ονομάζεται η διαφορά της πραγματικής τιμής της εξαρτημένης μεταβλητής από την προβλεπόμενη βάσει του γραμμικού μοντέλου ενώ το τυποποιημένο υπόλοιπο (standardized residual) είναι το αποτέλεσμα της διαίρεσης του υπολοίπου με την τυπική απόκλιση όλων των υπολοίπων.
Το τυποποιημένο κατά student υπόλοιπο (studentized residual) είναι το αποτέλεσμα της διαίρεσης του υπολοίπου με το τυπικό σφάλμα του.

Φανερά, και οι δύο αυτές τιμές είναι δυνατό να υπολογιστούν ύστερα από την ολοκλήρωση όλων των απαραίτητων υπολογισμών και τη δημιουργία του γραμμικού μοντέλου.

3 Απλή Γραμμική Παλινδρόμηση για εκτίμηση χρονοσειράς

H απλή γραμμική παλινδρόμηση μπορεί να χρησιμοποιηθεί για μοντελοποίηση χρονοσειρών ως προς το χρόνο. Ωστόσο, τα μοντέλα αυτής της ομάδας είναι από τη φύση τους απλά και χωρίς ιδιαίτερη αποτελεσματικότητα.

3.1 Παράδειγμα Ι

Στο dataframe a.df παρουσιάζεται η παραγωγή ζαχαρότευτλου (σε τόνους) σε μία περιοχή ανά έτος από το έτος 2010 (-3) έως το έτος 2016 (+3). Πρώτα προσαρμόζουμε γραμμικό μοντέλο πρόβλεψης της εξέλιξης της παραγωγής στο χρόνο.

years = c(-3, -2, -1, 0, 1, 2, 3)
production = c(40, 45, 46, 42, 47, 50, 46)
a.df = data.frame(years, production)

my.model = my_lm(a.df, production ~ years)
my.model$htmlreport

					95% C.I.
Linear regression results for model production ~ const + years
	B	SE	t	p	Lower	Upper
Constant	45.143	0.997	45.258	<0.001	42.579	47.707
years	1.036	0.499	2.077	0.092	-0.246	2.318

Coefficient of determination
R² = 0.463

Homoscedasticity Breusch–Pagan test
x²(1) = 0.288, p = 0.591. Homoscedasticity assumption is not rejected. Model seems to be valid.

	Statistic	Description
Normality tests of models’ residuals
Skewness & Kurtosis
Skewness	-0.19	Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
Kurtosis	1.315	Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η₀: This sample is from a normal distribution vs Η₁: Not the Η₀.
Shapiro–Wilk	W = 0.873, p = 0.199	This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
Lilliefors	D = 0.23, p = 0.321	The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.

sjPlot::plot_scatter(a.df, years, production, fit.line = 'lm')

3.2 Παράδειγμα ΙΙ

Στα πλαίσια της απλής γραμμικής παλινδρόμησης εντάσσεται και το AR(1) μοντέλο. Η προσαρμογή με τη συνάρτηση lm αποδίδει αντίστοιχο μοντέλο με το απότελεσμα της συνάρτησης Arima, με τις μικρές διαφορές στους υπολογισμούς να οφείλονται στο διαφορετικό αλγόριθμο (lm: ελάχιστα τετράγωνα vs Arima: MLE με φίλτρα Kalman).

Θα αξιοποιήσουμε τα δεδομένα uschange από τη βιβλιοθήκη fpp2.

library(fpp2)
Consumption = as.vector(uschange[,"Consumption"])
Consumption.lag1 = c(NA, head(Consumption, -1))
data.lm = data.frame(Consumption, Consumption.lag1)

3.2.1 Προσαρμογή ευθείας ελαχίστων τετραγώνων

Πρώτα, προσαρμόζουμε ένα μοντέλο απλής γραμμικής παλινδρόμησης με την εντολή my_lm που βασίζεται στη συνάρτηση lm της baseR.

my.model = my_lm(data.lm, Consumption ~ Consumption.lag1)
my.model$htmlreport

					95% C.I.
Linear regression results for model Consumption ~ const + Consumption.lag1
	B	SE	t	p	Lower	Upper
Constant	0.486	0.069	7.084	<0.001	0.351	0.621
Consumption.lag1	0.35	0.069	5.068	<0.001	0.214	0.486

Coefficient of determination
R² = 0.122

Homoscedasticity Breusch–Pagan test
x²(1) = 9.82, p = 0.002. Homoscedasticity assumption is rejected! Read model coefficients with caution.

	Statistic	Description
Normality tests of models’ residuals
Skewness & Kurtosis
Skewness	-0.481	Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
Kurtosis	5.106	Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η₀: This sample is from a normal distribution vs Η₁: Not the Η₀.
Shapiro–Wilk	W = 0.967, p = 0	This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
Lilliefors	D = 0.069, p = 0.031	The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.

sjPlot::plot_scatter(data.lm, Consumption.lag1, Consumption, fit.line = 'lm')

Παρατηρούμε ότι από τη δοκιμασία γραμμικής παλινδρόμησης προτείνεται η εξίσωση πρόβλεψης:

Consumption_t = 0.48593 + 0.34992 Consumption_t-1

Σημείωση
Στην baseR υπάρχει η συνάρτηση ar.ols η οποία προσαρμόζει μοντέλα της μορφής AR(p) σε μία χρονοσειρά με την μέθοδο των ελαχίστων τετραγώνων. Ενδεικτικά, η εφαρμοφή της στα ίδια δεδομένα θα έδινε το εξής αποτέλεσμα:

ar.ols(Consumption, order.max = 1, demean=F, intercept=T)

## 
## Call:
## ar.ols(x = Consumption, order.max = 1, demean = F, intercept = T)
## 
## Coefficients:
##      1  
## 0.3499  
## 
## Intercept: 0.4859 (0.06823) 
## 
## Order selected 1  sigma^2 estimated as  0.3768

Αυτή η συνάρτηση μπορεί να χρησιμοποιηθεί για την αυτόματη προσαρμογή του βέλτιστου AR μοντέλου σε μία χρονοσειρά, με την μέθοδο των ελαχίστων τετραγώνων, επιλέγοντας το βέλτιστο μοντέλο βάσει του δείκτη AIC. Ενδεικτικά, η εφαρμογή της για το διάνυσμα Consumption είναι η παρακάτω:

ar.ols(Consumption, demean=F, intercept=T)

## 
## Call:
## ar.ols(x = Consumption, demean = F, intercept = T)
## 
## Coefficients:
##       1        2        3        4        5        6        7        8  
##  0.2279   0.2072   0.1880  -0.0061   0.0040   0.0186  -0.0017  -0.0747  
##       9       10       11       12       13       14       15       16  
##  0.0065  -0.0044   0.0726  -0.0663  -0.0042  -0.1004   0.0645   0.0775  
##      17       18       19       20       21       22  
## -0.0500   0.0118  -0.0474  -0.0168  -0.0363   0.1634  
## 
## Intercept: 0.2677 (0.1322) 
## 
## Order selected 22  sigma^2 estimated as  0.2649

3.2.2 Προσαρμογή ευθείας με MLE

Στη συνέχεια, εφαρμόζουμε τη συνάρτηση Arima και προσαρμόζουμε ένα AR(1) μοντέλο πρόβλεψης.

Arima(Consumption, order = c(1, 0, 0))

## Series: Consumption 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1    mean
##       0.3481  0.7460
## s.e.  0.0683  0.0685
## 
## sigma^2 = 0.3789:  log likelihood = -173.67
## AIC=353.33   AICc=353.46   BIC=363.02

Βρίσκουμε, πως από τη συνάρτηση Arima, προτείνεται το μοντέλο

(Consumption_t - 0.746) = 0.3481(Consumption_t-1 - 0.746)

Εύκολα μπορούμε να επαληθεύσουμε ότι η δεύτερη εξίσωση είναι η ίδια με την πρώτη.

Σημείωση
Στη γλώσσα R χρειάζεται λίγη προσοχή για την ανάγνωση των μοντέλων στην περίπτωση όπου υπάρχει σταθερά. Ενδεικτικά, πληροφορίες υπάρχουν εδώ: https://otexts.com/fpp2/arima-r.html#understanding-constants-in-r ή και εδώ: https://www.stat.pitt.edu/stoffer/tsa2/Rissues.htm

4 Πολλαπλή Γραμμική Παλινδρόμηση (Multiple Linear Regression - MLR).

Η πολλαπλή γραμμική παλινδρόμηση είναι η γενίκευση της απλής γραμμικής παλινδρόμησης όταν για την πρόβλεψη της τιμής μίας εξαρτημένης συνεχούς μεταβλητής χρησιμοποιούνται περισσότερες από μία ανεξάρτητες μεταβλητές που ενδεχομένως το επηρεάζουν. Πρώτη προϋπόθεση για την έγκυρη ερμηνεία της μεθόδου είναι η κανονικότητα των υπολοίπων (residuals) του γραμμικού μοντέλου που θα προκύψει μετά το τέλος της διαδικασίας.

Επιπλέον, προϋποτίθεται η ομοσκεδαστικότητα των παρατηρήσεων δηλαδή η ομοιογένεια της διακύμανσης των τιμών της εξαρτημένης μεταβλητής ανάμεσα στις τιμές των ανεξάρτητων μεταβλητών. Τέλος, η εφαρμογή της μεθόδου έχει νόημα όταν η σχέση μεταξύ της εξαρτημένης και των ανεξαρτήτων είναι γραμμική.

Παράδειγμα
Ένας καθηγητής θέλει να δημιουργήσει γραμμικό μοντέλο πρόβλεψης της τελικής επίδοσης των φοιτητών του τμήματός του από την ενδιάμεση εξέταση προόδου, το φύλο και τις απουσίες από το μάθημα. Για το σκοπό αυτό θα αξιοποιήσει ιστορικά στοιχεία για 14 φοιτητές για τους οποίους κατέγραψε την επίδοση στην Πρόοδο (μεταβλητή test), την τελική επίδοση του φοιτητή (μεταβλητή exam), το φύλο του φοιτητή (μεταβλητή gender) και τις απουσίες του φοιτητή από τις διαλέξεις του μαθήματος (μεταβλητή absences).

test = c(5.5, 6, 7, 5, 8, 2, 9, 10, 2, 3, 4, 6.5, 8.5, 1)
exam = c(6.5, 6, 8, 7.5, 7, 4, 8, 10, 1, 5, 5, 6, 9, 5)
gender = c(0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1)
gender = factor(gender, levels = c(0, 1), labels = c("Γυναίκα", "Άνδρας"))
absences = c(2, 3, 1, 1, 2, 5, 0, 0, 7, 5, 4, 3, 0, 4)
exam.df = data.frame(test, exam, gender, absences)
Hmisc::label(exam.df$test) = 'Βαθμός προόδου'
Hmisc::label(exam.df$exam) = 'Βαθμός εξετάσεων'
Hmisc::label(exam.df$gender) = 'Φύλο'
Hmisc::label(exam.df$absences) = 'Απουσίες'

Με τις παρακάτω εντολές υπολογίζουμε τους συντελεστές συσχέτισης Pearson μεταξύ των τριών συνεχών μεταβλητών test, exam και absences:

my_cor_table(exam.df, c('exam', 'test', 'absences'))

	Βαθμός εξετάσεων	Βαθμός προόδου	Απουσίες
Pearson Correlation Coefficients (:p<.05, :p<.01, **:p<.001)
Βαθμός εξετάσεων	-	.853***	-.967***
Βαθμός προόδου	.853***	-	-.866***
Απουσίες	-.967***	-.866***	-

Τα αποτελέσματα μας επιτρέπουν να κάνουμε δύο παρατηρήσεις, μία ενθαρρυντική και μία αποθαρρυντική για τη συνέχιση της διαδικασίας.

Η ενθαρρυντική παρατήρηση είναι πως η εξαρτημένη μεταβλητή (Εξετάσεις) είναι ισχυρά γραμμικά συσχετισμένη με τις δύο συνεχείς ανεξάρτητες (Πρόοδος, Απουσίες) που θα συμμετάσχουν στο γραμμικό μοντέλο, γεγονός που επιβεβαιώνει την καταλληλότητα ενός γραμμικού μοντέλου του βαθμού στις εξετάσεις από το βαθμό στην Πρόοδο και το σύνολο των απουσιών:

• Συντελεστής συσχέτισης Εξετάσεις με Πρόοδο r = 0,853,
• Συντελεστής συσχέτισης Εξετάσεις με Απουσίες r = - 0,967,

Η αποθαρρυντική παρατήρηση είναι πως ο συντελεστής συσχέτισης μεταξύ των δύο ανεξάρτητων παρατηρήσεων Πρόοδο και Απουσίες που σκοπεύουμε να εισάγουμε στη γραμμική εξίσωση είναι ίσος με r = -0,866, δηλαδή οι δύο ανεξάρτητες παρατηρήσεις έχουν ισχυρή γραμμική σχέση μεταξύ τους, γεγονός που ονομάζεται συγγραμικότητα (collinearity) και δημιουργεί σοβαρό πρόβλημα στον υπολογισμό των συντελεστών του μοντέλου.

Στο σημείο αυτό, σε μία πραγματική έρευνα θα έπρεπε να επιλεχθεί για τη συνέχεια μόνο μία από τις δύο συγγραμικές ανεξάρτητες μεταβλητές, ωστόσο για παιδαγωγικούς λόγους θα συνεχίσουμε τη διαδικασία αφήνοντας την R με τους κατάλληλους δείκτες να μας ενημερώσει για τους πιθανούς κινδύνους.

Για τον υπολογισμό του γραμμικού μοντέλου αρκεί να εκτελέσουμε τις εντολές

my.model = my_lm(exam.df, exam ~ test + gender + absences)
my.model$htmlreport

					95% C.I.
Linear regression results for model exam ~ const + test + gender + absences
	B	SE	t	p	Lower	Upper
Constant	8.466	1.101	7.688	<0.001	6.012	10.919
Βαθμός προόδου	0.052	0.125	0.417	0.686	-0.226	0.33
Φύλο	0.389	0.47	0.827	0.428	-0.659	1.437
Απουσίες	-1.018	0.179	-5.696	<0.001	-1.416	-0.62

Coefficient of determination
R² = 0.94

Homoscedasticity Breusch–Pagan test
x²(1) = 2.197, p = 0.533. Homoscedasticity assumption is not rejected. Model seems to be valid.

	Statistic	Description
Normality tests of models’ residuals
Skewness & Kurtosis
Skewness	0.398	Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
Kurtosis	3.099	Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η₀: This sample is from a normal distribution vs Η₁: Not the Η₀.
Shapiro–Wilk	W = 0.91, p = 0.158	This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
Lilliefors	D = 0.225, p = 0.052	The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.

Η επιρροή των παραγόντων αναπαριστάται διαγραμματικά από το ακόλουθο διάγραμμα

library(effects)
fit = lm(exam ~ test + gender + absences, data = exam.df)
# fit = my.model$fit
plot(effects::allEffects(fit))

Διάθεσιμο είναι επιπλέον και το εξής διάγραμμα:

library(sjPlot)    
plot_model(fit, axis.labels=c("Απουσίες", "Φύλο", 'Πρόοδος'), show.values=TRUE, show.p=TRUE, title="Παράγοντες μοντέλου πρόβλεψης τελικού βαθμού")

Ερμηνεία του output

H συνολική αποτελεσματικότητα του γραμμικού μοντέλου περιγράφεται από το συντελεστή προσδιορισμού R square. Παρατηρούμε πως R² = 0,940 = 94% δηλαδή το 94% της συνολικής μεταβλητότητας της τελικής επίδοσης ενός φοιτητή ερμηνεύεται από την επίδοσή του στην Πρόοδο, το φύλο του και τις απουσίες του από το μάθημα στη διάρκεια του εξαμήνου. Ανάλογα, το υπόλοιπο 0,06 = 6% οφείλεται σε άλλους παράγοντες που δεν καταγράφονται στην έρευνά μας.

Η δοκιμασία της ανάλυσης διακύμανσης (ANOVA) μας διαβεβαιώνει πως το σύνολο της μεταβλητότητας του τελικού βαθμού εξέτασης που ερμηνεύεται από το γραμμικό μοντέλο που ορίσαμε είναι στατιστικά σημαντικός F(3, 10) = 51,8, p < 0,001.

Το γραμμικό μοντέλο πρόβλεψης παρουσιάζεται στον τμήμα του output «Coefficients».

Η εξίσωσης πρόβλεψης και κάποια σχόλια…
Στη στήλη Estimate παρουσιάζονται οι συντελεστές του μοντέλου. Η εξίσωση πρόβλεψης που προκύπτει από τα δεδομένα μας είναι:

Εξετάσεις = 8,466 + 0,052 * Πρόοδος + 0,389 * Φύλο - 1,018 *Απουσίες

Καθώς κάποια/ες από τις ανεξάρτητες μεταβλητές που συμμετέχουν στο γραμμικό μοντέλο δεν έχουν σημαντικά διαφορετικούς από το μηδέν συντελεστές (στον πίνακα Coefficients, το αντίστοιχο p είναι μεγαλύτερο από 0,05), το γραμμικό μοντέλο μπορεί να χρησιμοποιηθεί αποκλειστικά για την πρόβλεψη των τιμών της εξαρτημένης μεταβλητής. Με απλά λόγια, η διαδικασία προβλέπει με επιτυχία τις τιμές της εξαρτημένης αλλά δεν μπορεί να εξηγήσει τη σχέση της εξαρτημένης με κάθε μία από τις ανεξάρτητες.

Για παράδειγμα, η πρόβλεψη για την τελική επίδοση στις εξετάσεις για έναν φοιτητή που είχε επίδοση 7 στην πρόοδο, έχει 2 απουσίες και είναι αγόρι (gender = 1) είναι:

Εξετάσεις = 8,466 + 0,052 * 7+ 0,389 * 1 - 1,018 *2 ≈ 7,2

Αν είναι επιθυμητή η εξήγηση της επιρροής των ανεξάρτητων μεταβλητών στην εξαρτημένη και δεν είναι απαραίτητη η συμπερίληψη όλων των μεταβλητών στις ανεξάρτητες τότε ο ερευνητής μπορεί:

(α) να χρησιμοποιήσει κάποια σταδιακή διαδικασία δημιουργίας του γραμμικού μοντέλου (όπως π.χ. τη Stepwise Linear Regression) για να κρατήσει τις περισσότερο σημαντικές από τις ανεξάρτητες μεταβλητές.

(β) να εφαρμόσει ανάλυση σε κύριες συνιστώσες (Principal Component Analysis - PCA) και να ορίσει με τον τρόπο αυτό έναν ή περισσότερους παράγοντες (factors) από γραμμικούς συνδυασμούς των ισχυρά συσχετισμένων ανεξαρτήτων μεταβλητών οι οποίοι παράγοντες θα είναι εκ κατασκευής ασυσχέτιστοι και θα συμμετάσχουν στο γραμμικό μοντέλο στη θέση των αρχικών μεταβλητών.

4.1 Προϋποθέσεις εγκυρότητας

Το ιστόγραμμα των υπολοίπων παράγεται με την εξής εντολή:

my_histFrequency_with_normal_curve(fit$residuals, main = 'Ιστόγραμμα υπολοίπων')

Κάποια επιπλέον διαγράμματα με τα οποία αξιολογείται η ποιότητα του γραμμικού μοντέλου και που είναι απαραίτητα στην περίπτωση όπου υπάρχουν περισσότερες από μία επεξηγηματικές μεταβλητές δίνονται με την εντολή plot(fit) η οποία παράγει 4 διαγράμματα:

layout(matrix(c(1,2,3,4),2,2))
plot(fit)

4.1.1 Διάγραμμα Normal Q – Q

Το διάγραμμα Quantile - Quantile φέρνει σε αντιπαράθεση τις τιμές των υπολοίπων με αυτές που θα περιμέναμε να υπάρχουν αν τα υπόλοιπα ακολουθούν κανονική κατανομή. Επιθυμητό είναι τα σημεία του διαγράμματος να βρίσκονται πάνω στην ευθεία y = x, δηλαδή την ευθεία όπου κάθε της σημείο έχει ίδιες συντεταγμένες. Σημαντική απόκλιση των σημείων από την ευθεία υποδεικνύει μη κανονική κατανομή των υπολοίπων.

Μία καλύτερη ιδέα για τον τρόπο που “διαβάζεται” το διάγραμμα Quantile - Quantile δίνεται στο παρακάτω διάγραμμα: Πηγή: Glen_b (https://stats.stackexchange.com/users/805/glen-b), How to interpret a QQ plot, URL (version: 2021-09-04): https://stats.stackexchange.com/q/101290)

Για την καλύτερη κατανόηση του τρόπου που μετασχηματίζεται το Q-Q plot καθώς η κατανομή απομακρύνεται από την κανονική δίνεται εδώ.

Στην περίπτωσή μας, παρατηρούμε πως υπάρχουν τρεις παρατηρήσεις (στις γραμμές των δεδομένων 7, 8 και 10) οι οποίες δημιουργούν ένα μικρό πρόβλημα στην κανονικότητα των υπολοίπων.

4.1.2 Διάγραμμα Residuals vs Fitted

Το διάγραμμα Residuals vs Fitted απεικονίζει την διασπορά των υπολοίπων σε όλο το εύρος τιμών που προβλέπει το γραμμικό μοντέλο. Για να είναι έγκυρη η χρήση του γραμμικού μοντέλου πρέπει τα υπόλοιπα να έχουν την ίδια διασπορά σε όλο το μήκος τιμών του οριζόντιου άξονα. Για να γίνει πιο εύκολη η ανάγνωση βρίσκεται η τάση των σημείων με την μέθοδο της τοπικής παλινδρόμησης (Local Regression - LOWESS) και εμφανίζεται με κόκκινο χρώμα. Επιθυμητό είναι η γραμμή αυτή να είναι περίπου οριζόντια.

Μία καλύτερη ιδέα για τον τρόπο που “διαβάζεται” το διάγραμμα Residuals vs Fitted δίνεται στο παρακάτω διάγραμμα: Ενδεικτικά διαγράμματα Residuals vs Fitted Πηγή: Linear Models with R (Chapman & Hall/CRC Texts in Statistical Science)

Στο πρώτο διάγραμμα παρατηρούμε πως τα υπόλοιπα έχουν την ίδια κατανομή σε όλο το εύρος των τιμών πρόβλεψης του μοντέλου, γεγονός που υποδεικνύει πως σε όλο το εύρος των τιμών της ανεξάρτητης μεταβλητής το μοντέλο προσφέρει προβλέψεις με την ίδια “ποιότητα”.

Στο δεύτερο διάγραμμα, το εύρος των υπολοίπων αυξάνεται καθώς οι προσαρμοσμένες τιμές (ή ισοδύναμα το x) αυξάνει. Αυτό το φαινόμενο ονομάζεται ετεροσκεδαστικότητα και υποδεικνύει πως το γραμμικό μοντέλο δεν κάνει την ίδια καλή πρόβλεψη σε όλο το εύρος τιμών του x, ή ισοδύναμα πως υπάρχουν και άλλοι παράγοντες που πρέπει να εισαχθούν στο μοντέλο. Ειδικότερα, η ερμηνεία των συντελεστών του μοντέλου είναι επισφαλής.

Στο τρίτο διάγραμμα τα υπόλοιπα είναι ως επί το πλείστον αρνητικά όταν η προσαρμοσμένη τιμή είναι μικρή, θετικά όταν η προσαρμοσμένη τιμή είναι στη μέση και αρνητικά όταν η προσαρμοσμένη τιμή είναι μεγάλη. Αυτό σημαίνει πως οι πραγματικές τιμές της εξαρτημένης είναι κατά σειρά κάτω - πάνω - κάτω από την ευθεία πρόβλεψης υποδεικνύοντας μία μη γραμμική (πιθανώς δεύτεροβάθμια) σχέση με την ανεξάρτητη μεταβλητή. Η ερμηνεία των συντελεστών του μοντέλου και στην περίπτωση αυτή είναι επισφαλής.

Στην περίπτωσή μας, η προϋπόθεση της ομοσκεδαστικότητας μπορεί να γίνει αποδεκτή για τα δεδομένα μας καθώς όπως παρατηρείται από το αντίστοιχο διάγραμμα διασποράς (Residuals vs Fitted) τα σημεία έχουν ανάλογη κύμανση για όλο το εύρος τιμών της τιμής που προβλέπεται από το γραμμικό μοντέλο.

4.1.3 Διάγραμμα Residuals vs Leverage

Η απόσταση Cook δείχνει το πόσο πολύ θα μεταβαλλόταν τα υπόλοιπα των άλλων παρατηρήσεων αν απομακρυνόταν η ίδια η παρατήρηση. Όσο μεγαλύτερη είναι η τιμή της τόσο πιο ισχυρή είναι η επιρροή του σημείου στον υπολογισμό των συντελεστών του μοντέλου. Από το διάγραμμα Residuals vs Leverage (υπόλοιπα παρατήρησης vs ικανότητα μόχλευσης) παρατηρούμε ότι υπάρχουν δεν υπάρχουν σημεία με σημαντική ικανότητα μόχλευσης, δηλαδή σημεία με απόσταση Cook μεγαλύτερη από το 1.

4.1.4 Διάγραμμα Scale-Location

Το Scale-Location διάγραμμα είναι ακόμα ένα διάγραμμα με στόχο την ανίχνευση της ετεροσκεδαστικότητας. Στο διάγραμμα αυτό, συγκρίνονται οι τιμές που προβλέπονται από το μοντέλο με τα τυποποιημένα υπόλοιπα, έχοντας μεγαλύτερη αξιοπιστία από το Διάγραμμα Residuals vs Fitted καθώς το τελευταίο επηρεάζεται από την μορφή της κατανομής της εξαρτημένης μεταβλητής.

Και σε αυτό, επιθυμητό είναι η κόκκινη γραμμή (LOWESS) να είναι περίπου οριζόντια (το μέσο μέγεθος των τυποποιημένων υπολειμμάτων δεν αλλάζει πολύ ως συνάρτηση των προσαρμοσμένων τιμών) και το εύρος γύρω από την κόκκινη γραμμή να μην διαφέρει κατά μήκος των προσαρμοσμένων τιμών.

Στην περίπτωσή μας, από το διάγραμμα Scale-Location υποδεικνύεται ένα πρόβλημα ετεροσκεδαστικότητας που οφείλεται στις παρατηρήσεις 8, 9 και 10.

Περισσότερες λεπτομέρειες για το είδος των διαγραμμάτων που προκύπτει από την plot(fit) παρουσιάζονται εδώ ή εδώ

Από τις παραπάνω παρατηρήσεις, ο ερευνητής μπορεί να μπει στον πειρασμό να εξαιρέσει την 9η παρατήρηση από τα δεδομένα και να επαναλάβει την διαδικασία, ωστόσο είναι προς όφελος της χρηστικότητας του μοντέλου να μην το κάνουμε αυτό καθώς η συγκεκριμένη παρατήρηση αντιπροσωπεύει ένα πραγματικό ζεύγος παρατηρήσεων που μπορεί να εμφανιστεί ξανά στο μέλλον και η διαγραφή του θα επέβαλε την ερμηνεία του μοντέλου σε μικρότερο εύρος τιμών μειώνοντας έτσι την σημαντικότητά του.

Ανοχή (Tolerance) και συγγραμικότητα.

Από τις στήλες t, p, φανερώνεται μία μη επιθυμητή πραγματικότητα. Οι συντελεστές των μεταβλητών Φύλο και Πρόοδος δεν προκύπτουν να είναι στατιστικώς σημαντικοί (t(12) = 0,827, p = 0,428 και t(12) = 0,417, p = 0,686 αντίστοιχα), γεγονός που οφείλεται στην έντονη συγγραμικότητα μεταξύ των μεταβλητών Πρόοδος και Απουσίες που παρατηρήσαμε στην αρχή της διαδικασίας με το συντελεστή συσχέτισης του Pearson.

Η συγγραμικότητα των μεταβλητών μπορεί να μετρηθεί με το παράγοντα μεγέθυνσης της κύμανσης (Variance Inflation Factor - VIF). Στην R, αυτός είναι εύκολο να υπολογιστεί με την εντολή vif της βιβλιοθήκης συναρτήσεων car

library(car)
vif(fit)

##     test   gender absences 
## 3.994111 1.862616 4.799177

Η ανοχή (Tolerance) μίας ανεξάρτητης μεταβλητής είναι το ποσοστό της μεταβλητότητάς της που δεν εξηγείται από τις υπόλοιπες ανεξάρτητες μεταβλητές. Φανερά, όσο μεγαλύτερη είναι η τιμή αυτή για μία ανεξάρτητη μεταβλητή τόσο το καλύτερο για αυτήν καθώς αυτό σημαίνει πως αποτελεί ανεξάρτητη οντότητα στη γραμμική εξίσωση. Αντίθετα, μία μικρή τιμή σε αυτόν το δείκτη σημαίνει πως η μεταβλητή είναι περιττή και πως μπορεί να βγει από την ομάδα των ανεξαρτήτων μεταβλητών.

Ο παράγοντας πληθωρισμού της κύμανσης ορίζεται ως ο αντίστροφος της ανοχής (VIF = 1 / Tolerance) και ερμηνεύεται με αντίστροφο τρόπο. Δεν υπάρχει κάποιο γενικά αποδεκτό όριο για το μέγεθος του VIF που να καταδεικνύει σημαντικό πρόβλημα συγγραμικότητας, με πολλούς ερευνητές να χρησιμοποιούν το 10 ως το μέγιστό όριο του παράγοντα αυτού που αντιστοιχεί σε ελάχιστη τιμή 0,1 για την ανοχή του ίδιου παράγοντα.

4.2 Μέθοδος Stepwise σταδιακής εισαγωγής – εξαγωγής

Με την προϋπόθεση πως ο χρήστης δεν θεωρεί απαραίτητη τη συμμετοχή στο γραμμικό μοντέλο όλων των μεταβλητών που αρχικά ορίζει, είναι δυνατή η επιλογή μίας από τις βαθμιδωτές μεθόδους που αποφασίζουν τη συμπερίληψη ή απόρριψη κάθε μίας ανεξάρτητης μεταβλητής στο γραμμικό μοντέλο βάσει του αποτελέσματος ενός στατιστικού κριτηρίου. Είναι σημαντικό όμως ο χρήστης να συγκρατήσει πως όποιο και αν είναι το στατιστικό κριτήριο που θα χρησιμοποιηθεί, δεν είναι βέβαιο πως το τελικό αποτέλεσμα της διαδικασίας θα είναι ικανοποιητικό.

Στην R, υπάρχουν πολλές εντολές που υλοποιούν τη βαθμιαία δημιουργία μοντέλου. Στην ενότητα αυτή παρουσιάζεται η εντολή:

step(fit, direction="both") # ή direction = “forward” ή direction = “backward”

## Start:  AIC=-9.29
## exam ~ test + gender + absences
## 
##            Df Sum of Sq     RSS      AIC
## - test      1    0.0708  4.1432 -11.0462
## - gender    1    0.2782  4.3506 -10.3623
## <none>                   4.0724  -9.2875
## - absences  1   13.2123 17.2848   8.9507
## 
## Step:  AIC=-11.05
## exam ~ gender + absences
## 
##            Df Sum of Sq    RSS      AIC
## - gender    1     0.273  4.417 -12.1517
## <none>                   4.143 -11.0462
## + test      1     0.071  4.072  -9.2875
## - absences  1    38.065 42.208  19.4498
## 
## Step:  AIC=-12.15
## exam ~ absences
## 
##            Df Sum of Sq    RSS     AIC
## <none>                   4.417 -12.152
## + gender    1     0.273  4.143 -11.046
## + test      1     0.066  4.351 -10.362
## - absences  1    62.941 67.357  23.993

## 
## Call:
## lm(formula = exam ~ absences, data = exam.df)
## 
## Coefficients:
## (Intercept)     absences  
##       8.966       -1.014

Ερμηνεία του output
Η εντολή step, υπολογίζει το δείκτη AIC (Akaike Information Criterion) για όλα τα πιθανά γραμμικά μοντέλα, επιλέγοντας τελικά το μοντέλο που μειώνει στο μεγαλύτερο βαθμό την τιμή του. Με το κριτήριο αυτό, η διαδικασία τελικά προτείνει μοντέλο πρόβλεψης με μοναδική ανεξάρτητη μεταβλητή τις απουσίες (absences). Η εξίσωση πρόβλεψης είναι η:

Εξετάσεις = 8,966 - 1,014 * Απουσίες,

και η ερμηνεία των συντελεστών είναι άμεση: Ένας φοιτητής δίχως απουσίες στη διάρκεια των μαθημάτων αναμένεται να έχει βαθμό 8,966 ≈ 9 στην τελική εξέταση ενώ για κάθε απουσία που πραγματοποιείται αναμένεται μείωση στο βαθμό κατά 1,014 ≈ 1 μονάδα.

5 Πολλαπλή γραμμική παλινδρόμηση για την πρόβλεψη χρονοσειρών

Η πολλαπλή γραμμική παλινδρόμηση μπορεί να εφαρμοστεί για την εκτίμηση των συντελεστών μοντέλων της μορφής. Ωστόσο, εκτός από τις γενικές προϋποθέσεις της γραμμικής παλινδρόμησης (ομοσκεδάστικότητα, κανονικά υπόλοιπα) χρειάζεται προσοχή στις εξής δύο προϋποθέσεις:

Προϋπόθεση 1
Αν η εξαρτημένη και μία επεξηγηματική σειρά έχουν τάση, τότε οπωσδήποτε θα έχουν υψηλή συσχέτιση, ανεξάρτητα από το αν έχουν μεταξύ τους συσχέτιση, με αποτέλεσμα να προκύπτει ψευδώς σημαντική αποτελεσματικότητα του μοντέλου (υψηλό R², μικρά p, κλπ). Συνεπώς, μία προϋπόθεση εγκυρότητας της μεθόδου των ελαχίστων τετραγώνων για χρονοσειρές, είναι αυτή να εφαρμόζεται αφού πρώτα αφαιρεθεί η τάση από τις εμπλεκόμενες σειρές. Στην περίπτωση όπου αυτό συμβεί τότε η αφαίρεση της τάσης πρέπει να γίνει τόσο από την εξαρτημένη όσο και από την ανεξάρτηση σειρά, ώστε το μοντέλο που θα προκύψει να έχει νόημα.

Προϋπόθεση 2
Αν η εξαρτημένη και οι ανεξάρτητες μεταβλητές είναι χρονοσειρές, τότε είναι πολύ πιθανό τα υπόλοιπα του μοντέλου να έχουν και αυτά τη δομή μίας χρονοσειράς, δηλαδή να είναι εξαρτημένα, γεγονός που επηρεάζει την εγκυρότητα των υπολογισμών στη μέθοδο των ελαχίστων τετραγώνων, οδηγώντας πιο εύκολα σε στατιστικώς σημαντικούς συντελεστές.

Αν τα υπόλοιπα είναι συσχετισμένα, τότε μπορούμε:

α) Να εφαρμόσουμε ένα μοντέλο ARIMA μόνο για αυτά και το γραμμικό μοντέλο της επιλογής μας για τη χρονοσειρά (εφικτό με τη συνάρτηση Arima της βιβλιοθήκης forecast ή τη συνάρτηση arimax (ARIMA with eXplanatory variables) της βιβλιοθήκης TSA.). Στην περίπτωση αυτή πρέπει όλες οι μεταβλητές που συμμετέχουν να είναι στάσιμες (ενδεικτική αναφορά: https://otexts.com/fpp2/dynamic.html). Το μοντέλο που ελέγχεται είναι της μορφής:

y_t = β x_t + ε_t
ε_t = φ₁ ε_t-1 + … + φ_p ε_t-p - θ₁ w_t-1 - … - θ_q w_t-q, w_t ~ N(0, σ²).

β) Nα προσαρμόσουμε ένα ARIMA μοντέλο σε αυτά, έχοντας προσθέσει και την επεξηγηματική μεταβλητή. Η μέθοδος αυτή είναι εφικτή με τη συνάρτηση arimax (ARIMA with eXplanatory variables) της βιβλιοθήκης TSA. Στην περίπτωση αυτή, το μοντέλο που ελέγχεται είναι της μορφής:

y_t = β x_t + φ₁ y_t-1 + … + φ_p y_t-p - θ₁ ε_t-1 + … - θ_q ε_t-q, ε_t ~ N(0, σ²).

5.1 Περίπτωση γραμμικού μοντέλου με τάση και εποχικότητα

Αν επιθυμούμε να ανιχνεύσουμε την τάση με τη μέθοδο γραμμικής παλινδρόμησης, τότε αρκεί να προσθέσουμε έναν όρο της μορφής α*t στο μοντέλο μας.

5.1.1 Εντοπισμός τάσης

Δίνεται η χρονοσειρά 47 τιμών (όγκος απορριμμάτων Δήμου Κατερίνης από το μήνα Νοέμβριο του 2003 έως τον Σεπτέμβριο του 2007):

ts.value = c(1661980, 1357840, 1678840, 1859050, 1859540, 1814700, 2280190, 2240690, 2072170, 1872380, 1693170, 1857710, 1609670, 1523130, 1919710, 1879330, 1723990, 2160040, 2318030, 2308970, 2092730, 1975510, 1907070, 1904990, 1686480, 1511710, 1840750, 2083950, 1946310, 1917100, 2315390, 2393510, 2087530, 1961500, 1807630, 1877010, 1709370, 1621490, 1946690, 2010600, 2069190, 2031980, 2318730, 2451960, 2213610, 2196030, 1880220)

Δημιουργούμε ένα διάνυσμα όπου θα καταχωρηθεί η χρονική στιγμή:

time.var = c(1:length(ts.value))

Τώρα, αρκεί να εφαρμόσουμε γραμμικό μοντέλο πρόβλεψης με την εντολή:

katerini.df = data.frame(ts.value, time.var)
fit.katerini = my_lm(katerini.df, ts.value ~ time.var)
fit.katerini$htmlreport

					95% C.I.
Linear regression results for model ts.value ~ const + time.var
	B	SE	t	p	Lower	Upper
Constant	1787097.28	70167.333	25.469	<0.001	1645773.016	1928421.544
time.var	6610.459	2545.237	2.597	0.013	1484.089	11736.829

Coefficient of determination
R² = 0.13

Homoscedasticity Breusch–Pagan test
x²(1) = 0.305, p = 0.581. Homoscedasticity assumption is not rejected. Model seems to be valid.

	Statistic	Description
Normality tests of models’ residuals
Skewness & Kurtosis
Skewness	0.128	Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
Kurtosis	2.477	Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η₀: This sample is from a normal distribution vs Η₁: Not the Η₀.
Shapiro–Wilk	W = 0.968, p = 0.216	This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
Lilliefors	D = 0.104, p = 0.232	The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.

Πρόβλεψη μοντέλου:

predict(fit, newdata = data.frame(time.var = c(48:72)), interval = "confidence")

##         fit       lwr       upr
## 1  6.715980 6.0472520  7.384707
## 2  6.112968 5.5118272  6.714109
## 3  7.811683 7.2197776  8.403589
## 4  7.707763 6.8204771  8.595048
## 5  6.845880 6.0568233  7.634937
## 6  3.869601 3.2467040  4.492498
## 7  8.933367 8.2892405  9.577493
## 8  8.985327 8.2348642  9.735790
## 9  1.834074 0.8339787  2.834170
## 10 3.921561 3.3453627  4.497760
## 11 4.991285 4.4875948  5.494974
## 12 6.138948 5.4684695  6.809427
## 13 9.296159 8.2116121 10.380705
## 14 4.835404 3.8655561  5.805251

Μία γραφική αναπαράσταση μπορεί να γίνει με τον εξής κώδικα:

fit = tslm(ts(ts.value) ~ trend)
plot(forecast(fit, h=20), showgap = F)

Παρατηρούμε πως για κάθε έναν μήνα που περνάει, αναμένουμε αύξηση 6.610 κιλά στα απορρίματα του Δήμου Κατερίνης, αύξηση που είναι στατιστικώς σημαντική (p = 0.0127). Είναι φανερό πως πέρα από αυτήν την πληροφορία το μοντέλο δεν προσφέρει κάποια άλλη υπηρεσία καθώς όπως αναμενόταν τα υπόλοιπα είναι συσχετισμένα:

acf(fit.katerini$fit$residuals)

5.1.2 Εντοπισμός εποχικότητας

Γνωρίζοντας πως το φαινόμενο είναι περιοδικό με περίοδο 12 μπορούμε να προσαρμόσουμε 11 ψευδομεταβλητές οι οποίες θα ανιχνεύσουν διαφοροποίηση μεταξύ του πρώτου μήνα και όλων των υπολοίπων.

Ενδεικτικά, στο παράδειγμα των απορριμάτων του Δήμου Κατερίνης, αυτό μπορεί να γίνει με τον παρακάτω κώδικα:

Feb = rep(c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), length.out=length(ts.value))
Mar = rep(c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0), length.out=length(ts.value))
Apr = rep(c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0), length.out=length(ts.value))
May = rep(c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0), length.out=length(ts.value))
Jun = rep(c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0), length.out=length(ts.value))
Jul = rep(c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0), length.out=length(ts.value))
Aug = rep(c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0), length.out=length(ts.value))
Sep = rep(c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0), length.out=length(ts.value))
Oct = rep(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0), length.out=length(ts.value))
Nov = rep(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0), length.out=length(ts.value))
Dec = rep(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), length.out=length(ts.value))

katerini.df = data.frame(ts.value, time.var, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec)

fit.katerini = my_lm(katerini.df, ts.value ~ Feb+Mar+Apr+May+Jun+Jul+Aug+Sep+Oct+Nov+Dec)
fit.katerini$htmlreport

					95% C.I.
Linear regression results for model ts.value ~ const + Feb + Mar + Apr + May + Jun + Jul + Aug + Sep + Oct + Nov + Dec
	B	SE	t	p	Lower	Upper
Constant	1666875	51693.889	32.245	<0.001	1561930.826	1771819.174
Feb	-163332.5	73106.199	-2.234	0.032	-311745.974	-14919.026
Mar	179622.5	73106.199	2.457	0.019	31209.026	328035.974
Apr	291357.5	73106.199	3.985	<0.001	142944.026	439770.974
May	232882.5	73106.199	3.186	0.003	84469.026	381295.974
Jun	314080	73106.199	4.296	<0.001	165666.526	462493.474
Jul	641210	73106.199	8.771	<0.001	492796.526	789623.474
Aug	681907.5	73106.199	9.328	<0.001	533494.026	830320.974
Sep	449635	73106.199	6.15	<0.001	301221.526	598048.474
Oct	334480	73106.199	4.575	<0.001	186066.526	482893.474
Nov	155147.5	73106.199	2.122	0.041	6734.026	303560.974
Dec	213028.333	78963.72	2.698	0.011	52723.46	373333.207

Coefficient of determination
R² = 0.871

Homoscedasticity Breusch–Pagan test
x²(1) = 13.438, p = 0.266. Homoscedasticity assumption is not rejected. Model seems to be valid.

	Statistic	Description
Normality tests of models’ residuals
Skewness & Kurtosis
Skewness	0.023	Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
Kurtosis	2.68	Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η₀: This sample is from a normal distribution vs Η₁: Not the Η₀.
Shapiro–Wilk	W = 0.982, p = 0.665	This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
Lilliefors	D = 0.077, p = 0.684	The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.

Παρατηρούμε πως η περιοδική φύση της σειράς αντανακλάστηκε σε στατιστικά σημαντικούς συντελεστές που προσδιορίζουν και τη διαφορά μεταξύ των 11 μηνών με τον Ιανουάριο.

Ελέγχουμε τη συσχέτιση των υπολοίπων:

acf(fit.katerini$fit$residuals)

Παρατηρούμε πως υπάρχουν μερικές θετικές αυτοσυσχετίσεις. Αυτό οφείλεται στο γεγονός πως το μοντέλο που ελέγξαμε δεν θεραπεύει το πρόβλημα της τάσης της σειράς.

5.1.2.1 Ποιος είναι ο λόγος που λείπει ο Ιανουάριος από το μοντέλο;

Στην εύλογη απορία σχετικά με το λόγο που δεν υπάρχει ο Ιανουάριος στο μοντέλο, η απάντηση είναι πως σε αυτήν την περίπτωση θα είχαμε μία ομάδα αποκρισεων που θα είχαν τέλεια γραμμική σχέση:

Jan + Feb + Mar + Apr + May + Jun + Jul + Aug + Sep + Oct + Nov + Dec = 1.

Αυτό προκαλεί πρόβλημα στον υπολογισμό των συντελεστών με την μέθοδο των ελαχίστων τετραγώνων, γιατί η διαδικασία βασίζεται στην αντιστροφή του πίνακα X^TX, όπου ο Χ έχει διάσταση n x (k + 1) με n το πλήθος των παρατηρήσεων (47 στο παράδειγμα) και k το πλήθος των επεξηγηματικών μεταβλητών. Μαζί με τον Ιανουάριο, ο πίνακας αυτός θα είχε 13 στήλες που θα ήταν γραμμικώς εξαρτημένες, συνεπώς ούτε αυτός, ούτε ο X^TX θα μπορούσε να αντιστραφεί και η διαδικασία δεν ολοκληρώνεται.

Σχετική αναφορά: https://en.wikipedia.org/wiki/Multicollinearity#Perfect_multicollinearity.

Διαισθητικά, το πρόβλημα αναδεικνύεται από το γεγονός πως στη γραμμική παλινδρόνηση, ο συντελεστής κάθε μίας μεταβλητής δείχνει την αναμενόμενη επιρροή της μεταβλητής στην εξαρτημένη μεταβλητή όταν όλες οι υπόλοιπες μεταβλητές παραμένουν σταθερές.

Οπότε, αν υπήρχε μέσα και η μεταβλητή Jan, τότε ο συντελεστής της θα έπρεπε να δείχνει την επιρροή του Ιανουαρίου στην εξαρτημένη μεταβλητή, όταν όλες οι υπόλοιπες μεταβλητές είναι σταθερές, κάτι αδύνατον διότι:

Feb + Mar + Apr + May + Jun + Jul + Aug + Sep + Oct + Nov + Dec = 1 - Jan.

5.1.3 Εντοπισμός τάσης και εποχικότητας με την βιβλιοθήκη tslm

Η συνάρτηση tslm της ομώνυμης βιβλιοθήκης εφαρμόζει ευθεία ελαχίστων τετραγώνων όπως και η συνάρτηση lm, με τη διαφορά πως δίνει το δικαίωμα για χρήση των ψευδο-μεταβλητών trend και season που ανιχνεύουν αντιστοιχα την τάση και την εποχικόττηα μίας σειράς. Για να εφαρμοστεί η συνάρτηση αυτή, πρέπει τα δεδομένα να έχουν δηλωθεί στην μορφή χρονοσειράς με την εντολή ts και να έχει δηλωθεί η περίοδός τους.

ts.value = c(1661980, 1357840, 1678840, 1859050, 1859540, 1814700, 2280190, 2240690, 2072170, 1872380, 1693170, 1857710, 1609670, 1523130, 1919710, 1879330, 1723990, 2160040, 2318030, 2308970, 2092730, 1975510, 1907070, 1904990, 1686480, 1511710, 1840750, 2083950, 1946310, 1917100, 2315390, 2393510, 2087530, 1961500, 1807630, 1877010, 1709370, 1621490, 1946690, 2010600, 2069190, 2031980, 2318730, 2451960, 2213610, 2196030, 1880220)

katerini.trash = ts(ts.value, frequency = 12, start = c(2003, 11))
print(katerini.trash)

##          Jan     Feb     Mar     Apr     May     Jun     Jul     Aug     Sep
## 2003                                                                        
## 2004 1678840 1859050 1859540 1814700 2280190 2240690 2072170 1872380 1693170
## 2005 1919710 1879330 1723990 2160040 2318030 2308970 2092730 1975510 1907070
## 2006 1840750 2083950 1946310 1917100 2315390 2393510 2087530 1961500 1807630
## 2007 1946690 2010600 2069190 2031980 2318730 2451960 2213610 2196030 1880220
##          Oct     Nov     Dec
## 2003         1661980 1357840
## 2004 1857710 1609670 1523130
## 2005 1904990 1686480 1511710
## 2006 1877010 1709370 1621490
## 2007

Εφαρμόζουμε γραμμικό μοντέλο πρόβλεψης με την εντολή:

fit.katerini <- tslm(katerini.trash ~ trend + season)
summary(fit.katerini)

## 
## Call:
## tslm(formula = katerini.trash ~ trend + season)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -147916  -43771  -11961   36094  206937 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1749015.7    41977.0  41.666  < 2e-16 ***
## trend          4642.0      840.7   5.521 3.61e-06 ***
## season2      107093.0    53865.0   1.988   0.0549 .  
## season3       43976.0    53884.7   0.816   0.4201    
## season4      120531.5    53917.5   2.235   0.0321 *  
## season5      443019.5    53963.3   8.210 1.41e-09 ***
## season6      479075.1    54022.2   8.868 2.30e-10 ***
## season7      242160.6    54094.1   4.477 8.11e-05 ***
## season8      122363.6    54179.0   2.259   0.0304 *  
## season9      -61610.9    54276.8  -1.135   0.2643    
## season10      19479.9    58228.4   0.335   0.7400    
## season11    -170338.5    53884.7  -3.161   0.0033 ** 
## season12    -338313.0    53865.0  -6.281 3.73e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 76170 on 34 degrees of freedom
## Multiple R-squared:  0.932,  Adjusted R-squared:  0.9079 
## F-statistic: 38.81 on 12 and 34 DF,  p-value: 2.713e-16

Οπτικοποιούμε την προσαρμογή του μοντέλου με το παρακάτω γράφημα:

library(ggplot2)
autoplot(katerini.trash, series="Δεδομένα") +
  autolayer(fitted(fit.katerini), series="Εκτίμηση") +
  xlab("Μήνας") + ylab("Χιλιάδες Κιλά") +
  ggtitle("Μηνιαία Παραγωγή Απορριμάτων")

Παρατηρούμε τις προβλέψεις του μοντέλου, όπως παρακάτω:

plot(forecast(fit.katerini, h=20), showgap = F)

Χρήσιμη αναφορά: https://roznn.github.io/Forecasting-with-R/LinearRegressionForecasting.html

5.1.4 Παράδειγμα

Στο διάνυσμα ausbeer δίνεται η τριμηνιαία παραγωγή μπύρας στην Αυστραλία (σε μεγάλιτρα) από το 1ο τετράμηνο 1956: Q1 έως 2008: Q3. Θα αξιοποιήσουμε την εντολή tslm της ομώνυμης βιβλιοθήκης για να προσαρμόσουμε ένα μοντέλο πρόβλεψης.

library(fpp)
data(ausbeer)
my.plot.ts.acf(ausbeer)

Θα εστιάσουμε στο (φαινομενικά) στάσιμο μέρος της σειράς, μετά το 1992.

beer2 <- window(ausbeer, start=1992)
my.plot.ts.acf(beer2)

Προσαρμόζουμε το γραμμικό μοντέλο με τη συνάρτηση tslm.

fit.beer <- tslm(beer2 ~ trend + season)
summary(fit.beer)

## 
## Call:
## tslm(formula = beer2 ~ trend + season)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.916  -7.877  -0.070   7.594  21.494 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 442.78341    3.98067 111.233  < 2e-16 ***
## trend        -0.35886    0.07866  -4.562 2.45e-05 ***
## season2     -35.40585    4.26869  -8.294 1.22e-11 ***
## season3     -19.28229    4.27086  -4.515 2.89e-05 ***
## season4      72.79268    4.33485  16.792  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.44 on 62 degrees of freedom
## Multiple R-squared:  0.923,  Adjusted R-squared:  0.918 
## F-statistic: 185.8 on 4 and 62 DF,  p-value: < 2.2e-16

Οπτικοποιούμε την προσαρμογή του μοντέλου με το παρακάτω γράφημα:

library(ggplot2)
autoplot(beer2, series="Data") +
  autolayer(fitted(fit.beer), series="Fitted") +
  xlab("Year") + ylab("Megalitres") +
  ggtitle("Quarterly Beer Production")

cbind(Data=beer2, Fitted=fitted(fit.beer)) %>%
  as.data.frame() %>%
  ggplot(aes(x = Data, y = Fitted,
             colour = as.factor(cycle(beer2)))) +
  geom_point() +
  ylab("Fitted") + xlab("Actual values") +
  ggtitle("Quarterly beer production") +
  scale_colour_brewer(palette="Dark2", name="Quarter") +
  geom_abline(intercept=0, slope=1)

Οι προβλέψεις αναπαριστώνται στο παρακάτω διάγραμμα:

plot(forecast(fit.beer, h=20), showgap = F)

5.2 Περίπτωση AR(p)

Η μέθοδος των ελαχίστων τετραγώνων μπορεί να χρησιμοποιηθεί για την προσαρμογή μοντέλων που περιλαμβάνουν μόνο αυτοπαλινδρομικούς όρους.

5.2.1 Εφαρμογή lm και my_lm

Στο dataframe uschange της βιβλιοθήκης fpp2 παρουσιάζεται (μεταξύ άλλων) η μεταβολή στα έξοδα και στο εισόδημα, για τις ΗΠΑ ανά τετράμηνο από το έτος 1970 έως το έτος 2016.

library(fpp2)
head(uschange)

##         Consumption     Income Production   Savings Unemployment
## 1970 Q1   0.6159862  0.9722610 -2.4527003 4.8103115          0.9
## 1970 Q2   0.4603757  1.1690847 -0.5515251 7.2879923          0.5
## 1970 Q3   0.8767914  1.5532705 -0.3587079 7.2890131          0.5
## 1970 Q4  -0.2742451 -0.2552724 -2.1854549 0.9852296          0.7
## 1971 Q1   1.8973708  1.9871536  1.9097341 3.6577706         -0.1
## 1971 Q2   0.9119929  1.4473342  0.9015358 6.0513418         -0.1

Consumption = as.vector(uschange[,"Consumption"])
plot(Consumption, type = 'l')

acf(Consumption)

Από το acf διάγραμμα, προκύπτει πως ένα μοντέλο που θα μπορούσε να προσαρμοστεί στη χρονοσειρά καλό είναι να περιέχει 3 ΜΑ όρους. Ωστόσο, για παιδαγωγικούς λόγους προχωρούμε στην προσαρμογή ενός AR(2) μοντέλου.

Τις μετατοπισμένες εκδοχές της σειράς τις υπολογίζουμε με τις παρακάτω εντολές.

Consumption.lag1 = c(NA, head(Consumption, -1))
Consumption.lag2 = c(NA, head(Consumption.lag1, -1))
data.lm = data.frame(Consumption, Consumption.lag1, Consumption.lag2)

Τώρα μπορούμε να εφαρμόσουμε μοντέλο πρόβλεψης AR(2).

my.model = my_lm(data.lm, Consumption ~ Consumption.lag1 + Consumption.lag2)
my.model$htmlreport

					95% C.I.
Linear regression results for model Consumption ~ const + Consumption.lag1 + Consumption.lag2
	B	SE	t	p	Lower	Upper
Constant	0.382	0.076	5.029	<0.001	0.232	0.532
Consumption.lag1	0.273	0.072	3.775	<0.001	0.13	0.416
Consumption.lag2	0.219	0.072	3.023	0.003	0.076	0.361

Coefficient of determination
R² = 0.164

Homoscedasticity Breusch–Pagan test
x²(1) = 9.509, p = 0.009. Homoscedasticity assumption is rejected! Read model coefficients with caution.

	Statistic	Description
Normality tests of models’ residuals
Skewness & Kurtosis
Skewness	-0.503	Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
Kurtosis	5.185	Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η₀: This sample is from a normal distribution vs Η₁: Not the Η₀.
Shapiro–Wilk	W = 0.964, p = 0	This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
Lilliefors	D = 0.084, p = 0.003	The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.

Το μοντέλο που προτείνεται είναι το εξής:

Consumption_t = 0.382 + 0.273 Consumption_t-1 + 0.219 Consumption_t-2

Παρατηρούμε πως οι προϋποθέσεις της ομοσκεδαστικότητας και της κανονικότητας της κατανομής των υπολοίπων δεν ισχύουν, άρα το μοντέλο μας δεν προσαρμόζεται με αξιοπιστία στις παρατηρήσεις.

Η έλλειψη προσαρμογής επιβεβαιώνεται και από το διάγραμμα acf των υπολοίπων από το οποίο προκύπτει αυτοσυσχέτιση των υπολοίπων σε lag 3.

acf(my.model$fit$residuals, lag.max = 50)

5.2.2 Εφαρμογή ar.ols

Εναλλακτικά, το μοντέλο AR(2) μπορεί να γίνει προσαρμογή στις τιμές της Consumption με τη συνάρτηση ar.ols της baseR.

ar.ols(Consumption, order.max = 2, demean=F, intercept=T)

## 
## Call:
## ar.ols(x = Consumption, order.max = 2, demean = F, intercept = T)
## 
## Coefficients:
##      1       2  
## 0.2730  0.2187  
## 
## Intercept: 0.3818 (0.07531) 
## 
## Order selected 2  sigma^2 estimated as  0.3604

Σημείωση
Η συνάρτηση ar.ols έχει τη δυνατότητα να αποδόσει το βέλτιστο AR(p) μοντέλο αν εκτελεστεί χωρίς περιορισμό στο order.max.

5.2.3 Εφαρμογή ARIMA

Όπως αναμένεται, το αποτέλεσμα των ελαχίστων τετραγώνων δεν διαφέρει πολύ από το αποτέλεσμα της ARIMA.

library(forecast)
Arima(Consumption, order = c(2, 0, 0))

## Series: Consumption 
## ARIMA(2,0,0) with non-zero mean 
## 
## Coefficients:
##         ar1     ar2    mean
##       0.272  0.2166  0.7458
## s.e.  0.071  0.0710  0.0848
## 
## sigma^2 = 0.3628:  log likelihood = -169.13
## AIC=346.27   AICc=346.49   BIC=359.19

Βρίσκουμε, πως από τη συνάρτηση Arima, προτείνεται το μοντέλο

(Consumption_t - 0.746) = 0.272(Consumption_t-1 - 0.746) + 0.217(Consumption_t-2 - 0.746)

Εύκολα μπορούμε να επαληθεύσουμε ότι η δεύτερη εξίσωση είναι η ίδια με την πρώτη.Οι μικρές διαφορές οφείλονται στο γεγονός πως η lm και η ar.ols χρησιμοποιεί ελάχιστα τετράγωνα ενώ η Arima χρησιμοποιεί μέγιστη πιθανοφάνεια και φίλτρα Kalman.

6 Βιβλιογραφία

Μία εξαιρετική και πλήρης παρουσίαση πολλών μεθόδων στις χρονοσειρές είναι το

Hyndman, R.J., & Athanasopoulos, G. (2021) Forecasting: principles and practice, 3rd edition (Greek translation), OTexts: Melbourne, Australia. OTexts.com/fppgr. Accessed on 22/04/2024.

Το πιο ωραιο του χαρακτηριστικό: είναι ελεύθερα διαθεσιμο και μεταφρασμένο στην Ελληνική γλώσσα!

Ιστοσελίδα αναφοράς: https://otexts.com/fppgr/.

Επικοινωνία: Επαμεινώνδας Διαμαντοπουλος, Τμήμα ΗΜ/ΜΥ, Δ.Π.Θ. Email: epdiaman@ee.duth.gr. Τηλ: 25410 79757, 6944683327↩︎

Στοχαστική Ανάλυση Χρονοσειρών

Γραμμική Παλινδρόμηση και Χρονοσειρές

Επαμεινώνδας Διαμαντόπουλος¹

1 Ευθεία

2 Απλή Γραμμική Παλινδρόμηση (Simple Linear Regression - SLR).

2.1 Παρατηρήσεις

2.1.1 Τυποποιημένοι Συντελεστές (Standardized Coefficients)

2.1.2 Γραμμικό μοντέλο χωρίς σταθερό όρο

2.1.3 Η προϋπόθεση της κανονικότητας

2.1.4 Standardized και Studentized Residual Value

3 Απλή Γραμμική Παλινδρόμηση για εκτίμηση χρονοσειράς

3.1 Παράδειγμα Ι

3.2 Παράδειγμα ΙΙ

3.2.1 Προσαρμογή ευθείας ελαχίστων τετραγώνων

3.2.2 Προσαρμογή ευθείας με MLE

4 Πολλαπλή Γραμμική Παλινδρόμηση (Multiple Linear Regression - MLR).

4.1 Προϋποθέσεις εγκυρότητας

4.1.1 Διάγραμμα Normal Q – Q

4.1.2 Διάγραμμα Residuals vs Fitted

4.1.3 Διάγραμμα Residuals vs Leverage

4.1.4 Διάγραμμα Scale-Location

4.2 Μέθοδος Stepwise σταδιακής εισαγωγής – εξαγωγής

5 Πολλαπλή γραμμική παλινδρόμηση για την πρόβλεψη χρονοσειρών

5.1 Περίπτωση γραμμικού μοντέλου με τάση και εποχικότητα

5.1.1 Εντοπισμός τάσης

5.1.2 Εντοπισμός εποχικότητας

5.1.2.1 Ποιος είναι ο λόγος που λείπει ο Ιανουάριος από το μοντέλο;

5.1.3 Εντοπισμός τάσης και εποχικότητας με την βιβλιοθήκη tslm

5.1.4 Παράδειγμα

5.2 Περίπτωση AR(p)

5.2.1 Εφαρμογή lm και my_lm

5.2.2 Εφαρμογή ar.ols

5.2.3 Εφαρμογή ARIMA

6 Βιβλιογραφία

Στοχαστική Ανάλυση Χρονοσειρών

Γραμμική Παλινδρόμηση και Χρονοσειρές

Επαμεινώνδας Διαμαντόπουλος1

1 Ευθεία

2 Απλή Γραμμική Παλινδρόμηση (Simple Linear Regression - SLR).

2.1 Παρατηρήσεις

2.1.1 Τυποποιημένοι Συντελεστές (Standardized Coefficients)

2.1.2 Γραμμικό μοντέλο χωρίς σταθερό όρο

2.1.3 Η προϋπόθεση της κανονικότητας

2.1.4 Standardized και Studentized Residual Value

3 Απλή Γραμμική Παλινδρόμηση για εκτίμηση χρονοσειράς

3.1 Παράδειγμα Ι

3.2 Παράδειγμα ΙΙ

3.2.1 Προσαρμογή ευθείας ελαχίστων τετραγώνων

3.2.2 Προσαρμογή ευθείας με MLE

4 Πολλαπλή Γραμμική Παλινδρόμηση (Multiple Linear Regression - MLR).

4.1 Προϋποθέσεις εγκυρότητας

4.1.1 Διάγραμμα Normal Q – Q

4.1.2 Διάγραμμα Residuals vs Fitted

4.1.3 Διάγραμμα Residuals vs Leverage

4.1.4 Διάγραμμα Scale-Location

4.2 Μέθοδος Stepwise σταδιακής εισαγωγής – εξαγωγής

5 Πολλαπλή γραμμική παλινδρόμηση για την πρόβλεψη χρονοσειρών

5.1 Περίπτωση γραμμικού μοντέλου με τάση και εποχικότητα

5.1.1 Εντοπισμός τάσης

5.1.2 Εντοπισμός εποχικότητας

5.1.2.1 Ποιος είναι ο λόγος που λείπει ο Ιανουάριος από το μοντέλο;

5.1.3 Εντοπισμός τάσης και εποχικότητας με την βιβλιοθήκη tslm

5.1.4 Παράδειγμα

5.2 Περίπτωση AR(p)

5.2.1 Εφαρμογή lm και my_lm

5.2.2 Εφαρμογή ar.ols

5.2.3 Εφαρμογή ARIMA

6 Βιβλιογραφία

Επαμεινώνδας Διαμαντόπουλος¹