vtec_decision_tree_random_forest.py

# -*- coding: utf-8 -*-
"""VTEC_Decision Tree_Random Forest.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1Ja27JM5eodBHo5jXGvh9kOWoJNSYLA0N

### This example demonstrate application of machine learning methods for forecasting Vertical Total Electron Content (VTEC) in the Earth's ionosphere. VTEC data were extracted for single point in mid-latitude at 10E 40N from the Global Ionosphere Map (GIM) provided by CODE, Bern from https://cddis.nasa.gov/archive/gnss/products/ionex/.
### Input data contain indices and measurements of solar activity and solar-terestrial processes, i.e., solar radio flux F10.7, solar wind speed, Bz index of interplenetary geomagnetic field and Dst (disturbance storm time) index of geomagnetic field, downloaded from OMNIWeb Service from NASA at: https://omniweb.gsfc.nasa.gov/form/dx1.html.
### Learning algorithms tested are Decision Tree and ensemble learning of Random Forest.

### **The notebook is prepared by Randa Natras, Technical University of Munich: randa.natras@tum.de **

# **Load packages and data**
"""

import numpy as np
import pandas as pd
import seaborn as sns

from sklearn.tree import export_graphviz
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.model_selection import KFold
from sklearn.metrics import mean_squared_error, mean_absolute_error
from sklearn.model_selection import cross_val_score
from sklearn import metrics

import matplotlib as mpl
import matplotlib.dates as mdates 
import matplotlib.pyplot as plt

import datetime
import pydotplus
import plotly.express as px
from pydotplus import graph_from_dot_data
from IPython.display import Image

"""NumPy showes floating point (real) numbers in the scientific notation. We can turn that off."""

np.set_printoptions(suppress=True)

"""The data will be loaded into a pandas DataFrame, which are good tools for manipulating and displaying the data. Date-time column is assigned as the index column."""

data_df = pd.read_csv ('ML_VTEC_data.csv', index_col='Date-time')

"""# **Examining data**

It's good practice to check the data first.
"""

data_df.head()

data_df.columns

data_df.dtypes

# See the column data types and non-missing values
data_df.info()

data_df.describe()

"""Now we will create output to forecast VTEC 24 hours in future."""

data_df.reset_index(inplace=True)

#prepare VTEC output 24h ahead
VTEC_24h=data_df.loc[24::, ['VTEC 10E 40N']]

VTEC_24h.reset_index(inplace=True)
VTEC_24h=VTEC_24h.drop(['index'], axis=1)

data_df['VTEC 10E 40N (t+24h)'] = VTEC_24h
data_df = data_df.set_index('Date-time')

data_df

"""Check if data contains missing NaN values."""

data_df.isna().sum()

#creating output for 24h forecast resulted in NaN values. We will drop NaN now:
data_df=data_df.dropna()

"""# **Exploratory Data Analysis**

Exploratory Data Analysis is a process of analyzing data to summarize their main characteristics using statistical graphics and other data visualization methods in order to find trends, anomalies, patterns, and relationships within the data. The findings can be used to help us decide which input features to use.

## **Plotting data**
"""

fig, axes = plt.subplots(nrows=6, ncols=1, sharex=True)

data_df.plot(ax=axes[0], y='F10.7 index (sfu)', figsize=(10, 8))
data_df.plot(ax=axes[1], y='Solar wind speed (km/s)', figsize=(10, 8))
data_df.plot(ax=axes[2], y='Bz index (nT)', figsize=(10, 8))
data_df.plot(ax=axes[3], y='Dst index (nT)', figsize=(10, 8))
data_df.plot(ax=axes[4], y='VTEC 10E 40N (TECU)', figsize=(10, 8))
data_df.plot(ax=axes[5], y='VTEC 10E 40N (t+24h)', figsize=(10, 8))
plt.xticks(rotation=45, ha='right')
plt.rcParams.update({'font.size': 14})
plt.show()

"""Plot shows outliers in solar wind speed and Bz index. Values of 999.9 and 9999 usually denotes missing values or instrument failure. We can check how many of these values are contained in our dataset."""

data_df.isin([9999]).sum()

data_df.isin([999.9]).sum()

# Replace all 999.9 and 9999 with NaNs to make it easier to work with
data_df.replace(to_replace=999.9, value=np.nan, inplace=True)
data_df.replace(to_replace=9999, value=np.nan, inplace=True)

data_df.isin([9999]).sum()

data_df.isin([999.9]).sum()

"""Values NaN can be interpolated or dropped. Here we will drop all NaN values."""

data_df = data_df.dropna()

fig, axes = plt.subplots(nrows=5, ncols=1, sharex=True)

data_df.plot(ax=axes[0], y='F10.7 index (sfu)', figsize=(10, 8))
data_df.plot(ax=axes[1], y='Solar wind speed (km/s)', figsize=(10, 8))
data_df.plot(ax=axes[2], y='Bz index (nT)', figsize=(10, 8))
data_df.plot(ax=axes[3], y='Dst index (nT)', figsize=(10, 8))
data_df.plot(ax=axes[4], y='VTEC 10E 40N (TECU)', figsize=(10, 8))
plt.xticks(rotation=45, ha='right')

"""## **Distribution of data**

One important tool for checking our data is data visualization. The goal is to predict the VTEC so a reasonable place to start is examining the distribution of this variable. A histogram is a simple yet effective way to visualize the distribution of a single variable. It is easy to make using seaborn visualization library.
"""

# Histogram of the VTEC
sns.histplot(data=data_df, x='VTEC 10E 40N (TECU)')
sns.set(rc = {'figure.figsize':(5,3)})

# Historgram and kernel density estimate 
# complementary information about the shape of the distribution
sns.histplot(data=data_df, x='VTEC 10E 40N (TECU)', kde=True)
sns.set(rc = {'figure.figsize':(5,3)})

#Plot the distribution with a histogram and maximum likelihood gaussian distribution fit
from scipy.stats import norm
sns.distplot(data_df['VTEC 10E 40N (TECU)'], fit=norm, kde=False)
sns.set(rc = {'figure.figsize':(5,3)})

data_df_2=data_df.drop(['Hour of day', 'Day of year'], axis=1)

#histogram of input observation data
plt.style.use('default')
plt.rcParams.update({'font.size': 12})
data_df_2.hist(figsize=(10,10), xlabelsize=12, ylabelsize=12)

"""A matrix of scatter plots can be use to check how much one variable is affected by another and show the (linear) relationship between two sets of observations."""

plt.style.use('default')
sns.pairplot(data_df)

"""## Correlation analysis"""

#Correlation table
np.round(data_df.corr(), 2)

c = np.round(data_df.corr(), 2)
plt.figure(figsize=(15,15))
sns.heatmap(c, annot=True, vmin=-1, vmax=1, cmap='coolwarm', square=True)
plt.rcParams.update({'font.size':20})

# Find all correlations with the output and sort 
correlations =data_df.corr()['VTEC 10E 40N (t+24h)'].sort_values()
correlations

"""# **Preparing the data for training**

Data will be prepared as NumPy arrays that will be used by scikit-learn algorithms. X will be a matrix which has different datapoints in different rows and different input features in different columns. Labels (or targets) y will be a vector consisting of 'VTEC 10E 40N (t+24h)' column.

Cyclic continues features such as Day of year and Hour of day are sometimes encoded to sine and cosine features to perserve their cyclical nature. This is required for the learning algorithms that compute Euclidean distances in order to have two similar features nearby in feature space. However, the decision tree works differently by splitting the dataset from one region into two non-overlapping sub-regions based on a single splitting variable and the split point (more details in https://www.mdpi.com/2072-4292/14/15/3547). Therefore, no sin/cos encoding is performed in this example. Studies on cyclic data transformation (such as  https://scholarworks.calstate.edu/downloads/pv63g5147) show that decision trees do not benefit from sin/cos encoding, resulting in similar performance with or without cycle data transformation.
"""

X = data_df.drop(['VTEC 10E 40N (t+24h)'], axis=1).to_numpy()
y = data_df ['VTEC 10E 40N (t+24h)'].to_numpy()

X = data_df.drop(['VTEC 10E 40N', 'VTEC 10E 40N (t+24h)'], axis=1).to_numpy()
y = data_df ['VTEC 10E 40N'].to_numpy()

dates = pd.to_datetime(data_df.index)

#we will perform testing on for year 2019 (from 1 January to 30 December) which is 8736 data points
#test data need to be excluded from the training data
split_time = 8736
X_train = X [:-split_time]
y_train = y [:-split_time]
time_train = dates [:-split_time]

X_test= X [-split_time:]
y_test = y [-split_time:]
time_test = dates [-split_time:]

time_train

time_test

y_train.size

y_test.size

fig = plt.figure()
ax = fig.add_subplot()
ax.plot(time_train, y_train, label='Train VTEC')
ax.plot(time_test, y_test, label='Test VTEC')
plt.ylabel('VTEC (TECU)')
plt.xlabel('Time (UTC)')
plt.legend(loc="best", bbox_to_anchor=(0.5, 0.5, 0.5, 0.5))
plt.xticks(rotation=45, ha='right')
plt.show()
plt.rcParams.update({'font.size': 14})
plt.rcParams['figure.figsize'] = (12,5)

"""# Model Tuning

## Splitting the data

Data are splitted applying time series cross-validation splitting based on Kfold with 20 splits as in https://www.mdpi.com/2072-4292/14/15/3547.
"""

from sklearn.model_selection import TimeSeriesSplit
tscv = TimeSeriesSplit(n_splits=20)

"""## Decision Tree Regressor Tuning

We can use cross_val_score() along with KFold to evaluate the model for different hyperparameters. Here we are going to try different hyperparameter values and choose the ones for which we get the highest model score.

There are multiple hyperparameters such as max_depth, min_samples_split, min_samples_leaf etc. which affect the model performance. More information can be found, for example in https://www.mdpi.com/2072-4292/14/15/3547. Here we are going to tune the model based on 'max_depth'. We will try different values for max_depth values from 1 to 20. Depending on the final 'RMSE' score, we will choose the final value of max_depth.
"""

def rmse(score):
    rmse = np.sqrt(-score)
    print(f'rmse= {"{:.2f}".format(rmse)}')

max_depth = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20]

for val in max_depth:
    score = cross_val_score(DecisionTreeRegressor(max_depth= val), X_train, y_train, cv= tscv, scoring="neg_mean_squared_error")
    print(f'For max depth: {val}')
    rmse(score.mean())

"""## Random Forest Regressor Tuning

There are multiple hyperparameters like n_estimators, max_depth, max_features, min_samples_split etc. which affect the model performance. Here we will tune the 'n_estimators'. We will try different range of values for max_depth, max_features and number of estimators. Depending on the final 'RMSE' score, we will choose the value of these hyperparameters.
"""

max_depth = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 15, 20, 25, 50, None]

for count in max_depth:
    score = cross_val_score(RandomForestRegressor(max_depth = count), X_train, y_train, cv = tscv, scoring="neg_mean_squared_error")
    print(f'For max depth: {count}')
    rmse(score.mean())

max_features = [1, 2, 3, 4, 5, 6, 7]

for count in max_features:
    score = cross_val_score(RandomForestRegressor(max_features = count, max_depth = 7), X_train, y_train, cv = tscv, scoring="neg_mean_squared_error")
    print(f'For max_features: {count}')
    rmse(score.mean())

estimators = [10, 50, 100, 200, 300, 400, 500, 1000]

for count in estimators:
    score = cross_val_score(RandomForestRegressor(n_estimators = count, max_features = 3, max_depth = 7), X_train, y_train, cv = tscv, scoring="neg_mean_squared_error")
    print(f'For estimators: {count}')
    rmse(score.mean())

"""# **Decision Tree**

## Model training
"""

model_dtree = DecisionTreeRegressor(max_depth=4)

model_dtree=model_dtree.fit(X_train, y_train)
y_pred_train=model_dtree.predict(X_train)

print ('Metrics on train data')
print('Root mean square error (RMSE):', round(np.sqrt(mean_squared_error(y_train,y_pred_train)),2))
print("Explain variance score =", round(metrics.explained_variance_score(y_train,y_pred_train), 2)) 
print("R2 score =", round(metrics.r2_score(y_train,y_pred_train), 2))

dtree_train_rmse = np.sqrt(mean_squared_error(y_train,y_pred_train))

"""## **Visualizing the Decision Tree**

Now, let's visualize our Decision Tree. Let's export our model as a special kind of data, create a visual representation from that, generate a graph and show that graph as an image:
"""

import pydotplus
import Image
dot_data = 'tree.dot'
export_graphviz(model_dtree, out_file=dot_data,feature_names=data_df.columns[:-1], filled=True, rounded=True, impurity=False, special_characters=True)
graph=pydotplus.graph_from_dot_file(dot_data)
Image(graph.create_png(), unconfined=True)

import pydotplus
import plotly.express as px
from pydotplus import graph_from_dot_data
from IPython.display import Image

dot_data = 'tree.dot'
export_graphviz(model_dtree, out_file=dot_data,feature_names=data_df.columns[:-1], filled=True, rounded=True, impurity=True, special_characters=True)
graph=pydotplus.graph_from_dot_file(dot_data)
Image(graph.create_png(), unconfined=True)

model_dtree = DecisionTreeRegressor(max_depth=3)

model_dtree=model_dtree.fit(X_train, y_train)
y_pred_train=model_dtree.predict(X_train)

dot_data = 'tree.dot'
export_graphviz(model_dtree, out_file=dot_data,feature_names=data_df.columns[:-2], filled=True, rounded=True, impurity=True, special_characters=True)
graph=pydotplus.graph_from_dot_file(dot_data)
plt.rcParams.update({'font.size': 14})
plt.rcParams['figure.figsize'] = (12,4)
Image(graph.create_png(), unconfined=True)

dot_data = 'tree.dot'
export_graphviz(model_dtree, out_file=dot_data,feature_names=data_df.columns[:-2], filled=True, rounded=True, impurity=True, special_characters=True)
graph=pydotplus.graph_from_dot_file(dot_data)
Image(graph.create_png(), unconfined=True)

features = data_df.columns[:-2]

importance = model_dtree.feature_importances_
  
  # summarize feature importance
  for i,v in enumerate(importance):
    print('Feature: %0d, Relative importance: %.5f' % (i,v))

  # plot feature importance
  plt.barh([x_train for x_train in range(len(importance))], importance, align='center')
  plt.yticks(range(len(features)), features)
  plt.xticks((np.arange(0.0, 1.1, 0.2)))
  plt.rcParams ['figure.figsize'] = [5, 3]
  plt.rcParams.update({'font.size': 14})
  plt.xlabel('Relative importance')
  plt.show()

importance = model_dtree.feature_importances_
  
  # summarize feature importance
  for i,v in enumerate(importance):
    print('Feature: %0d, Relative importance: %.5f' % (i,v))

  # plot feature importance
  plt.barh([x_train for x_train in range(len(importance))], importance, align='center')
  plt.yticks(range(len(features)), features)
  plt.xticks((np.arange(0.0, 1.1, 0.2)))
  plt.rcParams ['figure.figsize'] = [5, 3]
  plt.rcParams.update({'font.size': 14})
  plt.xlabel('Relative importance')
  plt.show()

"""## **Testing the model**

Finally, let's evaluate our model on the test data to get a final accuracy performance score for our model.
"""

y_pred_test=model_dtree.predict (X_test)

print ('Metrics on test data:')
print('Root mean squared error (RMSE):', round(np.sqrt(metrics.mean_squared_error(y_test,y_pred_test)),2))
print("Explain variance score =", round(metrics.explained_variance_score(y_test,y_pred_test), 2)) 
print("R2 score =", round(metrics.r2_score(y_test,y_pred_test), 2))

dtree_test_rmse = np.sqrt(mean_squared_error(y_test,y_pred_test))

#plot for year 2019: January - December 2019
plt.plot(time_test, y_test, 'tab:orange', label='Ground truth VTEC')
plt.plot(time_test, y_pred_test, 'tab:green', label='VTEC 24h forecast, DT')
plt.ylabel('VTEC (TECU)')
plt.xlabel('Time (UTC)')
plt.xticks(rotation=45, ha='right')
plt.yticks((np.arange(0.0, 31.0, 5.0)))
plt.legend(loc="best", bbox_to_anchor=(0.5, 0.5, 0.5, 0.5))
plt.show()
plt.rcParams.update({'font.size': 14})
plt.rcParams['figure.figsize'] = (12,4)

#plot for 8 days: 23 - 30 December 2019
plt.plot(time_test[-190:], y_test[-190:], 'tab:orange',label='Ground truth VTEC')
plt.plot(time_test[-190:], y_pred_test[-190:], 'tab:green', label='VTEC 24h forecast, DT')
plt.ylabel('VTEC (TECU)')
plt.xlabel('Time (UTC)')
plt.xticks(rotation=45, ha='right')
plt.yticks((np.arange(0.0, 16.0, 5.0)))
plt.legend(loc="best", bbox_to_anchor=(0.5, 0, 0.5, 0.5))
plt.show()
plt.rcParams.update({'font.size': 14})

"""# **Random Forest**

Let's try to improve the score using an ensemble of decision trees in a form of Random forest. The example below represents Random forest model with hyperparameters that minimize RMSE, which we found previously.

## **Model training**
"""

model_rfr = RandomForestRegressor(max_features = 3, max_depth = 7, n_estimators = 500)

model_rfr = model_rfr.fit(X_train, y_train)
y_pred_train_rf = model_rfr .predict(X_train)

print ('Metrics on train data')
print('Root mean square error (RMSE):', round(np.sqrt(mean_squared_error(y_train,y_pred_train_rf)),2))
print("Explain variance score =", round(metrics.explained_variance_score(y_train,y_pred_train_rf), 2)) 
print("R2 score =", round(metrics.r2_score(y_train,y_pred_train_rf), 2))

rfr_train_rmse = np.sqrt(mean_squared_error(y_train, y_pred_train_rf))

"""## **Feature importance**

Ensembles of decision trees, like Random forest can be used to calculate a feature importance score. It is useful to estimate the relative importance of input features when developing predictive models.
"""

features = data_df.columns[:-1]

importance = model_rfr.feature_importances_
  
  # summarize feature importance
  for i,v in enumerate(importance):
    print('Feature: %0d, Score: %.5f' % (i,v))

  # plot feature importance
  plt.barh([x_train for x_train in range(len(importance))], importance, align='center')
  plt.yticks(range(len(features)), features)
  plt.xticks((np.arange(0.0, 0.8, 0.2)))
  plt.rcParams ['figure.figsize'] = [9, 4]
  plt.rcParams.update({'font.size': 14})
  plt.xlabel('Relative importance')
  plt.show()

"""## **Testing the model**

Finally, let's evaluate our model on the test data to get a final accuracy performance score for our model.
"""

y_pred_test_rf = model_rfr.predict (X_test)

print ('Metrics on test data:')
print('Root mean squared error (RMSE):', round(np.sqrt(metrics.mean_squared_error(y_test,y_pred_test_rf)),2))
print("Explain variance score =", round(metrics.explained_variance_score(y_test,y_pred_test_rf), 2)) 
print("R2 score =", round(metrics.r2_score(y_test,y_pred_test_rf), 2))

rfr_test_rmse = np.sqrt(mean_squared_error(y_test,y_pred_test_rf))

plt.plot(time_test, y_test, 'tab:orange',label='Ground truth VTEC')
plt.plot(time_test, y_pred_test_rf, 'tab:green', label='VTEC 24h forecast, RF')
plt.ylabel('VTEC (TECU)')
plt.xlabel('Time (UTC)')
plt.yticks((np.arange(0.0, 31.0, 5.0)))
plt.xticks(rotation=45, ha='right')
plt.legend(loc="best", bbox_to_anchor=(0.5, 0.5, 0.5, 0.5))
plt.show()
plt.rcParams.update({'font.size': 14})
plt.rcParams['figure.figsize'] = (12,4)

#forecast for 8 days at the end of December 2019
plt.plot(time_test[-190:], y_test[-190:], 'tab:orange',label='Ground truth VTEC')
plt.plot(time_test[-190:], y_pred_test_rf[-190:], 'tab:green', label='VTEC 24h forecast, RF')
plt.ylabel('VTEC (TECU)')
plt.xlabel('Time (UTC)')
plt.yticks((np.arange(4.0, 13.0, 2.0)))
plt.legend(loc="best", bbox_to_anchor=(0.5, 0, 0.5, 0.5))
plt.xticks(rotation=45, ha='right')
plt.yticks((np.arange(0.0, 16.0, 5.0)))
plt.show()
plt.rcParams.update({'font.size': 14})

"""# **Comparison: Decision Tree and Random Forest**

We will calculate differences between ground truth VTEC and VTEC forecasts by Decision Tree and Random Forest models. Using the seaborn library, we will then compare histogram plots between two models and see where most of the differences are located.
"""

diff_dtree = y_test - y_pred_test
diff_rf = y_test - y_pred_test_rf

diff_rf = y_test - y_pred_test_rf

sns.histplot(diff_dtree, log_scale=False, element="step", fill=False, label='Decision Tree')
sns.histplot(diff_rf, log_scale=False, element="step", fill=False, label='Random Forest')
plt.legend(loc="best", bbox_to_anchor=(0.5, 0.5, 0.5, 0.5))
plt.xlabel('Differences between ground truth and forecast (TECU)')
plt.rcParams ['figure.figsize'] = [7.0, 4.0]

sns.histplot(diff_dtree, log_scale=False, element="step", fill=False, label='Decision Tree', stat="density")
sns.histplot(diff_rf, log_scale=False, element="step", fill=False, label='Random Forest', stat="density")
plt.legend(loc="best", bbox_to_anchor=(0.5, 0.5, 0.5, 0.5))
plt.xlabel('Differences between ground truth and forecast (TECU)')
plt.rcParams ['figure.figsize'] = [7.0, 4.0]

x = ['Train', 'Test']
y_dt = [dtree_train_rmse, dtree_test_rmse]
y_rf = [rfr_train_rmse, rfr_test_rmse]
df = pd.DataFrame({'Decision Tree': y_dt,
                   'Random Forest': y_rf}, index=x)
ax = df.plot.bar(ylabel='RMSE (TECU)')
plt.xticks(rotation=0, ha='right')
ax.legend(loc="best", bbox_to_anchor=(1, 0.8, 0, 0.5))
ax.set_yticks((np.arange(0.0, 2.1, 0.5)))