Australia Weather Forecast

Introduction

Η παρούσα εργασία εκπονήθηκε στο πλαίσιο του μαθήματος Επιχειρηματικής Αναλυτικής, το οποίο διδάσκεται στο 8ο εξάμηνο του Τμήματος Εφαρμοσμένης Πληροφορικής του Πανεπιστημίου Μακεδονίας. Στόχος της εργασίας είναι η εφαρμογή στατιστικής ανάλυσης και η αξιοποίηση της μεθόδου της λογιστικής παλινδρόμησης σε ένα επιλεγμένο σύνολο δεδομένων.

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

Τα δεδομένα προέρχονται από την ιστοσελίδα της Αυστραλιανής Μετεωρολογικής Υπηρεσίας (Bureau of Meteorology) και καλύπτουν περίπου δέκα έτη μετεωρολογικών μετρήσεων. Κάθε εγγραφή αντιπροσωπεύει τις καιρικές συνθήκες μιας συγκεκριμένης ημέρας και τοπικής περιοχής.

Ο στόχος του συνόλου δεδομένων είναι να προβλέψει αν θα βρέξει την επόμενη ημέρα. Η μεταβλητή-στόχος (RainTomorrow) είναι δυαδική (Ναι/Όχι) και υποδεικνύει αν η βροχόπτωση της επόμενης ημέρας θα είναι 1 χιλιοστό ή περισσότερο.

Πεδία που περιλαμβάνονται:

Ημερομηνία παρατήρησης
Θερμοκρασίες (ελάχιστη, μέγιστη)
Ποσότητα βροχής
Διεύθυνση και ταχύτητα ανέμου
Επίπεδα υγρασίας και ατμοσφαιρικής πίεσης
Ημερήσια πρόβλεψη βροχής (RainToday)
Μεταβλητή πρόβλεψης (RainTomorrow)

Η πηγή και οι ορισμοί έχουν βασιστεί στις ανακοινώσεις και στα δεδομένα που παρέχει επίσημα η Αυστραλιανή Μετεωρολογική Υπηρεσία μέσω της πλατφόρμας Climate Data Online.

Πνευματικά δικαιώματα: Commonwealth of Australia 2010, Bureau of Meteorology.

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

Data Overview

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

library(readr)
library(knitr)

df <- read.csv(file = "../weatherAUS.csv")
kable(df[1:5, ], caption = "Οπτικοποίηση αρχής δεδομένων")

Οπτικοποίηση αρχής δεδομένων
Date	Location	MinTemp	MaxTemp	Rainfall	Evaporation	Sunshine	WindGustDir	WindGustSpeed	WindDir9am	WindDir3pm	WindSpeed9am	WindSpeed3pm	Humidity9am	Humidity3pm	Pressure9am	Pressure3pm	Cloud9am	Cloud3pm	Temp9am	Temp3pm	RainToday	RainTomorrow
2008-12-01	Albury	13.4	22.9	0.6	NA	NA	W	44	W	WNW	20	24	71	22	1007.7	1007.1	8	NA	16.9	21.8	No	No
2008-12-02	Albury	7.4	25.1	0.0	NA	NA	WNW	44	NNW	WSW	4	22	44	25	1010.6	1007.8	NA	NA	17.2	24.3	No	No
2008-12-03	Albury	12.9	25.7	0.0	NA	NA	WSW	46	W	WSW	19	26	38	30	1007.6	1008.7	NA	2	21.0	23.2	No	No
2008-12-04	Albury	9.2	28.0	0.0	NA	NA	NE	24	SE	E	11	9	45	16	1017.6	1012.8	NA	NA	18.1	26.5	No	No
2008-12-05	Albury	17.5	32.3	1.0	NA	NA	W	41	ENE	NW	7	20	82	33	1010.8	1006.0	7	8	17.8	29.7	No	No

Στήλη	Περιγραφή
`Date`	Ημερομηνία παρατήρησης
`Location`	Τοποθεσία/Όνομα μετεωρολογικού σταθμού
`MinTemp`	Ελάχιστη θερμοκρασία της ημέρας (°C)
`MaxTemp`	Μέγιστη θερμοκρασία της ημέρας (°C)
`Rainfall`	Ποσότητα βροχής που καταγράφηκε (mm)
`Evaporation`	Ποσότητα εξάτμισης (mm) – μέτρηση με Class A Pan
`Sunshine`	Ώρες ηλιοφάνειας κατά τη διάρκεια της ημέρας
`WindGustDir`	Διεύθυνση του ισχυρότερου ριπιδίου ανέμου
`WindGustSpeed`	Ταχύτητα του ισχυρότερου ριπιδίου ανέμου (km/h)
`WindDir9am`	Διεύθυνση ανέμου στις 9 π.μ.
`WindDir3pm`	Διεύθυνση ανέμου στις 3 μ.μ.
`WindSpeed9am`	Ταχύτητα ανέμου στις 9 π.μ. (km/h)
`WindSpeed3pm`	Ταχύτητα ανέμου στις 3 μ.μ. (km/h)
`Humidity9am`	Υγρασία (%) στις 9 π.μ.
`Humidity3pm`	Υγρασία (%) στις 3 μ.μ.
`Pressure9am`	Ατμοσφαιρική πίεση (hPa) στις 9 π.μ.
`Pressure3pm`	Ατμοσφαιρική πίεση (hPa) στις 3 μ.μ.
`Cloud9am`	Κάλυψη νεφών στις 9 π.μ. (κλίμακα 0–8)
`Cloud3pm`	Κάλυψη νεφών στις 3 μ.μ. (κλίμακα 0–8)
`Temp9am`	Θερμοκρασία στις 9 π.μ. (°C)
`Temp3pm`	Θερμοκρασία στις 3 μ.μ. (°C)
`RainToday`	Έβρεξε σήμερα; Ναι/Όχι (>=1mm)
`RainTomorrow`	Θα βρέξει αύριο; Ναι/Όχι (>=1mm) – μεταβλητή στόχος

Data Preprocessing

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

summary(df)

##      Date             Location            MinTemp         MaxTemp     
##  Length:145460      Length:145460      Min.   :-8.50   Min.   :-4.80  
##  Class :character   Class :character   1st Qu.: 7.60   1st Qu.:17.90  
##  Mode  :character   Mode  :character   Median :12.00   Median :22.60  
##                                        Mean   :12.19   Mean   :23.22  
##                                        3rd Qu.:16.90   3rd Qu.:28.20  
##                                        Max.   :33.90   Max.   :48.10  
##                                        NA's   :1485    NA's   :1261   
##     Rainfall        Evaporation        Sunshine     WindGustDir       
##  Min.   :  0.000   Min.   :  0.00   Min.   : 0.00   Length:145460     
##  1st Qu.:  0.000   1st Qu.:  2.60   1st Qu.: 4.80   Class :character  
##  Median :  0.000   Median :  4.80   Median : 8.40   Mode  :character  
##  Mean   :  2.361   Mean   :  5.47   Mean   : 7.61                     
##  3rd Qu.:  0.800   3rd Qu.:  7.40   3rd Qu.:10.60                     
##  Max.   :371.000   Max.   :145.00   Max.   :14.50                     
##  NA's   :3261      NA's   :62790    NA's   :69835                     
##  WindGustSpeed     WindDir9am         WindDir3pm         WindSpeed9am   
##  Min.   :  6.00   Length:145460      Length:145460      Min.   :  0.00  
##  1st Qu.: 31.00   Class :character   Class :character   1st Qu.:  7.00  
##  Median : 39.00   Mode  :character   Mode  :character   Median : 13.00  
##  Mean   : 40.03                                         Mean   : 14.04  
##  3rd Qu.: 48.00                                         3rd Qu.: 19.00  
##  Max.   :135.00                                         Max.   :130.00  
##  NA's   :10263                                          NA's   :1767    
##   WindSpeed3pm    Humidity9am      Humidity3pm      Pressure9am    
##  Min.   : 0.00   Min.   :  0.00   Min.   :  0.00   Min.   : 980.5  
##  1st Qu.:13.00   1st Qu.: 57.00   1st Qu.: 37.00   1st Qu.:1012.9  
##  Median :19.00   Median : 70.00   Median : 52.00   Median :1017.6  
##  Mean   :18.66   Mean   : 68.88   Mean   : 51.54   Mean   :1017.6  
##  3rd Qu.:24.00   3rd Qu.: 83.00   3rd Qu.: 66.00   3rd Qu.:1022.4  
##  Max.   :87.00   Max.   :100.00   Max.   :100.00   Max.   :1041.0  
##  NA's   :3062    NA's   :2654     NA's   :4507     NA's   :15065   
##   Pressure3pm        Cloud9am        Cloud3pm        Temp9am     
##  Min.   : 977.1   Min.   :0.00    Min.   :0.00    Min.   :-7.20  
##  1st Qu.:1010.4   1st Qu.:1.00    1st Qu.:2.00    1st Qu.:12.30  
##  Median :1015.2   Median :5.00    Median :5.00    Median :16.70  
##  Mean   :1015.3   Mean   :4.45    Mean   :4.51    Mean   :16.99  
##  3rd Qu.:1020.0   3rd Qu.:7.00    3rd Qu.:7.00    3rd Qu.:21.60  
##  Max.   :1039.6   Max.   :9.00    Max.   :9.00    Max.   :40.20  
##  NA's   :15028    NA's   :55888   NA's   :59358   NA's   :1767   
##     Temp3pm       RainToday         RainTomorrow      
##  Min.   :-5.40   Length:145460      Length:145460     
##  1st Qu.:16.60   Class :character   Class :character  
##  Median :21.10   Mode  :character   Mode  :character  
##  Mean   :21.68                                        
##  3rd Qu.:26.40                                        
##  Max.   :46.70                                        
##  NA's   :3609

Empty Values Check

missing_val <- sapply(df, function(x) {
  paste(
    round(sum(is.na(x)) / length(x) * 100, digits = 2),
    "%"
  )
})

df_missing <- data.frame(
  Missing_Values = missing_val
)

kable(df_missing)

	Missing_Values
Date	0 %
Location	0 %
MinTemp	1.02 %
MaxTemp	0.87 %
Rainfall	2.24 %
Evaporation	43.17 %
Sunshine	48.01 %
WindGustDir	7.1 %
WindGustSpeed	7.06 %
WindDir9am	7.26 %
WindDir3pm	2.91 %
WindSpeed9am	1.21 %
WindSpeed3pm	2.11 %
Humidity9am	1.82 %
Humidity3pm	3.1 %
Pressure9am	10.36 %
Pressure3pm	10.33 %
Cloud9am	38.42 %
Cloud3pm	40.81 %
Temp9am	1.21 %
Temp3pm	2.48 %
RainToday	2.24 %
RainTomorrow	2.25 %

Όπως εμφανίζεται στον παραπάνω πίνακα, υπάρχουν πεδία που περιλαμβάνουν κενές τιμές, αυτό προκαλεί πρόβλημα κατά την εκπαίδευση του μοντέλου οπότε θα πρέπει να αφαίρεθούνε

df_clean <- na.omit(df)

print(paste(
  "Before:", nrow(df),
  "After omit:", nrow(df_clean),
  "Diff:", nrow(df) - nrow(df_clean), "removed"
))

## [1] "Before: 145460 After omit: 56420 Diff: 89040 removed"

Correlation Matrix

library(corrplot)

## corrplot 0.95 loaded

# Calculate correlation matrix
cor_matrix <- cor(
  df_clean[, sapply(df_clean, is.numeric)],
  method = "pearson"
)

# Visualize
corrplot(cor_matrix,
  method = "color",
  type = "upper",
  tl.col = "black", # text label color
  tl.srt = 45, # text label rotation
  tl.cex = 0.8, # text label size (default is 1)
  addCoef.col = "black", # show correlation coefficients
  number.cex = 0.7 # coefficient number size
)

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

Σε αυτήν την περίπτωση μπορούμε να αφαιρέσουμε τα ακόλουθα:

Temp3pm - MaxTemp
Temp9am - MinTemp
Pressure9am - Pressure3pm

Data Visualization

library(ggplot2)

ggplot(df_clean, aes(x = factor(Cloud9am), y = Sunshine)) +
  geom_boxplot() +
  labs(x = "Clouds 9am", y = "Sunshine", title = "Cloudy - Sunshine boxplot")

Το διάγραμμα παρουσιάζει τη σχέση μεταξύ της νεφοκάλυψης στις 9 π.μ. (Clouds 9am) και της διάρκειας ηλιοφάνειας (Sunshine) κατά τη διάρκεια της ημέρας. Χρησιμοποιείται διάγραμμα τύπου boxplot για να δείξει την κατανομή των τιμών ηλιοφάνειας για κάθε επίπεδο νεφοκάλυψης (από 0 έως 8).

Ο άξονας Χ (Clouds 9am): δείχνει την τιμή νεφοκάλυψης, με την κλίμακα από 0 (καθόλου σύννεφα) έως 8 (πλήρης νεφοκάλυψη). Ο άξονας Υ (Sunshine): παρουσιάζει τις ώρες ηλιοφάνειας. Με συμπέρασμα:

Όσο αυξάνεται η νεφοκάλυψη, η ηλιοφάνεια μειώνεται σταθερά.

Για τιμές νεφοκάλυψης 0 έως 2, η διάμεσος ηλιοφάνειας είναι υψηλή (~10 ώρες).

Αντίθετα, για τιμές 7 και 8, η διάμεσος πέφτει σε πολύ χαμηλά επίπεδα (~3-5 ώρες).

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

ggplot(df_clean, aes(x = Humidity3pm, y = Humidity9am)) +
  geom_point() +
  labs(x = "Humidity 3pm", y = "Humidity 9am", title = "Humidity at 3pm and 9am")

Το παραπάνω διάγραμμα διασποράς (scatter plot) απεικονίζει τη σχέση μεταξύ της υγρασίας στις 3 μ.μ. (άξονας Χ) και της υγρασίας στις 9 π.μ. (άξονας Υ) για πλήθος ημερήσιων μετεωρολογικών παρατηρήσεων. Κάθε κουκίδα στο γράφημα αντιπροσωπεύει μία ημέρα μέτρησης σε έναν συγκεκριμένο σταθμό.

Παρατηρούμε ότι υπάρχει μια θετική συσχέτιση μεταξύ των δύο μεταβλητών: γενικά, όταν η υγρασία είναι υψηλή στις 3 μ.μ., είναι επίσης υψηλή και στις 9 π.μ. Αντίστοιχα, όταν η υγρασία είναι χαμηλή το πρωί, τείνει να είναι χαμηλή και το απόγευμα. Ωστόσο, η κατανομή των σημείων δεν είναι απολύτως γραμμική – παρατηρείται μια “τριγωνική” κατανομή των δεδομένων, με πολλές τιμές υγρασίας στις 9 π.μ. να βρίσκονται κοντά στο μέγιστο (100%), ακόμη και όταν η απογευματινή υγρασία είναι μέτρια ή χαμηλή.

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

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

ggplot(df_clean, aes(x = WindGustSpeed, y = Rainfall)) +
  geom_density2d() +
  labs(x = "Wind Gust Speed", y = "Rainfall", title = "geom_density2d")

Το παραπάνω διάγραμμα παρουσιάζει ένα δισδιάστατο διάγραμμα πυκνότητας (2D density plot) μεταξύ των μεταβλητών ένταση ριπής ανέμου (Wind Gust Speed) και ύψος βροχόπτωσης (Rainfall). Η συγκεκριμένη απεικόνιση, μέσω των ισοϋψών καμπυλών (contour lines), αποτυπώνει τις περιοχές στο επίπεδο των δύο μεταβλητών όπου συγκεντρώνεται μεγαλύτερη πυκνότητα παρατηρήσεων.

Αναλύοντας το γράφημα, διαπιστώνουμε ότι η μεγαλύτερη συγκέντρωση παρατηρήσεων εντοπίζεται γύρω από τις μεσαίες τιμές της έντασης των ριπών ανέμου (περίπου 30–40 km/h) και χαμηλά επίπεδα βροχόπτωσης (0 έως 1 mm). Οι καμπύλες πυκνότητας αποκαλύπτουν ότι οι περισσότερες ημερήσιες παρατηρήσεις αφορούν περιπτώσεις όπου οι ριπές ανέμου δεν είναι ιδιαίτερα έντονες, και η βροχόπτωση παραμένει περιορισμένη. Αυτό είναι σύμφωνο με τις συνήθεις καιρικές συνθήκες σε πολλές περιοχές όπου οι ακραίες καιρικές καταστάσεις, όπως ισχυρές καταιγίδες, είναι λιγότερο συχνές.

Καθώς αυξάνεται η ένταση του ανέμου πέραν των 40–50 km/h, η πυκνότητα των σημείων μειώνεται, γεγονός που υποδηλώνει ότι αυτές οι τιμές είναι πιο σπάνιες. Παρόμοια παρατηρείται και για τις υψηλές τιμές βροχόπτωσης (>2 mm), οι οποίες εμφανίζονται σποραδικά και δεν συγκεντρώνουν μεγάλο όγκο παρατηρήσεων.

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

ggplot(df_clean, aes(x = Evaporation, y = Temp3pm)) +
  geom_point() +
  labs(x = "Evaporation", y = "Temp 3pm", title = "Evaporation - Temp")

Το διάγραμμα διασποράς απεικονίζει τη σχέση μεταξύ της εξάτμισης (“Evaporation”) και της θερμοκρασίας στις 3 μ.μ. (Temp 3pm). Παρατηρώντας το σύνολο των σημείων, διαπιστώνεται μια θετική συσχέτιση μεταξύ των δύο μεταβλητών. Αυτό σημαίνει ότι, γενικά, καθώς αυξάνεται η εξάτμιση, τείνει να αυξάνεται και η θερμοκρασία το απόγευμα.

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

ggplot(df_clean, aes(x = Location, fill = RainTomorrow)) +
  geom_bar(position = "fill") +
  coord_flip() +
  labs(title = "Πιθανότητα Βροχής ανά Τοποθεσία", y = "Ποσοστό", x = "Περιοχή")

Αυτό το οριζόντιο ραβδόγραμμα απεικονίζει την πιθανότητα βροχής την επόμενη ημέρα για μια ποικιλία τοποθεσιών. Στον κάθετο άξονα παρατίθενται οι διάφορες γεωγραφικές περιοχές (“Περιοχή”), ενώ ο οριζόντιος άξονας αντιπροσωπεύει το ποσοστό ή την πιθανότητα (“Ποσοστό”), κυμαινόμενο από 0 έως 1 (ή 0% έως 100%). Κάθε ράβδος για κάθε τοποθεσία είναι διαχωρισμένη σε δύο χρώματα, που αντιπροσωπεύουν την πρόβλεψη για τη βροχή αύριο (“RainTomorrow”): το γαλάζιο τμήμα αντιστοιχεί στην πιθανότητα “Ναι” (θα βρέξει), ενώ το κόκκινο τμήμα στην πιθανότητα “Όχι” (δεν θα βρέξει).

Εξετάζοντας τις ράβδους, μπορούμε εύκολα να συγκρίνουμε την πιθανότητα βροχής μεταξύ των διαφόρων τοποθεσιών. Ορισμένες περιοχές, όπως η Cairns και η Darwin, εμφανίζουν σημαντικά μεγαλύτερο γαλάζιο τμήμα, υποδεικνύοντας υψηλότερη πιθανότητα βροχής την επόμενη ημέρα σε σύγκριση με άλλες. Αντίθετα, περιοχές όπως η Woomera και η Williamtown έχουν πολύ μικρότερο γαλάζιο τμήμα και κυριαρχείται η κόκκινη περιοχή, υποδηλώνοντας χαμηλή πιθανότητα βροχής. Το διάγραμμα καθιστά σαφές ότι η πιθανότητα βροχής την επόμενη ημέρα ποικίλλει σημαντικά από τη μία τοποθεσία στην άλλη, παρέχοντας μια οπτική σύνοψη των προβλέψεων για κάθε καταγεγραμμένη περιοχή.

ggplot(df_clean, aes(x = RainTomorrow, y = Temp3pm, fill = RainTomorrow)) +
  geom_boxplot() +
  labs(title = "Θερμοκρασία στις 3μμ σε σχέση με τη Βροχή")

ggplot(df_clean, aes(x = Temp3pm, y = Humidity3pm)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm") +
  labs(title = "Θερμοκρασία vs Υγρασία στις 3μμ")

df_city <- df_clean[df_clean$Location == "Williamtown", ]

ggplot(df_city, aes(x = Temp3pm, y = Humidity3pm)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm") +
  labs(title = "Θερμοκρασία vs Υγρασία στις 3μμ στο William")

Αυτό το διάγραμμα διασποράς φέρει τον τίτλο “Θερμοκρασία vs Υγρασία στις 3μμ” και αναπαριστά τη σχέση μεταξύ της θερμοκρασίας και της σχετικής υγρασίας που καταγράφηκαν στις 3 το απόγευμα. Ο οριζόντιος άξονας αναπαριστά τη θερμοκρασία στις 3 μ.μ. (“Temp3pm”), ενώ ο κάθετος άξονας αναπαριστά την υγρασία στις 3 μ.μ. (“Humidity3pm”). Τα σημεία στο διάγραμμα αντιπροσωπεύουν μεμονωμένες μετρήσεις των δύο αυτών μεταβλητών.

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

Values Transformation

df_city$RainTomorrow <- ifelse(df_city$RainTomorrow == "Yes", 1, 0)
df_city$RainToday <- ifelse(df_city$RainToday == "Yes", 1, 0)

df_city$WindDir3pm <- as.factor(df_city$WindDir3pm)
df_city$WindDir9am <- as.factor(df_city$WindDir9am)
df_city$WindGustDir <- as.factor(df_city$WindGustDir)
df_city$Location <- as.factor(df_city$Location)

library(lubridate)

## Warning: package 'lubridate' was built under R version 4.4.3

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

df_city$Date <- as.Date(df_city$Date)
df_city$Date <-
  year(df_city$Date) + (yday(df_city$Date) - 1) /
    ifelse(leap_year(df_city$Date), 366, 365)

kable(df_city[5:10, ])

	Date	Location	MinTemp	MaxTemp	Rainfall	Evaporation	Sunshine	WindGustDir	WindGustSpeed	WindDir9am	WindDir3pm	WindSpeed9am	WindSpeed3pm	Humidity9am	Humidity3pm	Pressure9am	Pressure3pm	Cloud9am	Cloud3pm	Temp9am	Temp3pm	RainToday	RainTomorrow
39543	2009.011	Williamtown	14.8	36.0	0.0	8.6	11.9	ENE	37	NNW	ENE	15	26	67	39	1013.6	1008.7	2	3	22.5	33.8	0	0
39544	2009.014	Williamtown	18.6	33.7	0.0	9.2	13.0	ENE	39	SSE	ESE	13	19	68	49	1011.1	1008.1	1	1	26.3	29.2	0	0
39545	2009.016	Williamtown	19.7	35.5	0.0	11.0	5.3	NE	30	S	ENE	15	19	81	40	1009.9	1005.8	4	7	24.1	33.8	0	0
39546	2009.019	Williamtown	20.7	23.8	0.0	7.8	0.1	SSW	48	SSW	S	35	24	83	79	1010.3	1011.9	7	8	21.7	21.8	0	1
39547	2009.022	Williamtown	17.6	23.0	1.4	8.6	1.9	S	43	SE	SSE	24	28	61	51	1018.7	1017.6	8	6	19.1	21.9	1	0
39548	2009.025	Williamtown	13.3	24.9	0.0	4.2	10.1	ESE	44	NE	ESE	7	30	57	46	1016.2	1012.8	6	6	20.6	22.9	0	0

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

Ένα από τα βασικά βήματα είναι η μετατροπή των δυαδικών μεταβλητών, όπως οι RainTomorrow και RainToday, από τις αρχικές κατηγορικές τιμές (“Yes”, “No”) σε αριθμητική δυαδική μορφή (1, 0). Η λογιστική παλινδρόμηση μοντελοποιεί την πιθανότητα να συμβεί ένα συγκεκριμένο γεγονός, το οποίο στην περίπτωση αυτή είναι η βροχή. Η μαθηματική διατύπωση του μοντέλου βασίζεται σε μια δυαδική εξαρτημένη μεταβλητή, όπου η μία τιμή (συνήθως το 1) αντιπροσωπεύει την “επιτυχία” ή την εκδήλωση του γεγονότος που μας ενδιαφέρει να προβλέψουμε (π.χ., βροχή), και η άλλη τιμή (το 0) αντιπροσωπεύει την “αποτυχία” ή τη μη εκδήλωση του γεγονότος. Επομένως, η κωδικοποίηση των μεταβλητών πρόβλεψης βροχής ως 0 και 1 είναι απαραίτητη προϋπόθεση για την εκτίμηση των παραμέτρων του μοντέλου και την ερμηνεία των αποτελεσμάτων ως πιθανότητες.

Επιπλέον, μεταβλητές που αντιπροσωπεύουν κατηγορίες χωρίς εγγενή διάταξη, όπως οι διευθύνσεις ανέμου (WindDir3pm, WindDir9am, WindGustDir) και η Location, μετατράπηκαν σε παράγοντες (factor). Αυτή η μετατροπή είναι κρίσιμη διότι γνωστοποιεί στο στατιστικό λογισμικό ότι οι τιμές αυτές δεν είναι αριθμητικές με ποσοτική σημασία, αλλά διακριτές κατηγορίες. Όταν ένας παράγοντας εισάγεται σε ένα μοντέλο παλινδρόμησης στην R, το λογισμικό δημιουργεί αυτόματα ένα σύνολο από δυαδικές (dummy) μεταβλητές. Κάθε μία από αυτές τις dummy μεταβλητές αντιστοιχεί σε μία από τις κατηγορίες του παράγοντα (με μία κατηγορία να ορίζεται ως αναφορά) και λαμβάνει τιμή 1 εάν η παρατήρηση ανήκει στην αντίστοιχη κατηγορία και 0 διαφορετικά. Με αυτόν τον τρόπο, το μοντέλο λογιστικής παλινδρόμησης μπορεί να εκτιμήσει την ξεχωριστή επίδραση κάθε κατηγορίας (π.χ., την επίδραση της βόρειας κατεύθυνσης ανέμου σε σύγκριση με την κατεύθυνση αναφοράς) στην πιθανότητα του υπό εξέταση γεγονότος (π.χ., βροχή).

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

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

Train Test Split

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

library(caTools)
set.seed(183)

split <- sample.split(df_city$RainTomorrow, SplitRatio = 0.65)

train <- subset(df_city, split == TRUE)
test <- subset(df_city, split == FALSE)

print(paste(
  "Αναλογία Train / Test:", nrow(train), "/", nrow(test), "=",
  round(nrow(train) / nrow(test), 3)
))

## [1] "Αναλογία Train / Test: 779 / 419 = 1.859"

Πριν από κάθε εκπαίδευση και αξιολόγηση του μοντέλου, το σύνολο δεδομένων διαχωρίζεται σε δύο κατηγορίες. Η πρώτη κατηγορία ονομάζεται Train set και χρησιμοποιείται για την εκπαίδευση του μοντέλου. Επιπλέον, σε ορισμένες περιπτώσεις, το Train set μπορεί να χρησιμοποιηθεί και ως validation set όταν δουλεύουμε με πιο σύνθετα μοντέλα, όπως τα νευρωνικά δίκτυα, για να αξιολογήσουμε την απόδοση του μοντέλου κατά τη διάρκεια της εκπαίδευσης.

Η δεύτερη κατηγορία δεδομένων είναι το Test set, το οποίο χρησιμεύει για την τελική αξιολόγηση του μοντέλου αφού έχει ολοκληρωθεί η εκπαίδευση. Το Test set πρέπει να παραμείνει “αόρατο” από το μοντέλο κατά τη διάρκεια της εκπαίδευσης, και καμία προεπεξεργασία (όπως κανονικοποίηση με τη μέθοδο min-max) δεν πρέπει να εφαρμοστεί σε αυτό το σύνολο. Αν γίνει οποιαδήποτε επεξεργασία ή η χρήση των δεδομένων του Test set πριν από την τελική αξιολόγηση, αυτό μπορεί να οδηγήσει σε data leakage, δηλαδή την αποκατάσταση κρίσιμων πληροφοριών που επηρεάζουν την εκπαίδευση του μοντέλου, καθιστώντας την αξιολόγηση αναξιόπιστη. Ο σωστός διαχωρισμός και η προστασία του Test set είναι κρίσιμος για τη σωστή και αξιόπιστη εκτίμηση της απόδοσης του μοντέλου.

Logistic Regression

Fit Train Function

show_scores <- function(preds, df){
# ROC AUC function
  rocr_preds <- prediction(preds, df$RainTomorrow)

  print("Test Results:")
  table(df$RainTomorrow, preds > 0.5)

  predicted_binary <- preds > 0.5
  TP <- sum(df$RainTomorrow == 1 & predicted_binary)
  TN <- sum(df$RainTomorrow == 0 & !predicted_binary)

  FP <- sum(df$RainTomorrow == 0 & predicted_binary)
  FN <- sum(df$RainTomorrow == 1 & !predicted_binary)

  rocr_perd <- performance(rocr_preds, "tpr", "fpr")

  plot(rocr_perd,
    colorize = TRUE,
    print.cutoffs.at = seq(0, 1, 0.1),
    text.adj = c(-0.2, 1.7),
    main = "ROC curve"
  )

  score <- data.frame(
    AUC = as.numeric(performance(rocr_preds, "auc")@y.values),
    Accuracy = (TP + TN) / (TP + TN + FP + FN),
    Precision = TP / (TP + FP),
    Recall = TP / (TP + FN),
    F1 = (2 * TP) / (2 * TP + FP + FN),
    MCC = (TP * TN - FP * FN) / sqrt((TP + FP) * (TP + FN) * (TN + FP) * (TN + FN))
  )

  kable(score)
}

library(ROCR)

fit_and_predict <- function(cols, train, test) {
  print("Training model...")
  model <- glm(as.formula(colnames), data = train, family = "binomial")
  print(summary(model))

  print("Predicting Test set...")
  preds <- predict(model, type = "response", newdata = test)
  rocr_preds <- prediction(preds, test$RainTomorrow)

  show_scores(preds, test)
}

Characteristics

## [1] "Date + Location + MinTemp + MaxTemp + Rainfall + Evaporation + Sunshine + WindGustDir + WindGustSpeed + WindDir9am + WindDir3pm + WindSpeed9am + WindSpeed3pm + Humidity9am + Humidity3pm + Pressure9am + Pressure3pm + Cloud9am + Cloud3pm + Temp9am + Temp3pm + RainToday + RainTomorrow"

Fit Test 1

colnames <- "RainTomorrow ~ Date + MinTemp + MaxTemp + Rainfall + Evaporation + Sunshine + WindGustDir + WindGustSpeed + WindDir9am + WindDir3pm + WindSpeed9am + WindSpeed3pm + Humidity9am + Humidity3pm + Pressure9am + Pressure3pm + Cloud9am + Cloud3pm + Temp9am + Temp3pm + RainToday"
fit_and_predict(colnames, train, test)

## [1] "Training model..."
## 
## Call:
## glm(formula = as.formula(colnames), family = "binomial", data = train)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.280e+02  1.051e+02  -1.219  0.22299    
## Date            7.483e-02  5.222e-02   1.433  0.15190    
## MinTemp        -9.757e-02  6.833e-02  -1.428  0.15331    
## MaxTemp         2.200e-01  1.031e-01   2.135  0.03280 *  
## Rainfall        3.879e-02  1.941e-02   1.999  0.04565 *  
## Evaporation     2.056e-04  3.724e-02   0.006  0.99559    
## Sunshine       -2.864e-01  6.643e-02  -4.310 1.63e-05 ***
## WindGustDirENE  3.632e-01  6.476e-01   0.561  0.57488    
## WindGustDirESE -4.493e-01  6.396e-01  -0.702  0.48243    
## WindGustDirN   -1.069e-01  1.291e+00  -0.083  0.93400    
## WindGustDirNE  -4.589e-01  7.907e-01  -0.580  0.56167    
## WindGustDirNNE  2.021e-01  9.867e-01   0.205  0.83772    
## WindGustDirNNW -2.766e-01  1.170e+00  -0.236  0.81313    
## WindGustDirNW  -9.184e-01  7.956e-01  -1.154  0.24836    
## WindGustDirS   -4.054e-01  6.473e-01  -0.626  0.53107    
## WindGustDirSE  -7.935e-01  6.267e-01  -1.266  0.20540    
## WindGustDirSSE -7.528e-01  6.660e-01  -1.130  0.25836    
## WindGustDirSSW -1.645e-01  6.681e-01  -0.246  0.80546    
## WindGustDirSW   3.937e-01  8.920e-01   0.441  0.65890    
## WindGustDirW   -5.533e-01  8.391e-01  -0.659  0.50958    
## WindGustDirWNW -6.428e-01  6.624e-01  -0.970  0.33185    
## WindGustDirWSW  8.204e-03  1.017e+00   0.008  0.99357    
## WindGustSpeed   6.338e-02  1.594e-02   3.976 7.01e-05 ***
## WindDir9amENE   1.746e-01  1.054e+00   0.166  0.86845    
## WindDir9amESE  -1.090e-01  1.023e+00  -0.107  0.91511    
## WindDir9amN     1.187e-02  1.031e+00   0.012  0.99081    
## WindDir9amNE   -5.535e-01  1.065e+00  -0.520  0.60335    
## WindDir9amNNE   8.316e-01  9.952e-01   0.836  0.40336    
## WindDir9amNNW  -2.583e-01  1.040e+00  -0.248  0.80380    
## WindDir9amNW   -3.605e-01  9.685e-01  -0.372  0.70972    
## WindDir9amS    -5.415e-01  1.009e+00  -0.537  0.59156    
## WindDir9amSE    8.854e-01  9.690e-01   0.914  0.36089    
## WindDir9amSSE  -7.463e-01  1.078e+00  -0.692  0.48889    
## WindDir9amSSW   1.946e-01  1.019e+00   0.191  0.84855    
## WindDir9amSW   -7.593e-02  1.010e+00  -0.075  0.94008    
## WindDir9amW    -6.531e-01  9.932e-01  -0.658  0.51080    
## WindDir9amWNW  -3.511e-01  9.527e-01  -0.369  0.71246    
## WindDir9amWSW   4.040e-01  9.972e-01   0.405  0.68538    
## WindDir3pmENE   4.744e-01  6.638e-01   0.715  0.47486    
## WindDir3pmESE   1.225e+00  5.983e-01   2.048  0.04061 *  
## WindDir3pmN    -6.474e-01  1.081e+00  -0.599  0.54929    
## WindDir3pmNE    6.626e-01  7.451e-01   0.889  0.37388    
## WindDir3pmNNE   9.564e-01  1.010e+00   0.947  0.34361    
## WindDir3pmNNW   2.407e-01  9.159e-01   0.263  0.79273    
## WindDir3pmNW    2.091e-02  8.504e-01   0.025  0.98038    
## WindDir3pmS     1.018e+00  6.999e-01   1.455  0.14573    
## WindDir3pmSE    1.220e+00  6.093e-01   2.002  0.04533 *  
## WindDir3pmSSE   1.092e+00  6.877e-01   1.588  0.11221    
## WindDir3pmSSW   1.846e+00  8.234e-01   2.242  0.02493 *  
## WindDir3pmSW    1.557e+00  8.866e-01   1.756  0.07908 .  
## WindDir3pmW     9.064e-02  9.230e-01   0.098  0.92178    
## WindDir3pmWNW  -2.528e-01  7.908e-01  -0.320  0.74919    
## WindDir3pmWSW   2.277e+00  1.000e+00   2.276  0.02284 *  
## WindSpeed9am    2.531e-02  2.235e-02   1.133  0.25740    
## WindSpeed3pm   -2.554e-02  2.386e-02  -1.070  0.28441    
## Humidity9am     5.817e-03  1.645e-02   0.354  0.72353    
## Humidity3pm     5.009e-02  1.663e-02   3.013  0.00259 ** 
## Pressure9am     2.031e-02  9.727e-02   0.209  0.83458    
## Pressure3pm    -4.838e-02  9.358e-02  -0.517  0.60518    
## Cloud9am       -6.200e-02  6.931e-02  -0.895  0.37101    
## Cloud3pm       -8.269e-03  7.594e-02  -0.109  0.91329    
## Temp9am        -1.609e-01  1.113e-01  -1.446  0.14817    
## Temp3pm        -8.053e-03  1.176e-01  -0.068  0.94540    
## RainToday       1.954e-02  3.335e-01   0.059  0.95329    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 904.05  on 778  degrees of freedom
## Residual deviance: 567.81  on 715  degrees of freedom
## AIC: 695.81
## 
## Number of Fisher Scoring iterations: 5
## 
## [1] "Predicting Test set..."
## [1] "Test Results:"

AUC	Accuracy	Precision	Recall	F1	MCC
0.8200035	0.7947494	0.6666667	0.4642857	0.5473684	0.431593

Αφαιρούμε τις στήλες με το υψηλότερο \(Pr(>|z|)\): WindGustDir, WindDir9am, Humidity9am, Pressure9am, Pressure3pm

Fit Test 2

colnames <- "RainTomorrow ~ Date + MinTemp + MaxTemp + Rainfall + Evaporation + Sunshine + WindGustSpeed + WindDir3pm + WindSpeed9am + WindSpeed3pm + Humidity3pm + Cloud9am + Cloud3pm + Temp9am + RainToday"

fit_and_predict(colnames, train, test)

## [1] "Training model..."
## 
## Call:
## glm(formula = as.formula(colnames), family = "binomial", data = train)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.163e+02  9.878e+01  -1.178 0.238983    
## Date           5.437e-02  4.910e-02   1.107 0.268117    
## MinTemp       -7.820e-02  5.380e-02  -1.453 0.146113    
## MaxTemp        2.183e-01  6.159e-02   3.544 0.000394 ***
## Rainfall       4.462e-02  1.862e-02   2.396 0.016556 *  
## Evaporation    4.367e-03  3.532e-02   0.124 0.901610    
## Sunshine      -2.622e-01  6.197e-02  -4.231 2.33e-05 ***
## WindGustSpeed  7.086e-02  1.417e-02   5.002 5.68e-07 ***
## WindDir3pmENE  6.549e-01  6.083e-01   1.077 0.281661    
## WindDir3pmESE  9.567e-01  5.364e-01   1.784 0.074485 .  
## WindDir3pmN   -3.993e-01  9.289e-01  -0.430 0.667278    
## WindDir3pmNE   9.003e-01  6.655e-01   1.353 0.176104    
## WindDir3pmNNE  8.572e-01  9.437e-01   0.908 0.363686    
## WindDir3pmNNW  5.257e-02  8.214e-01   0.064 0.948975    
## WindDir3pmNW  -1.435e-01  7.349e-01  -0.195 0.845142    
## WindDir3pmS    6.545e-01  5.598e-01   1.169 0.242356    
## WindDir3pmSE   8.494e-01  5.233e-01   1.623 0.104575    
## WindDir3pmSSE  4.472e-01  5.718e-01   0.782 0.434135    
## WindDir3pmSSW  1.498e+00  6.956e-01   2.154 0.031252 *  
## WindDir3pmSW   1.367e+00  7.625e-01   1.793 0.072945 .  
## WindDir3pmW   -2.819e-01  8.286e-01  -0.340 0.733652    
## WindDir3pmWNW -6.500e-01  6.637e-01  -0.979 0.327373    
## WindDir3pmWSW  2.090e+00  8.785e-01   2.378 0.017388 *  
## WindSpeed9am   1.436e-02  1.819e-02   0.789 0.429843    
## WindSpeed3pm  -2.661e-02  2.122e-02  -1.254 0.210014    
## Humidity3pm    5.303e-02  1.068e-02   4.966 6.83e-07 ***
## Cloud9am      -2.669e-02  6.359e-02  -0.420 0.674707    
## Cloud3pm      -1.420e-02  6.939e-02  -0.205 0.837843    
## Temp9am       -1.564e-01  7.891e-02  -1.982 0.047472 *  
## RainToday      6.143e-02  3.034e-01   0.202 0.839562    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 904.05  on 778  degrees of freedom
## Residual deviance: 593.30  on 749  degrees of freedom
## AIC: 653.3
## 
## Number of Fisher Scoring iterations: 5
## 
## [1] "Predicting Test set..."
## [1] "Test Results:"

AUC	Accuracy	Precision	Recall	F1	MCC
0.8378025	0.8114558	0.7088608	0.5	0.5863874	0.4809471

Παρατηρούμε πως η απόδοση αυξήθηκε σημαντικά σε όλες τις μετρηκές!

Αφαιρούμε τα χαρακτηριστικά: WindGustDir, WindDir9am, Humidity9am, Pressure9am, Pressure3pm, Temp3pm

Fit Test 3

colnames <- "RainTomorrow ~ MinTemp + MaxTemp + Rainfall + Sunshine + WindGustSpeed + WindDir3pm + Humidity3pm + Temp9am"

fit_and_predict(colnames, train, test)

## [1] "Training model..."
## 
## Call:
## glm(formula = as.formula(colnames), family = "binomial", data = train)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -7.198958   1.200590  -5.996 2.02e-09 ***
## MinTemp       -0.082063   0.047168  -1.740  0.08189 .  
## MaxTemp        0.216089   0.054149   3.991 6.59e-05 ***
## Rainfall       0.046688   0.015649   2.984  0.00285 ** 
## Sunshine      -0.243384   0.042744  -5.694 1.24e-08 ***
## WindGustSpeed  0.063493   0.010814   5.872 4.32e-09 ***
## WindDir3pmENE  0.633659   0.595367   1.064  0.28719    
## WindDir3pmESE  0.869584   0.521332   1.668  0.09532 .  
## WindDir3pmN   -0.363574   0.927964  -0.392  0.69521    
## WindDir3pmNE   0.987772   0.654498   1.509  0.13125    
## WindDir3pmNNE  0.994427   0.930722   1.068  0.28532    
## WindDir3pmNNW  0.273464   0.812700   0.336  0.73650    
## WindDir3pmNW  -0.084178   0.728621  -0.116  0.90802    
## WindDir3pmS    0.605038   0.534720   1.132  0.25784    
## WindDir3pmSE   0.837232   0.507058   1.651  0.09871 .  
## WindDir3pmSSE  0.473300   0.550237   0.860  0.38969    
## WindDir3pmSSW  1.534677   0.680600   2.255  0.02414 *  
## WindDir3pmSW   1.540699   0.754626   2.042  0.04118 *  
## WindDir3pmW   -0.089492   0.807823  -0.111  0.91179    
## WindDir3pmWNW -0.573863   0.645398  -0.889  0.37392    
## WindDir3pmWSW  2.273997   0.868617   2.618  0.00885 ** 
## Humidity3pm    0.053470   0.009993   5.351 8.76e-08 ***
## Temp9am       -0.156199   0.071743  -2.177  0.02947 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 904.05  on 778  degrees of freedom
## Residual deviance: 597.10  on 756  degrees of freedom
## AIC: 643.1
## 
## Number of Fisher Scoring iterations: 5
## 
## [1] "Predicting Test set..."
## [1] "Test Results:"

AUC	Accuracy	Precision	Recall	F1	MCC
0.8408853	0.8114558	0.7142857	0.4910714	0.5820106	0.4792463

Η απόδοση έπεσε.

Scale Fit Test 4

Έστω να δοκιμάσουμε με min-max scale.

scale_with_categoricals <- function(df) {
  numeric_cols <- sapply(df, is.numeric)
  numeric_cols["RainTomorrow"] <- FALSE
  df_scaled_numeric <- as.data.frame(scale(df[, numeric_cols]))
  df_categorical <- df[, !numeric_cols, drop = FALSE]
  df_final <- cbind(df_scaled_numeric, df_categorical)
  df_final <- df_final[, colnames(df)]

  return(df_final)
}

colnames <- "RainTomorrow ~ Date + MinTemp + MaxTemp + Rainfall + Sunshine + WindGustSpeed + WindDir3pm + WindSpeed9am + WindSpeed3pm + Humidity3pm + Temp9am + RainToday"

scaled_train <- scale_with_categoricals(train)
scaled_test <- scale_with_categoricals(test)
fit_and_predict(colnames, scaled_train, scaled_test)

## [1] "Training model..."
## 
## Call:
## glm(formula = as.formula(colnames), family = "binomial", data = train)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.92922    0.45994  -4.194 2.73e-05 ***
## Date           0.12324    0.11027   1.118  0.26371    
## MinTemp       -0.43626    0.25810  -1.690  0.09098 .  
## MaxTemp        1.18830    0.33688   3.527  0.00042 ***
## Rainfall       0.30537    0.12712   2.402  0.01629 *  
## Sunshine      -0.92135    0.16790  -5.487 4.08e-08 ***
## WindGustSpeed  0.89920    0.17734   5.071 3.97e-07 ***
## WindDir3pmENE  0.64936    0.60655   1.071  0.28436    
## WindDir3pmESE  0.95656    0.53520   1.787  0.07389 .  
## WindDir3pmN   -0.40674    0.92595  -0.439  0.66046    
## WindDir3pmNE   0.89735    0.66321   1.353  0.17604    
## WindDir3pmNNE  0.84192    0.94157   0.894  0.37123    
## WindDir3pmNNW  0.06831    0.81855   0.083  0.93350    
## WindDir3pmNW  -0.14457    0.73494  -0.197  0.84405    
## WindDir3pmS    0.66957    0.55759   1.201  0.22982    
## WindDir3pmSE   0.84997    0.52166   1.629  0.10324    
## WindDir3pmSSE  0.45645    0.56781   0.804  0.42147    
## WindDir3pmSSW  1.52247    0.69054   2.205  0.02747 *  
## WindDir3pmSW   1.37986    0.76135   1.812  0.06993 .  
## WindDir3pmW   -0.28554    0.82740  -0.345  0.73001    
## WindDir3pmWNW -0.64347    0.66191  -0.972  0.33098    
## WindDir3pmWSW  2.10450    0.87664   2.401  0.01637 *  
## WindSpeed9am   0.12020    0.14413   0.834  0.40432    
## WindSpeed3pm  -0.22978    0.17491  -1.314  0.18896    
## Humidity3pm    0.99307    0.19796   5.017 5.26e-07 ***
## Temp9am       -0.73199    0.36345  -2.014  0.04401 *  
## RainToday      0.02293    0.13386   0.171  0.86402    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 904.05  on 778  degrees of freedom
## Residual deviance: 593.54  on 752  degrees of freedom
## AIC: 647.54
## 
## Number of Fisher Scoring iterations: 5
## 
## [1] "Predicting Test set..."
## [1] "Test Results:"

AUC	Accuracy	Precision	Recall	F1	MCC
0.8387913	0.8138425	0.7073171	0.5178571	0.5979381	0.4904495

PCA and Neural Network

library(nnet)

dim_fit_predict <- function(cols, train, test) {
  pca_train <- prcomp(
    train[, setdiff(names(train), "RainTomorrow")],
    center = TRUE, scale. = TRUE
  )

  pca_test <- predict(
    pca_train,
    newdata = test[, setdiff(names(test), "RainTomorrow")]
  )

  #print(paste("PCA on Train set:", summary(pca_train)))

  biplot(pca_train, scale = 0)
  plot(pca_train, type = "lines", main = "Scree Plot")

  nn_model <- nnet::nnet(
    x = pca_train$x,
    y = train$RainTomorrow,
    size = 10,
    decay = 1e-4,
    maxit = 300
  )

  nn_preds <- predict(nn_model, pca_test)
  
  show_scores(nn_preds, test)
}

Bonus with NN and PCA

colnames <- "RainTomorrow ~ Date + MinTemp + MaxTemp + Rainfall + Evaporation + Sunshine + WindGustDir + WindGustSpeed + WindDir9am + WindDir3pm + WindSpeed9am + WindSpeed3pm + Humidity9am + Humidity3pm + Pressure9am + Pressure3pm + Cloud9am + Cloud3pm + Temp9am + Temp3pm + RainToday"

dim_fit_predict(
  colnames,
  scaled_train[, sapply(scaled_train, is.numeric)],
  scaled_test[, sapply(scaled_test, is.numeric)]
)

## # weights:  201
## initial  value 143.486091 
## iter  10 value 91.676442
## iter  20 value 65.695160
## iter  30 value 47.227915
## iter  40 value 37.439418
## iter  50 value 35.814082
## iter  60 value 34.730769
## iter  70 value 33.202574
## iter  80 value 32.651937
## iter  90 value 32.410193
## iter 100 value 31.357815
## iter 110 value 31.168381
## iter 120 value 31.068061
## iter 130 value 30.999504
## iter 140 value 30.915752
## iter 150 value 30.844718
## iter 160 value 30.781911
## iter 170 value 30.736843
## iter 180 value 30.685508
## iter 190 value 30.609785
## iter 200 value 30.558699
## iter 210 value 30.035934
## iter 220 value 28.940281
## iter 230 value 28.562348
## iter 240 value 28.523324
## iter 250 value 28.500411
## iter 260 value 28.483036
## iter 270 value 28.470065
## iter 280 value 28.459342
## iter 290 value 28.449187
## iter 300 value 28.440986
## final  value 28.440986 
## stopped after 300 iterations
## [1] "Test Results:"

AUC	Accuracy	Precision	Recall	F1	MCC
0.7894369	0.7637232	0.5656566	0.5	0.5308057	0.3749812

Conclusion

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

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

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

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