Applicazioni pratiche di machine learning/Previsioni sul livello del mare

Indice del libro

Caricamento librerie

modifica
 library(dplyr)
 library(ggplot2)
 library(caret)
 library(reshape2)
 library(forecast)

Parte 1: Dati

modifica

I dati annuali presi dal sito: http://www.psmsl.org/data/obtaining/ contengono serie temporali di misurazioni del livello del mare in varie stazioni del mondo . Tali misurazioni sono in millimetri rispetto a 7 metri al di sotto del livello del mare definito livello "RLR" che sta per "REVISED LOCAL REFERENCE", in modo da evitare la possibilità di ottenere valori negativi. In questa analisi si considerano soltanto le serie temporali relative alle seguenti 2 stazioni:

  • Venezia punta della salute dall'anno 1906 al 2000
  • Venezia II dall'anno 2002 al 2015

Si vuole valutare se nei prossimi anni cioè dal 2020 in poi ci sarà ancora un progressivo innalzamento del livello del mare nelle 2 stazioni di Venezia.


Caricamento dei dati:

 #Dati relativi alle 2 stazioni:

 Venezia_punta_della_salute_168 <- read.csv("Venezia_punta_della_salute_168.csv", sep = ";")
 Venezia_II_2100 <- read.csv("Venezia_II_2100.csv", sep = ";")

 #Ho preparato il file stazioni.csv contenente i dati relativi alle 2 stazioni

 stazioni <- read.csv("stazioni.csv", sep = ";")
 str(stazioni)
'data.frame':	97 obs. of  3 variables:
$ Anno                      : int  1909 1911 1914 1915 1916 1917 1918 1919 1920 1923 ...
$ Venezia.punta.della.salute: int  6904 6933 6938 7007 6998 6983 6930 6985 6925 6953 ...
$ Venezia.II                : int  NA NA NA NA NA NA NA NA NA NA ...


 #Elimino dai dataset i dati mancanti identificati da -99999
 Venezia_punta_della_salute_168<-  Venezia_punta_della_salute_168[-which(Venezia_punta_della_salute_168$Altezza==-99999),]
 stazioni<- stazioni[-which(stazioni$Venezia.punta.della.salute==-99999),]

 #Ottengo un nuovo dataframe contenente la variabile Stazioni
 df <- melt(stazioni,id="Anno", na.rm = TRUE)
 names(df)[2]<-"Stazioni"
 str(df)
'data.frame':	97 obs. of  3 variables:
$ Anno    : int  1909 1911 1914 1915 1916 1917 1918 1919 1920 1923 ...
$ Stazioni: Factor w/ 2 levels "Venezia.punta.della.salute",..: 1 1 1 1 1 1 1 1 1 1 ...
$ value   : int  6904 6933 6938 7007 6998 6983 6930 6985 6925 6953 ...

Parte 2: Esplorazione dei dati

modifica

Il sommario delle 2 stazioni "Venezia punta della salute" e "Venezia II" con minimo, massimo, media, mediana ecc. è il seguente:

 print("Venezia punta della salute")
 summary(Venezia_punta_della_salute_168[,c(1,2)])
 print("Venezia II")
 summary(Venezia_II_2100[,c(1,2)])
[1] "Venezia punta della salute"
     Anno         Altezza    
Min.   :1909   Min.   :6900  
1st Qu.:1935   1st Qu.:6994  
Median :1957   Median :7060  
Mean   :1957   Mean   :7051  
3rd Qu.:1980   3rd Qu.:7114  
Max.   :2000   Max.   :7171  

[1] "Venezia II"
     Anno         Altezza    
Min.   :2002   Min.   :6857  
1st Qu.:2005   1st Qu.:6931  
Median :2008   Median :6971  
Mean   :2008   Mean   :6979  
3rd Qu.:2012   3rd Qu.:7023  
Max.   :2015   Max.   :7108  
 df %>%
  ggplot(aes(Anno,value, colour=Stazioni)) +
  geom_line() +
  scale_x_continuous(breaks=seq(1906,2016,10))+
  ylab("Altezza annuale del mare in millimetri") +
  xlab("Anno") +
  ggtitle("Innalzamento annuale del mare in 2 stazioni di Venezia dal 1909 al 2015",subtitle = "in millimetri rispetto a 7m al di sotto del livello del mare")
 


Parte 3 : Modellizzazione e previsione

modifica

Usando l' algoritmo Bagged Mars per creare un modello previsionale nel dataset relativo alla stazione Venezia punta della salute si ottengono previsioni sulla crescita del livello del mare con un errore RMSE di soli 34,9 mm e un Rsquared del 78,77% che rappresenta la variabilità dell'altezza espressa dal modello.

 model <-  train( Altezza ~ Anno, data=Venezia_punta_della_salute_168, method="bagEarth")
 model

 for (y in 2020:2040) {
  
  new_year <- data.frame(Anno=y)
  print("_________________________________________________________________")
  pred <-predict(model, new_year)
  print(paste("Nell'anno ",y," la predetta altezza annuale sarà ",round(pred[1],2), "."))
  
 }
Loading required package: earth
Loading required package: Formula
Loading required package: plotmo
Loading required package: plotrix
Loading required package: TeachingDemos
Bagged MARS 
83 samples
1 predictor
No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 83, 83, 83, 83, 83, 83, ... 
Resampling results across tuning parameters:
 nprune  RMSE      Rsquared   MAE     
 2       35.65377  0.7769987  28.94041
 5       34.96663  0.7859512  27.79646
 8       34.90546  0.7877331  27.67771
Tuning parameter 'degree' was held constant at a value of 1
RMSE was used to select the optimal model using
the smallest value.
The final values used for the model were nprune = 8
and degree = 1.
[1] "_________________________________________________________________"
[1] "Nell'anno  2020  la predetta altezza annuale sarà  7300.11 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2021  la predetta altezza annuale sarà  7307.28 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2022  la predetta altezza annuale sarà  7314.45 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2023  la predetta altezza annuale sarà  7321.62 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2024  la predetta altezza annuale sarà  7328.79 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2025  la predetta altezza annuale sarà  7335.97 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2026  la predetta altezza annuale sarà  7343.14 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2027  la predetta altezza annuale sarà  7350.31 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2028  la predetta altezza annuale sarà  7357.48 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2029  la predetta altezza annuale sarà  7364.65 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2030  la predetta altezza annuale sarà  7371.82 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2031  la predetta altezza annuale sarà  7378.99 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2032  la predetta altezza annuale sarà  7386.16 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2033  la predetta altezza annuale sarà  7393.33 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2034  la predetta altezza annuale sarà  7400.5 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2035  la predetta altezza annuale sarà  7407.67 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2036  la predetta altezza annuale sarà  7414.84 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2037  la predetta altezza annuale sarà  7422.01 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2038  la predetta altezza annuale sarà  7429.18 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2039  la predetta altezza annuale sarà  7436.36 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2040  la predetta altezza annuale sarà  7443.53 ."

Usando l' algoritmo Bagged Mars per creare un modello previsionale nel dataset relativo alla stazione "Venezia II" si ottengono previsioni sulla crescita del livello del mare con un errore RMSE di soli 49,25 mm e un Rsquared del 75,27% che rappresenta la variabilità dell'altezza espressa dal modello.

 model <-  train( Altezza ~ Anno,data=Venezia_II_2100, method="bagEarth")
 model

 for (y in 2020:2040) {
  
  new_year <- data.frame(Anno=y)
  print("_________________________________________________________________")
  pred <-predict(model, new_year)
  print(paste("Nell'anno ",y," la predetta altezza annuale sarà ",round(pred[1],2), "."))
  
 }
Bagged MARS 
14 samples
1 predictor
No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 14, 14, 14, 14, 14, 14, ... 
Resampling results across tuning parameters:
 nprune  RMSE      Rsquared   MAE     
 2       50.10059  0.7493955  40.03839
 3       49.25730  0.7527612  39.48663
 4       50.57717  0.7470795  40.38399
Tuning parameter 'degree' was held constant at a value of 1
RMSE was used to select the optimal model using
the smallest value.
The final values used for the model were nprune = 3
and degree = 1.
[1] "_________________________________________________________________"
[1] "Nell'anno  2020  la predetta altezza annuale sarà  7106.03 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2021  la predetta altezza annuale sarà  7115.94 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2022  la predetta altezza annuale sarà  7125.85 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2023  la predetta altezza annuale sarà  7135.76 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2024  la predetta altezza annuale sarà  7145.67 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2025  la predetta altezza annuale sarà  7155.58 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2026  la predetta altezza annuale sarà  7165.49 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2027  la predetta altezza annuale sarà  7175.4 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2028  la predetta altezza annuale sarà  7185.31 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2029  la predetta altezza annuale sarà  7195.22 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2030  la predetta altezza annuale sarà  7205.13 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2031  la predetta altezza annuale sarà  7215.04 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2032  la predetta altezza annuale sarà  7224.95 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2033  la predetta altezza annuale sarà  7234.86 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2034  la predetta altezza annuale sarà  7244.77 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2035  la predetta altezza annuale sarà  7254.68 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2036  la predetta altezza annuale sarà  7264.59 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2037  la predetta altezza annuale sarà  7274.5 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2038  la predetta altezza annuale sarà  7284.41 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2039  la predetta altezza annuale sarà  7294.32 ."
[1] "_________________________________________________________________"
[1] "Nell'anno  2040  la predetta altezza annuale sarà  7304.23 ." 


Usando il modello ets della libreria forecast si ottiene un trend crescente dell'altezza per la stazione "Venezia punta della salute" ed un plateau per la stazione "Venezia II" con un intervallo di confidenza all'80% e al 95% rappresentato nei grafici con le aree grigio chiaro e scuro:

 serie <- ts(Venezia_punta_della_salute_168$Altezza, start = 2000-nrow(Venezia_punta_della_salute_168)+1,end = 2000)
 fore<-forecast(serie)
 plot(fore)
 as.data.frame(fore)

 serie <- ts(Venezia_II_2100$Altezza, start = 2002,end = 2015)
 fore<-forecast(serie)
 plot(fore)
 as.data.frame(fore)  
 }


Stazione Venezia punta della salute:

 
     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
2001       7145.584 7098.337 7192.831 7073.326 7217.842
2002       7148.031 7099.695 7196.367 7074.108 7221.955
2003       7150.479 7101.077 7199.881 7074.925 7226.033
2004       7152.926 7102.480 7203.372 7075.775 7230.077
2005       7155.373 7103.903 7206.844 7076.656 7234.091
2006       7157.821 7105.345 7210.296 7077.566 7238.075
2007       7160.268 7106.805 7213.731 7078.504 7242.032
2008       7162.715 7108.282 7217.148 7079.467 7245.963
2009       7165.163 7109.776 7220.550 7080.455 7249.870
2010       7167.610 7111.284 7223.936 7081.467 7253.753


Stazione Venezia II:

 
     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
2016       7039.444 6963.251 7115.637 6922.917 7155.971
2017       7039.444 6950.737 7128.151 6903.778 7175.110
2018       7039.444 6939.782 7139.106 6887.024 7191.864
2019       7039.444 6929.917 7148.971 6871.937 7206.952
2020       7039.444 6920.870 7158.019 6858.100 7220.788
2021       7039.444 6912.465 7166.423 6845.246 7233.642
2022       7039.444 6904.583 7174.305 6833.192 7245.696
2023       7039.444 6897.137 7181.751 6821.804 7257.084
2024       7039.444 6890.061 7188.827 6810.982 7267.906
2025       7039.444 6883.306 7195.583 6800.651 7278.238