Applicazioni pratiche di machine learning/Previsioni sul livello del mare
Caricamento librerie
modifica library(dplyr)
library(ggplot2)
library(caret)
library(reshape2)
library(forecast)
Parte 1: Dati
modificaI 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
modificaIl 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
modificaUsando 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