1. AI Automation: Data Science Tasks, Streamlined

AI is now a powerhouse for automating the once labor-intensive aspects of data science.

• Smarter Data Prep: AI tools, including AutoML platforms like Google AutoML and generative AI, are automating significant chunks of data cleaning, preprocessing, and even feature engineering. Generative AI can create synthetic data to bolster limited datasets, improving model accuracy and reducing bias.
• Effortless Model Building: AI excels at automated model selection and hyperparameter tuning, sifting through countless options faster than any human.
• AI-Assisted Coding: Tools like OpenAI’s Codex and GitHub Copilot generate code snippets, speeding up development and allowing data scientists to focus on higher-level problem-solving.
• Automated Insights: AI analytics can now perform root cause analysis and generate reports with natural language narratives, making complex data understandable to a wider audience. This automation frees up data scientists from 30–50% of routine tasks, allowing them to tackle more strategic challenges.

2. Beyond Algorithms: The New Premium on Human Skills

As AI handles the mechanics, the data scientist’s role elevates. The new focus is on:

• Strategic Problem Formulation: Identifying which business problems to solve with AI and why becomes paramount.
• Critical Thinking & Interpretation: Evaluating AI outputs, understanding their limitations, and providing nuanced contextual interpretation is key.
• Compelling Data Storytelling: Translating complex AI-driven insights into clear, actionable narratives for non-technical stakeholders is crucial.
• Collaboration & Domain Expertise: Working across teams and leveraging deep industry knowledge ensures AI solutions are relevant and impactful. These “soft skills” are becoming the hard differentiators in the AI era.

3. The AI-Augmented Toolkit: Essential New Competencies

While foundational skills (Python, SQL, R, statistics) remain vital , new technical skills are essential for the AI-augmented data scientist:

• Prompt Engineering: Skillfully guiding LLMs and generative AI tools to produce desired outputs is a new cornerstone.
• MLOps Understanding: Familiarity with deploying, monitoring, and managing models in production is non-negotiable.
• AI Model Interpretability (XAI): Using tools like LIME and SHAP to understand and explain “black box” AI decisions is critical for trust, fairness, and debugging.
• Advanced AI/ML Techniques: Expertise in deep learning, NLP, and frameworks like PyTorch and TensorFlow is increasingly demanded.
• Cloud Computing: Proficiency with AWS, Azure, or GCP is standard for handling the scale of modern AI.

4. The Evolving Mandate: Driving Value, Guarding Ethics

The data scientist’s mission is expanding:

• Delivering Measurable Business Value: The focus is shifting from technical outputs to driving tangible business growth and innovation through AI. Companies with clear AI strategies, often shaped by data scientists, see significantly better ROI.
• Championing Ethical AI: Proactively addressing fairness, transparency, bias, and accountability in AI systems is now a core responsibility. This includes being aware of AI’s potential to perpetuate societal biases or generate misleading information.
• Shaping Data & AI Governance: Data scientists are increasingly involved in defining data quality standards, AI-readiness, and ethical AI frameworks.

5. The Data Scientist of Tomorrow: Intelligent Collaboration

The future is about human-AI collaboration, not competition. This synergy is creating new and evolving roles:

• Lead/Principal Data Scientist: Focusing on AI strategy, team leadership, and ethical oversight.
• AI Ethicist: Specializing in guiding responsible AI development and deployment.
• MLOps Engineer: Ensuring the smooth operationalization of machine learning models.
• AI Solutions Architect: Designing end-to-end AI systems for specific business needs. Continuous learning and adaptability are crucial, as some analyses suggest up to 40% of roles may need redefinition due to AI.

Conclusion: Embrace the AI Metamorphosis

The silent transformation of the data scientist is an evolution, not an extinction. AI amplifies human capabilities, allowing data scientists to focus on more strategic, impactful, and ethically complex work. To thrive, data scientists must:

• Commit to Lifelong Learning: Stay ahead of rapid AI advancements.
• Cultivate Human-Centric Skills: Double down on strategy, critical thinking, communication, and ethics
• Become Strategic Influencers: Guide AI adoption responsibly and drive real business value.

The AI era calls for data scientists who are not just technically proficient but are also insightful strategists, compelling communicators, and ethical guardians. By embracing this metamorphosis, data scientists will continue to be indispensable in shaping an AI-powered future.

Thanks for reading this article. You can access the detailed codes of the project and other projects on my Github account or Kaggle account. Happy coding!

Please feel free to contact me if you need any further information.

Visit us at DataDrivenInvestor.com

Subscribe to DDIntel here.

Join our creator ecosystem here.

DDI Official Telegram Channel: https://t.me/+tafUp6ecEys4YjQ1

Follow us on LinkedIn, Twitter, YouTube, and Facebook.

AI’s Silent Revolution: How Data Scientists Are Evolving (and Thriving) was originally published in DataDrivenInvestor on Medium, where people are continuing the conversation by highlighting and responding to this story.

Let’s implement it in python.

About Dataset

Context

This dataset is originally from the National Institute of Diabetes and Digestive and Kidney Diseases. The objective of the dataset is to diagnostically predict whether or not a patient has diabetes, based on certain diagnostic measurements included in the dataset. Several constraints were placed on the selection of these instances from a larger database. In particular, all patients here are females at least 21 years old of Pima Indian heritage.

Content

The datasets consists of several medical predictor variables and one target variable, Outcome. Predictor variables includes the number of pregnancies the patient has had, their BMI, insulin level, age, and so on.

Can you build a machine learning model to accurately predict whether or not the patients in the dataset have diabetes or not?

1. Importing the libraries

import joblib
import pandas as pd
import seaborn as sns
from matplotlib import pyplot as plt
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier, VotingClassifier, AdaBoostClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_validate, GridSearchCV
from sklearn.neighbors import KNeighborsClassifier
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier
from sklearn.preprocessing import StandardScaler

from xgboost import XGBClassifier
from lightgbm import LGBMClassifier
from catboost import CatBoostClassifier
pd.set_option('display.max_columns', None)
pd.set_option('display.width', 500)
import warnings
warnings.filterwarnings('ignore')

2. Importing the dataset

df = pd.read_csv('/kaggle/input/pima-indians-diabetes-database/diabetes.csv')

# Checking the dataset
df.head()

Pregnancies Glucose BloodPressure SkinThickness Insulin BMI DiabetesPedigreeFunction Age Outcome
0               6       148         72 35 0 33.6 0.627 50 1
1               1       85          66 29 0 26.6 0.351 31 0
2               8       183         64 0 0 23.3 0.672 32 1
3               1       89          66 23 94 28.1 0.167 21 0
4               0       137         40 35 168 43.1 2.288 33 1

3. Exploratory Data Analysis

This is a basic function that provides a high-level overview of a dataframe. The function “check_df” takes two arguments: a dataframe and the number of head and tail rows to be displayed (default value is 5).

The function prints the shape of the dataframe, the data types of each column, the first five and last five rows of the dataframe, the count of missing values in each column, and the quantiles of each column. The output provides a quick way to check for basic issues such as missing values, data types, and distributions of the columns in the dataframe.

def check_df(dataframe, head=5):
    print("##################### Shape #####################")
    print(dataframe.shape)
    print("##################### Types #####################")
    print(dataframe.dtypes)
    print("##################### Head #####################")
    print(dataframe.head(head))
    print("##################### Tail #####################")
    print(dataframe.tail(head))
    print("##################### NA #####################")
    print(dataframe.isnull().sum())
    print("##################### Quantiles #####################")
    print(dataframe.quantile([0, 0.05, 0.50, 0.95, 0.99, 1]).T)
check_df(df)

##################### Shape #####################
(768, 9)
##################### Types #####################
Pregnancies                   int64
Glucose                       int64
BloodPressure                 int64
SkinThickness                 int64
Insulin                       int64
BMI                         float64
DiabetesPedigreeFunction    float64
Age                           int64
Outcome                       int64
dtype: object
##################### Head #####################
   Pregnancies  Glucose  BloodPressure  SkinThickness  Insulin   BMI  DiabetesPedigreeFunction  Age  Outcome
0            6      148             72             35        0  33.6                     0.627   50        1
1            1       85             66             29        0  26.6                     0.351   31        0
2            8      183             64              0        0  23.3                     0.672   32        1
3            1       89             66             23       94  28.1                     0.167   21        0
4            0      137             40             35      168  43.1                     2.288   33        1
##################### Tail #####################
     Pregnancies  Glucose  BloodPressure  SkinThickness  Insulin   BMI  DiabetesPedigreeFunction  Age  Outcome
763           10      101             76             48      180  32.9                     0.171   63        0
764            2      122             70             27        0  36.8                     0.340   27        0
765            5      121             72             23      112  26.2                     0.245   30        0
766            1      126             60              0        0  30.1                     0.349   47        1
767            1       93             70             31        0  30.4                     0.315   23        0
##################### NA #####################
Pregnancies                 0
Glucose                     0
BloodPressure               0
SkinThickness               0
Insulin                     0
BMI                         0
DiabetesPedigreeFunction    0
Age                         0
Outcome                     0
dtype: int64
##################### Quantiles #####################
                            0.00      0.05      0.50       0.95       0.99    1.00
Pregnancies                0.000   0.00000    3.0000   10.00000   13.00000   17.00
Glucose                    0.000  79.00000  117.0000  181.00000  196.00000  199.00
BloodPressure              0.000  38.70000   72.0000   90.00000  106.00000  122.00
SkinThickness              0.000   0.00000   23.0000   44.00000   51.33000   99.00
Insulin                    0.000   0.00000   30.5000  293.00000  519.90000  846.00
BMI                        0.000  21.80000   32.0000   44.39500   50.75900   67.10
DiabetesPedigreeFunction   0.078   0.14035    0.3725    1.13285    1.69833    2.42
Age                       21.000  21.00000   29.0000   58.00000   67.00000   81.00
Outcome                    0.000   0.00000    0.0000    1.00000    1.00000    1.00

Pandas Profiling

pandas-profiling primary goal is to provide a one-line Exploratory Data Analysis (EDA) experience in a consistent and fast solution. Like pandas df.describe() function, that is so handy, pandas-profiling delivers an extended analysis of a DataFrame while alllowing the data analysis to be exported in different formats such as html and json.

!pip install pandas-profiling
import pandas_profiling
pandas_profiling.ProfileReport(df)

Correlations

Spearman’s ρ

The Spearman’s rank correlation coefficient (ρ) is a measure of monotonic correlation between two variables, and is therefore better in catching nonlinear monotonic correlations than Pearson’s r. It’s value lies between -1 and +1, -1 indicating total negative monotonic correlation, 0 indicating no monotonic correlation and 1 indicating total positive monotonic correlation.

To calculate ρ for two variables X and Y, one divides the covariance of the rank variables of X and Y by the product of their standard deviations.

Skimpy

A light weight tool for creating summary statistics from dataframes.

skimpy is a light weight tool that provides summary statistics about variables in data frames within the console or your interactive Python window. Think of it as a super-charged version of df.describe(). You can find the documentation here.

You can install the latest release of skimpy via pip from PyPI:

!pip install skimpy

skim a dataframe and produce summary statistics within the console using:

from skimpy import skim
skim(df)

4. Data Preprocessing & Feature Engineering

1. outlier_thresholds: This function calculates the lower and upper limits for identifying outliers in a given column of a dataframe. It takes in the dataframe and the name of the col_name to calculate the limits for. The q1 and q3 parameters define the quartiles to be used for calculating the lower and upper limits and their default values are 0.25 and 0.75 respectively.
2. replace_with_thresholds: This function replaces the values in a given column of the dataframe that fall outside the lower and upper limits calculated in outlier_thresholds with the corresponding limit. It takes in the dataframe and the name of the column to replace values in as variable.
3. check_outlier: This function checks if there are any outliers in a given column of the dataframe. It returns True if there are any outliers, and False otherwise.
4. one_hot_encoder: This function applies one-hot encoding to a set of categorical columns in a dataframe. It takes in the dataframe and a list of categorical_cols to be one-hot encoded, as well as a boolean flag drop_first that determines whether or not to drop the first column of each categorical variable.
5. df.head(): This line displays the first 5 rows of the df dataframe, including the columns and their values in each row.

def outlier_thresholds(dataframe, col_name, q1=0.25, q3=0.75):
    quartile1 = dataframe[col_name].quantile(q1)
    quartile3 = dataframe[col_name].quantile(q3)
    interquantile_range = quartile3 - quartile1
    up_limit = quartile3 + 1.5 * interquantile_range
    low_limit = quartile1 - 1.5 * interquantile_range
    return low_limit, up_limit

def replace_with_thresholds(dataframe, variable):
    low_limit, up_limit = outlier_thresholds(dataframe, variable)
    dataframe.loc[(dataframe[variable] < low_limit), variable] = low_limit
    dataframe.loc[(dataframe[variable] > up_limit), variable] = up_limit

def check_outlier(dataframe, col_name, q1=0.25, q3=0.75):
    low_limit, up_limit = outlier_thresholds(dataframe, col_name, q1, q3)
    if dataframe[(dataframe[col_name] > up_limit) | (dataframe[col_name] < low_limit)].any(axis=None):
        return True
    else:
        return False

def one_hot_encoder(dataframe, categorical_cols, drop_first=False):
    dataframe = pd.get_dummies(dataframe, columns=categorical_cols, drop_first=drop_first)
    return dataframe

df.head()

df.columns = [col.upper() for col in df.columns]

df['NEW_GLUCOSE_CAT'] = pd.cut(x=df['GLUCOSE'], bins=[-1, 139, 200], labels=["normal", "prediabetes"])

# Age
df.loc[(df['AGE'] < 35), "NEW_AGE_CAT"] = 'young'
df.loc[(df['AGE'] >= 35) & (df['AGE'] <= 55), "NEW_AGE_CAT"] = 'middleage'
df.loc[(df['AGE'] > 55), "NEW_AGE_CAT"] = 'old'

# BMI
df['NEW_BMI_RANGE'] = pd.cut(x=df['BMI'], bins=[-1, 18.5, 24.9, 29.9, 100],
                             labels=["underweight", "healty", "overweight", "obese"])

# BloodPressure
df['NEW_BLOODPRESSURE'] = pd.cut(x=df['BLOODPRESSURE'], bins=[-1, 79, 89, 123], labels=["normal", "hs1", "hs2"])

check_df(df)

##################### Shape #####################
(768, 13)
##################### Types #####################
PREGNANCIES                    int64
GLUCOSE                        int64
BLOODPRESSURE                  int64
SKINTHICKNESS                  int64
INSULIN                        int64
BMI                          float64
DIABETESPEDIGREEFUNCTION     float64
AGE                            int64
OUTCOME                        int64
NEW_GLUCOSE_CAT             category
NEW_AGE_CAT                   object
NEW_BMI_RANGE               category
NEW_BLOODPRESSURE           category
dtype: object
##################### Head #####################
   PREGNANCIES  GLUCOSE  BLOODPRESSURE  SKINTHICKNESS  INSULIN   BMI  DIABETESPEDIGREEFUNCTION  AGE  OUTCOME NEW_GLUCOSE_CAT NEW_AGE_CAT NEW_BMI_RANGE NEW_BLOODPRESSURE
0            6      148             72             35        0  33.6                     0.627   50        1     prediabetes   middleage         obese            normal
1            1       85             66             29        0  26.6                     0.351   31        0          normal       young    overweight            normal
2            8      183             64              0        0  23.3                     0.672   32        1     prediabetes       young        healty            normal
3            1       89             66             23       94  28.1                     0.167   21        0          normal       young    overweight            normal
4            0      137             40             35      168  43.1                     2.288   33        1          normal       young         obese            normal
##################### Tail #####################
     PREGNANCIES  GLUCOSE  BLOODPRESSURE  SKINTHICKNESS  INSULIN   BMI  DIABETESPEDIGREEFUNCTION  AGE  OUTCOME NEW_GLUCOSE_CAT NEW_AGE_CAT NEW_BMI_RANGE NEW_BLOODPRESSURE
763           10      101             76             48      180  32.9                     0.171   63        0          normal         old         obese            normal
764            2      122             70             27        0  36.8                     0.340   27        0          normal       young         obese            normal
765            5      121             72             23      112  26.2                     0.245   30        0          normal       young    overweight            normal
766            1      126             60              0        0  30.1                     0.349   47        1          normal   middleage         obese            normal
767            1       93             70             31        0  30.4                     0.315   23        0          normal       young         obese            normal
##################### NA #####################
PREGNANCIES                 0
GLUCOSE                     0
BLOODPRESSURE               0
SKINTHICKNESS               0
INSULIN                     0
BMI                         0
DIABETESPEDIGREEFUNCTION    0
AGE                         0
OUTCOME                     0
NEW_GLUCOSE_CAT             0
NEW_AGE_CAT                 0
NEW_BMI_RANGE               0
NEW_BLOODPRESSURE           0
dtype: int64
##################### Quantiles #####################
                            0.00      0.05      0.50       0.95       0.99    1.00
PREGNANCIES                0.000   0.00000    3.0000   10.00000   13.00000   17.00
GLUCOSE                    0.000  79.00000  117.0000  181.00000  196.00000  199.00
BLOODPRESSURE              0.000  38.70000   72.0000   90.00000  106.00000  122.00
SKINTHICKNESS              0.000   0.00000   23.0000   44.00000   51.33000   99.00
INSULIN                    0.000   0.00000   30.5000  293.00000  519.90000  846.00
BMI                        0.000  21.80000   32.0000   44.39500   50.75900   67.10
DIABETESPEDIGREEFUNCTION   0.078   0.14035    0.3725    1.13285    1.69833    2.42
AGE                       21.000  21.00000   29.0000   58.00000   67.00000   81.00
OUTCOME                    0.000   0.00000    0.0000    1.00000    1.00000    1.00

def grab_col_names(dataframe, cat_th=10, car_th=20):

    # cat_cols, cat_but_car
    cat_cols = [col for col in dataframe.columns if dataframe[col].dtypes == "O"]
    num_but_cat = [col for col in dataframe.columns if dataframe[col].nunique() < cat_th and
                   dataframe[col].dtypes != "O"]
    cat_but_car = [col for col in dataframe.columns if dataframe[col].nunique() > car_th and
                   dataframe[col].dtypes == "O"]
    cat_cols = cat_cols + num_but_cat
    cat_cols = [col for col in cat_cols if col not in cat_but_car]

    # num_cols
    num_cols = [col for col in dataframe.columns if dataframe[col].dtypes != "O"]
    num_cols = [col for col in num_cols if col not in num_but_cat]

    return cat_cols, num_cols, cat_but_car

check_df(df)

cat_cols, num_cols, cat_but_car = grab_col_names(df, cat_th=5, car_th=20)
cat_cols = [col for col in cat_cols if "OUTCOME" not in col]

for col in num_cols:
    print(col, check_outlier(df, col, 0.05, 0.95))

replace_with_thresholds(df, "INSULIN")

PREGNANCIES False
GLUCOSE False
BLOODPRESSURE False
SKINTHICKNESS False
INSULIN True
BMI False
DIABETESPEDIGREEFUNCTION False
AGE False

cat_cols, num_cols, cat_but_car = grab_col_names(df, cat_th=5, car_th=20)

# Examination of categorical variables
for col in cat_cols:
    cat_summary(df, col)

NEW_AGE_CAT      Ratio
young              488  63.541667
middleage          230  29.947917
old                 50   6.510417
##########################################
   OUTCOME      Ratio
0      500  65.104167
1      268  34.895833
##########################################
             NEW_GLUCOSE_CAT      Ratio
normal                   571  74.348958
prediabetes              197  25.651042
##########################################
             NEW_BMI_RANGE      Ratio
obese                  472  61.458333
overweight             179  23.307292
healty                 102  13.281250
underweight             15   1.953125
##########################################
        NEW_BLOODPRESSURE      Ratio
normal                563  73.307292
hs1                   145  18.880208
hs2                    60   7.812500

# Examination of numerical variables
df[num_cols].describe().T

Standardization

The code provided is used to prepare data for a machine learning model. It is processing a dataframe, “df,” which contains information about patients who have been diagnosed with diabetes. The first step is to standardize the numerical columns in the dataframe using the StandardScaler from the scikit-learn library. Then the target variable, “OUTCOME”, is separated from the rest of the data, and stored in the “y” variable. The input variables, without the outcome, are stored in the “X” variable.

The “diabetes_data_prep” function takes in a dataframe as an argument and performs various data preparation steps on the input data. It starts by converting the column names of the dataframe to uppercase. Then it creates new columns in the dataframe by categorizing some of the numerical columns, such as “GLUCOSE,” “AGE,” “BMI,” and “BLOODPRESSURE.” These categorical columns are then one-hot encoded and the resulting dataframe is standardized. Finally, the function returns the input data, X, and the target variable, y.

# Standardization
X_scaled = StandardScaler().fit_transform(df[num_cols])
df[num_cols] = pd.DataFrame(X_scaled, columns=df[num_cols].columns)

y = df["OUTCOME"]
X = df.drop(["OUTCOME"], axis=1)

check_df(X)

def diabetes_data_prep(dataframe):
    dataframe.columns = [col.upper() for col in dataframe.columns]

    # Glucose
    dataframe['NEW_GLUCOSE_CAT'] = pd.cut(x=dataframe['GLUCOSE'], bins=[-1, 139, 200], labels=["normal", "prediabetes"])

    # Age
    dataframe.loc[(dataframe['AGE'] < 35), "NEW_AGE_CAT"] = 'young'
    dataframe.loc[(dataframe['AGE'] >= 35) & (dataframe['AGE'] <= 55), "NEW_AGE_CAT"] = 'middleage'
    dataframe.loc[(dataframe['AGE'] > 55), "NEW_AGE_CAT"] = 'old'

    # BMI
    dataframe['NEW_BMI_RANGE'] = pd.cut(x=dataframe['BMI'], bins=[-1, 18.5, 24.9, 29.9, 100],
                                        labels=["underweight", "healty", "overweight", "obese"])

    # BloodPressure
    dataframe['NEW_BLOODPRESSURE'] = pd.cut(x=dataframe['BLOODPRESSURE'], bins=[-1, 79, 89, 123],
                                            labels=["normal", "hs1", "hs2"])

    cat_cols, num_cols, cat_but_car = grab_col_names(dataframe, cat_th=5, car_th=20)

    cat_cols = [col for col in cat_cols if "OUTCOME" not in col]

    df = one_hot_encoder(dataframe, cat_cols, drop_first=True)

    df.columns = [col.upper() for col in df.columns]

    cat_cols, num_cols, cat_but_car = grab_col_names(df, cat_th=5, car_th=20)

    cat_cols = [col for col in cat_cols if "OUTCOME" not in col]

    replace_with_thresholds(df, "INSULIN")

    X_scaled = StandardScaler().fit_transform(df[num_cols])
    df[num_cols] = pd.DataFrame(X_scaled, columns=df[num_cols].columns)

    y = df["OUTCOME"]
    X = df.drop(["OUTCOME"], axis=1)

    return X, y

check_df(df)

X, y = diabetes_data_prep(df)

check_df(X)

5. Base Models

def base_models(X, y, scoring="roc_auc"):
    print("Base Models....")
    classifiers = [('LR', LogisticRegression()),
                   ('KNN', KNeighborsClassifier()),
                   ("SVC", SVC()),
                   ("CART", DecisionTreeClassifier()),
                   ("RF", RandomForestClassifier()),
                   ('Adaboost', AdaBoostClassifier()),
                   ('GBM', GradientBoostingClassifier()),
                   ('XGBoost', XGBClassifier(use_label_encoder=False, eval_metric='logloss')),
                   ('LightGBM', LGBMClassifier()),
                   ('CatBoost', CatBoostClassifier(verbose=False))]

    for name, classifier in classifiers:
        cv_results = cross_validate(classifier, X, y, cv=3, scoring=scoring)
        print(f"{scoring}: {round(cv_results['test_score'].mean(), 4)} ({name}) ")

base_models(X, y, scoring="accuracy")

Base Models....
accuracy: 0.7734 (LR) 
accuracy: 0.7409 (KNN) 
accuracy: 0.776 (SVC) 
accuracy: 0.6797 (CART) 
accuracy: 0.7643 (RF) 
accuracy: 0.7695 (Adaboost) 
accuracy: 0.7682 (GBM) 
accuracy: 0.7422 (XGBoost) 
accuracy: 0.7292 (LightGBM) 
accuracy: 0.7539 (CatBoost)

6. Automated Hyperparameter Optimization

knn_params = {"n_neighbors": range(2, 50)}

cart_params = {'max_depth': range(1, 20),
               "min_samples_split": range(2, 30)}

rf_params = {"max_depth": [8, 15, None],
             "max_features": [5, 7, "auto"],
             "min_samples_split": [15, 20],
             "n_estimators": [200, 300]}

xgboost_params = {"learning_rate": [0.1, 0.01],
                  "max_depth": [5, 8],
                  "n_estimators": [100, 200]}

lightgbm_params = {"learning_rate": [0.01, 0.1],
                   "n_estimators": [300, 500]}


classifiers = [('KNN', KNeighborsClassifier(), knn_params),
               ("CART", DecisionTreeClassifier(), cart_params),
               ("RF", RandomForestClassifier(), rf_params),
               ('XGBoost', XGBClassifier(use_label_encoder=False, eval_metric='logloss'), xgboost_params),
               ('LightGBM', LGBMClassifier(), lightgbm_params)]

def hyperparameter_optimization(X, y, cv=3, scoring="roc_auc"):
    print("Hyperparameter Optimization....")
    best_models = {}
    for name, classifier, params in classifiers:
        print(f"########## {name} ##########")
        cv_results = cross_validate(classifier, X, y, cv=cv, scoring=scoring)
        print(f"{scoring} (Before): {round(cv_results['test_score'].mean(), 4)}")

        gs_best = GridSearchCV(classifier, params, cv=cv, n_jobs=-1, verbose=False).fit(X, y)
        final_model = classifier.set_params(**gs_best.best_params_)

        cv_results = cross_validate(final_model, X, y, cv=cv, scoring=scoring)
        print(f"{scoring} (After): {round(cv_results['test_score'].mean(), 4)}")
        print(f"{name} best params: {gs_best.best_params_}", end="\n\n")
        best_models[name] = final_model
    return best_models

best_models = hyperparameter_optimization(X, y)

Hyperparameter Optimization....
########## KNN ##########
roc_auc (Before): 0.7671
roc_auc (After): 0.8026
KNN best params: {'n_neighbors': 17}

########## CART ##########
roc_auc (Before): 0.6392
roc_auc (After): 0.7943
CART best params: {'max_depth': 6, 'min_samples_split': 23}

########## RF ##########
roc_auc (Before): 0.8268
roc_auc (After): 0.8337
RF best params: {'max_depth': 8, 'max_features': 'auto', 'min_samples_split': 20, 'n_estimators': 200}

########## XGBoost ##########
roc_auc (Before): 0.7978
roc_auc (After): 0.8064
XGBoost best params: {'learning_rate': 0.1, 'max_depth': 5, 'n_estimators': 200}

########## LightGBM ##########
roc_auc (Before): 0.8009
roc_auc (After): 0.8139
LightGBM best params: {'learning_rate': 0.01, 'n_estimators': 300}

5. Stacking & Ensemble Learning

Ensemble techniques are the methods that use multiple learning algorithms or models to produce one optimal predictive model. The model produced has better performance than the base learners taken alone.

This code uses the VotingClassifier class from scikit-learn to build an ensemble classifier. It uses the best models found in the hyperparameter optimization step, which are ‘KNN’, ‘RF’, and ‘LightGBM’, as the individual models in the ensemble. The ensemble uses a “soft” voting method, where the predictions of the individual models are averaged. The function then fits the voting classifier on the input features X and target variable y. Finally, the function performs 3-fold cross-validation on the voting classifier and calculates the mean accuracy, F1 score, and ROC-AUC. The results are then printed to the console.

def voting_classifier(best_models, X, y):
    print("Voting Classifier...")

    voting_clf = VotingClassifier(estimators=[('KNN', best_models["KNN"]),
                                              ('RF', best_models["RF"]),
                                              ('LightGBM', best_models["LightGBM"])],
                                  voting='soft').fit(X, y)

    cv_results = cross_validate(voting_clf, X, y, cv=3, scoring=["accuracy", "f1", "roc_auc"])
    print(f"Accuracy: {cv_results['test_accuracy'].mean()}")
    print(f"F1Score: {cv_results['test_f1'].mean()}")
    print(f"ROC_AUC: {cv_results['test_roc_auc'].mean()}")
    return voting_clf

voting_clf = voting_classifier(best_models, X, y)

Voting Classifier...
Accuracy: 0.7643229166666666
F1Score: 0.6162162162162161
ROC_AUC: 0.8291571182800231

7. Prediction for a New Observation

X.columns
random_user = X.sample(1, random_state=45)
voting_clf.predict(random_user)

joblib.dump(voting_clf, "voting_clf2.pkl")

new_model = joblib.load("voting_clf2.pkl")
new_model.predict(random_user)

array([1])

This code first gets the column names of the dataframe X by calling X.columns. Then, it samples a random row from X and assigns it to the variable random_user. The predict method is called on the voting_classifier object voting_clf using the random user data as input. The trained classifier model voting_clf is then saved using the joblib.dump function, with the file name “voting_clf2.pkl”. Finally, the saved model is loaded back into memory using joblib.load and the predict method is called again on the loaded model using the same random user data, which returns the predicted class.

def main():
    df = pd.read_csv("/kaggle/input/pima-indians-diabetes-database/diabetes.csv")
    X, y = diabetes_data_prep(df)
    base_models(X, y)
    best_models = hyperparameter_optimization(X, y)
    voting_clf = voting_classifier(best_models, X, y)
    joblib.dump(voting_clf, "voting_clf.pkl")
    return voting_clf

if __name__ == "__main__":
    print("Process started")
    main()

Process started
Base Models....
roc_auc: 0.8409 (LR) 
roc_auc: 0.791 (KNN) 
roc_auc: 0.8355 (SVC) 
roc_auc: 0.6717 (CART) 
roc_auc: 0.8246 (RF) 
roc_auc: 0.8196 (Adaboost) 
roc_auc: 0.8244 (GBM) 
roc_auc: 0.8015 (XGBoost) 
roc_auc: 0.807 (LightGBM) 
roc_auc: 0.8349 (CatBoost) 
Hyperparameter Optimization....
########## KNN ##########
roc_auc (Before): 0.8191
roc_auc (After): 0.8211
KNN best params: {'n_neighbors': 20}

########## CART ##########
roc_auc (Before): 0.7943
roc_auc (After): 0.7943
CART best params: {'max_depth': 6, 'min_samples_split': 23}

########## RF ##########
roc_auc (Before): 0.8397
roc_auc (After): 0.8351
RF best params: {'max_depth': 8, 'max_features': 5, 'min_samples_split': 15, 'n_estimators': 200}

########## XGBoost ##########
roc_auc (Before): 0.8067
roc_auc (After): 0.8179
XGBoost best params: {'learning_rate': 0.01, 'max_depth': 5, 'n_estimators': 200}

########## LightGBM ##########
roc_auc (Before): 0.8185
roc_auc (After): 0.8185
LightGBM best params: {'learning_rate': 0.01, 'n_estimators': 300}

Voting Classifier...
Accuracy: 0.7708333333333334
F1Score: 0.6371707388656541
ROC_AUC: 0.8366793164869816

The code provided is for training a machine learning model using the Pima Indians Diabetes dataset. The code performs the following steps:

Reads the diabetes dataset from the specified file path into a Pandas dataframe Performs data preparation by extracting the feature matrix (X) and target vector (y) Runs a base model selection process using various machine learning algorithms Performs hyperparameter optimization to select the best models among the base models Creates a voting classifier using the best models, which makes predictions by aggregating the outputs of the individual classifiers Saves the voting classifier to a file using the joblib library Finally, the function returns the trained voting classifier. The code also checks if the script is being executed as the main program, in which case it will run the main function and print the message "Process started".

Conclusion

In conclusion, the comprehensive analysis of Pima Indians Diabetes provides a comprehensive understanding of the relationship between different variables and the incidence of diabetes. The predictive models built in this study can be used to classify patients as diabetic or non-diabetic and to guide the development of more effective treatment and prevention strategies.

You can find all the details of the project in my kaggle notebook.

🔔 If you’d like to get in touch, you can find me on Linkedin, drop me an email or follow me 👇

• Follow me: LinkedIn | Kaggle | Github

References

1. https://www.miuul.com/
2. https://aeturrell.github.io/skimpy/
3. https://github.com/AutoViML/AutoViz

Subscribe to DDIntel Here.

Visit our website here: https://www.datadriveninvestor.com

Join our network here: https://datadriveninvestor.com/collaborate

A Comprehensive Analysis of Pima Indians Diabetes was originally published in DataDrivenInvestor on Medium, where people are continuing the conversation by highlighting and responding to this story.

In this article, we will implement the application on the data set on the Kaggle platform. This competition is provided as a way to explore different time series techniques on a relatively simple and clean dataset. You are given 5 years of store-item sales data, and asked to predict 3 months of sales for 50 different items at 10 different stores.

Dataset Description

The objective of this competition is to predict 3 months of item-level sales data at different store locations.

File descriptions

• train.csv — Training data
• test.csv — Test data (Note: the Public/Private split is time based)
• sample_submission.csv — a sample submission file in the correct format

Data fields

• date — Date of the sale data. There are no holiday effects or store closures.
• store — Store ID
• item — Item ID
• sales — Number of items sold at a particular store on a particular date.

Import Libraries

import time    # To time our operations
import numpy as np  # linear algebra
import pandas as pd     # data processing, CSV file I/O (e.g. pd.read_csv)
from matplotlib import pyplot as plt        # To visualize
import seaborn as sns       # To visualize
import lightgbm as lgb      # To train our model
import warnings     # To ignore some warnings that we don't care
warnings.filterwarnings('ignore')   
pd.set_option('display.max_columns', None)  # Unlimited columns
pd.set_option('display.width', 500)     # Width of the display in the notebook

Loading the data

train = pd.read_csv('/kaggle/input/demand-forecasting-kernels-only/train.csv', parse_dates=['date'])    # parse_dates: Converts the date column to datetime format.
test = pd.read_csv('/kaggle/input/demand-forecasting-kernels-only/test.csv', parse_dates=['date'])
sample_sub = pd.read_csv('/kaggle/input/demand-forecasting-kernels-only/sample_submission.csv')

• We combined train and test data with concat.
• One of the reasons we combine train and test data is to avoid class differences and size differences.

df = pd.concat([train, test], sort=False)

EDA (Exploratory Data Analysis)

def check_df(dataframe, head=5):
    print("##################### Shape #####################")
    print(dataframe.shape)
    print("##################### Types #####################")
    print(dataframe.dtypes)
    print("##################### Head #####################")
    print(dataframe.head(head))
    print("##################### Tail #####################")
    print(dataframe.tail(head))
    print("##################### NA #####################")
    print(dataframe.isnull().sum())
    print("##################### Quantiles #####################")
    print(dataframe.quantile([0, 0.05, 0.50, 0.95, 0.99, 1]).T)

This function is used to perform basic exploratory data analysis (EDA) on a pandas dataframe. The function prints the shape of the dataframe, the data types of the columns, the first and last few rows of the dataframe, the number of missing values (NA) in each column, and the quantiles of the columns. The optional head parameter allows you to specify how many rows to display when printing the first and last few rows of the dataframe. By providing this information, the function helps in understanding the structure, contents and quality of the data.

df["date"].min(), df["date"].max()            # First and last dates in the dataset

(Timestamp('2013-01-01 00:00:00'), Timestamp('2018-03-31 00:00:00'))

check_df(df)

##################### Shape #####################
(958000, 5)
##################### Types #####################
date     datetime64[ns]
store             int64
item              int64
sales           float64
id              float64
dtype: object
##################### Head #####################
        date  store  item  sales  id
0 2013-01-01      1     1   13.0 NaN
1 2013-01-02      1     1   11.0 NaN
2 2013-01-03      1     1   14.0 NaN
3 2013-01-04      1     1   13.0 NaN
4 2013-01-05      1     1   10.0 NaN
##################### Tail #####################
            date  store  item  sales       id
44995 2018-03-27     10    50    NaN  44995.0
44996 2018-03-28     10    50    NaN  44996.0
44997 2018-03-29     10    50    NaN  44997.0
44998 2018-03-30     10    50    NaN  44998.0
44999 2018-03-31     10    50    NaN  44999.0
##################### NA #####################
date          0
store         0
item          0
sales     45000
id       913000
dtype: int64
##################### Quantiles #####################
       0.00     0.05     0.50      0.95      0.99     1.00
store   1.0     1.00      5.5     10.00     10.00     10.0
item    1.0     3.00     25.5     48.00     50.00     50.0
sales   0.0    16.00     47.0    107.00    135.00    231.0
id      0.0  2249.95  22499.5  42749.05  44549.01  44999.0

df[["store"]].nunique()      # Number of stores

store    10
dtype: int64

df[["item"]].nunique()       # Number of products

item    50
dtype: int64

df.groupby(["store"])["item"].nunique()    # Number of products in store breakdown

store
1     50
2     50
3     50
4     50
5     50
6     50
7     50
8     50
9     50
10    50
Name: item, dtype: int64

df.groupby(["store", "item"]).agg({"sales": ["sum"]})     # Total sales in store and product breakdown

df.groupby(["store", "item"]).agg({"sales": ["sum", "mean", "median", "std"]})          # Total sales, average sales, median sales and standard deviation in store and product breakdown

Feature Engineering

def create_date_features(df):
    df['month'] = df.date.dt.month
    df['day_of_month'] = df.date.dt.day
    df['day_of_year'] = df.date.dt.dayofyear
    df['week_of_year'] = df.date.dt.weekofyear
    df['day_of_week'] = df.date.dt.dayofweek
    df['year'] = df.date.dt.year
    df["is_wknd"] = df.date.dt.weekday // 4
    df['is_month_start'] = df.date.dt.is_month_start.astype(int)
    df['is_month_end'] = df.date.dt.is_month_end.astype(int)
    return df

df = create_date_features(df)

This function create_date_features creates various date-related features from the ‘date’ column of a given dataframe df. It returns the modified dataframe with the following additional columns:

• ‘month’: month of the year (1–12)
• ‘day_of_month’: day of the month (1–31)
• ‘day_of_year’: day of the year (1–365 or 1–366 in a leap year)
• ‘week_of_year’: week of the year (1–52)
• ‘day_of_week’: day of the week (0–6, where Monday is 0 and Sunday is 6)
• ‘year’: year
• ‘is_wknd’: a binary variable indicating whether the day is a weekend (0 or 1)
• ‘is_month_start’: a binary variable indicating whether the day is the first day of a month (0 or 1)
• ‘is_month_end’: a binary variable indicating whether the day is the last day of a month (0 or 1)

df.groupby(["store", "item", "month"]).agg({"sales": ["sum", "mean", "median", "std"]})

Autoviz

One of the key features of Autoviz is its ability to handle large and complex datasets. It can quickly identify patterns and relationships within the data, making it an ideal tool for data exploration and discovery. Additionally, Autoviz is highly customizable, allowing users to specify the types of visualizations they want to generate, as well as the color schemes and layouts of the charts.

pip install autoviz

from autoviz.AutoViz_Class import AutoViz_Class
AV = AutoViz_Class()

%matplotlib inline 
# Generate visualizations
AV.AutoViz("/kaggle/input/demand-forecasting-kernels-only/test.csv")

Random Noise

def random_noise(dataframe):
    return np.random.normal(scale=1.6, size=(len(dataframe),))

Lag/Shifted Features

Lag or shifted features are a commonly used feature engineering technique for time series forecasting. Lag features involve shifting the time series data by a specified number of time steps and using the shifted values as input features for the forecasting model.

For example, in a time series forecasting problem, the target variable may be the sales of a product in a given time period. To include information from past time periods, we can create lag features by shifting the sales data by one or more time steps and using the shifted values as input features. This provides the model with additional information about past trends and patterns that may be useful for making predictions.

Lag features can be created for multiple lags, such as a lag of 1, 2, 3, etc. The choice of how many lags to include and which lags to include depends on the characteristics of the time series data and the goals of the forecasting. It is also possible to create lag features for other related variables, such as prices, competitor data, etc.

Including lag features can be a powerful tool for improving the performance of time series forecasting models, but it can also lead to overfitting if not used carefully. It is important to carefully evaluate the impact of lag features on the model performance and to use appropriate methods, such as cross-validation, to avoid overfitting.

df.sort_values(by=['store', 'item', 'date'], axis=0, inplace=True)

pd.DataFrame({"sales": df["sales"].values[0:10],    # first 10 values in sales column
              "lag1": df["sales"].shift(1).values[0:10],    # The previous values of the first 10 values in the sales column
              "lag2": df["sales"].shift(2).values[0:10],    # The two previous values of the first 10 values in the sales column
              "lag3": df["sales"].shift(3).values[0:10],    # three previous values of the first 10 values in the sales column
              "lag4": df["sales"].shift(4).values[0:10]})   # four previous values of the first 10 values in the sales column

df.groupby(["store", "item"])['sales'].head()0         13.0
1         11.0
2         14.0
3         13.0
4         10.0
          ... 
911174    33.0
911175    37.0
911176    46.0
911177    51.0
911178    41.0
Name: sales, Length: 2500, dtype: float64

0         13.0
1         11.0
2         14.0
3         13.0
4         10.0
          ... 
911174    33.0
911175    37.0
911176    46.0
911177    51.0
911178    41.0
Name: sales, Length: 2500, dtype: float64

df.groupby(["store", "item"])['sales'].transform(lambda x: x.shift(1))      # Replaces the values in the sales column with the previous values.

0         NaN
1        13.0
2        11.0
3        14.0
4        13.0
         ... 
44995     NaN
44996     NaN
44997     NaN
44998     NaN
44999     NaN
Name: sales, Length: 958000, dtype: float64

def lag_features(dataframe, lags):
    for lag in lags:
        dataframe['sales_lag_' + str(lag)] = dataframe.groupby(["store", "item"])['sales'].transform(
            lambda x: x.shift(lag)) + random_noise(dataframe)
    return dataframe

df = lag_features(df, [91, 98, 105, 112, 119, 126, 182, 364, 546, 728])

check_df(df)

The lag_features function is a feature engineering function that creates lag features for a given time series dataframe. The function takes in a dataframe and a list of lags, and for each lag in the list, it creates a new feature in the dataframe that is the sales data shifted by the specified number of time steps. The new feature is created using the shift function in pandas, which shifts the values in the sales column by the specified number of time steps. The shifted values are then grouped by the store and item columns using the groupby function and transformed using a lambda function that adds random noise to the shifted values. The check_df function is used to inspect the resulting dataframe to ensure that the lag features have been created correctly. The exact purpose of this function is not clear from the code provided, but it may be used to check the shape, data types, or values in the dataframe.

Rolling Mean Features

Rolling mean features are another commonly used feature engineering technique for time series forecasting. Rolling mean features involve calculating the mean of the time series data over a specified window size and using the mean values as input features for the forecasting model. The window size can be specified in terms of the number of time steps or in terms of a time span, such as a number of days or weeks. The choice of window size depends on the characteristics of the time series data and the goals of the forecasting.

Including rolling mean features can be a useful tool for improving the performance of time series forecasting models, but it is important to use appropriate methods to avoid overfitting. It is also important to carefully evaluate the impact of the rolling mean feature on the model performance and to consider alternative feature engineering techniques, such as lag features, as well.

pd.DataFrame({"sales": df["sales"].values[0:10],
              "roll2": df["sales"].rolling(window=2).mean().values[0:10],
              "roll3": df["sales"].rolling(window=3).mean().values[0:10],
              "roll5": df["sales"].rolling(window=5).mean().values[0:10]})

pd.DataFrame({"sales": df["sales"].values[0:10],
              "roll2": df["sales"].shift(1).rolling(window=2).mean().values[0:10],
              "roll3": df["sales"].shift(1).rolling(window=3).mean().values[0:10],
              "roll5": df["sales"].shift(1).rolling(window=5).mean().values[0:10]})


def roll_mean_features(dataframe, windows):
    for window in windows:
        dataframe['sales_roll_mean_' + str(window)] = dataframe.groupby(["store", "item"])['sales']. \
                                                          transform(
            lambda x: x.shift(1).rolling(window=window, min_periods=10, win_type="triang").mean()) + random_noise(
            dataframe)
    return dataframe


df = roll_mean_features(df, [365, 546])

Exponentially Weighted Mean Features

Exponentially weighted mean features are similar to rolling mean features, but instead of using a fixed window size to calculate the mean, they use an exponentially decaying weighting scheme to give more importance to recent observations and less importance to older observations.

The decay factor in the weighting scheme can be specified to control the importance given to recent observations. A smaller decay factor will give more importance to recent observations, while a larger decay factor will give more importance to older observations. The choice of decay factor depends on the characteristics of the time series data and the goals of the forecasting.

Including exponentially weighted mean features can be a useful tool for improving the performance of time series forecasting models, but it is important to use appropriate methods to avoid overfitting. It is also important to carefully evaluate the impact of the exponentially weighted mean feature on the model performance and to consider alternative feature engineering techniques, such as lag features or rolling mean features, as well.

pd.DataFrame({"sales": df["sales"].values[0:10],
              "roll2": df["sales"].shift(1).rolling(window=2).mean().values[0:10],
              "ewm099": df["sales"].shift(1).ewm(alpha=0.99).mean().values[0:10],
              "ewm095": df["sales"].shift(1).ewm(alpha=0.95).mean().values[0:10],
              "ewm07": df["sales"].shift(1).ewm(alpha=0.7).mean().values[0:10],
              "ewm02": df["sales"].shift(1).ewm(alpha=0.1).mean().values[0:10]})

def ewm_features(dataframe, alphas, lags):
    for alpha in alphas:
        for lag in lags:
            dataframe['sales_ewm_alpha_' + str(alpha).replace(".", "") + "_lag_" + str(lag)] = \
                dataframe.groupby(["store", "item"])['sales'].transform(lambda x: x.shift(lag).ewm(alpha=alpha).mean())
    return dataframe

alphas = [0.95, 0.9, 0.8, 0.7, 0.5]
lags = [91, 98, 105, 112, 180, 270, 365, 546, 728]

df = ewm_features(df, alphas, lags)
check_df(df)

One-Hot Encoding

The code uses the pd.get_dummies function from the pandas library to perform one-hot encoding on several columns in the dataframe df. One-hot encoding is a common preprocessing step in machine learning, and it is used to convert categorical variables into numerical representations that can be used as input features for models.

The columns parameter of pd.get_dummies is used to specify which columns in the dataframe should be one-hot encoded. In this case, the store, item, day_of_week, and month columns are selected for one-hot encoding.

One-hot encoding creates a new binary column for each unique category in the specified columns. For example, if the store column has three unique categories, then three new columns will be created, one for each category, with a value of 1 or 0 indicating whether or not that category is present in the data.

The purpose of one-hot encoding is to allow categorical variables to be used as input features for machine learning models, which typically expect numerical inputs. By one-hot encoding the categorical variables, the information contained in the categories is transformed into a numerical representation that can be used by the models.

The check_df function is used to inspect the resulting dataframe to ensure that the one-hot encoding has been applied correctly. The exact purpose of this function is not clear from the code provided, but it may be used to check the shape, data types, or values in the dataframe.

df = pd.get_dummies(df, columns=['store', 'item', 'day_of_week', 'month'])
check_df(df)

Converting sales to log(1+sales)

df['sales'] = np.log1p(df["sales"].values)
check_df(df)

The code above applies a log transformation to the sales data in the dataframe df. The log transformation is performed using the np.log1p function from the NumPy library, which calculates the logarithm of (1 + x) for an input array x.

The log transformation is a common preprocessing step for time series data, as it can help to stabilize the variance in the data and make it more suitable for modeling. By transforming the sales data into the log-space, the influence of large fluctuations in the sales data is reduced, and the relationships between the features and the target variable become clearer.

Model

# MAE: mean absolute error
# MAPE: mean absolute percentage error
# SMAPE: Symmetric mean absolute percentage error (adjusted MAPE)


def smape(preds, target):
    n = len(preds)
    masked_arr = ~((preds == 0) & (target == 0))
    preds, target = preds[masked_arr], target[masked_arr]
    num = np.abs(preds - target)
    denom = np.abs(preds) + np.abs(target)
    smape_val = (200 * np.sum(num / denom)) / n
    return smape_val

def lgbm_smape(preds, train_data):
    labels = train_data.get_label()
    smape_val = smape(np.expm1(preds), np.expm1(labels))
    return 'SMAPE', smape_val, False

The code above defines a custom metric for evaluating the performance of machine learning models, called the symmetric mean absolute percentage error (SMAPE). The SMAPE metric is used to measure the percentage error between predicted values and actual target values.

The smape function takes two arrays as input, preds and target, which represent the predicted values and actual target values, respectively. The function first checks for cases where both the predicted value and the target value are zero and masks these cases from further calculation.

Next, the absolute difference between the predicted and target values is calculated, and this is divided by the sum of the absolute values of the predicted and target values. The SMAPE value is calculated as the average of these ratios across all data points. The SMAPE value is expressed as a percentage, with a value of 0 indicating a perfect prediction and a value greater than 0 indicating an error.

The SMAPE metric is a symmetric metric, meaning that it treats over-prediction and under-prediction equally. This makes it useful for evaluating the performance of models on data that may have both positive and negative values, as well as for data with a skewed distribution, where one type of prediction error is more common than the other.

Time-Based Validation Sets

train
test

# Train set until the beginning of 2017 (until the end of 2016).
train = df.loc[(df["date"] < "2017-01-01"), :]
# First 3 months of 2017 validation set.
val = df.loc[(df["date"] >= "2017-01-01") & (df["date"] < "2017-04-01"), :]

cols = [col for col in train.columns if col not in ['date', 'id', "sales", "year"]]

Y_train = train['sales']
X_train = train[cols]

Y_val = val['sales']
X_val = val[cols]

Y_train.shape, X_train.shape, Y_val.shape, X_val.shape

((730500,), (730500, 142), (45000,), (45000, 142))

Time Series Model with LightGBM

# LightGBM parameters
lgb_params = {'num_leaves': 10,
              'learning_rate': 0.02,
              'feature_fraction': 0.8,
              'max_depth': 5,
              'verbose': 0,
              'num_boost_round': 1000,
              'early_stopping_rounds': 200,
              'nthread': -1}

# metric mae: l1, absolute loss, mean_absolute_error, regression_l1
# mse: l2, square loss, mean_squared_error, mse, regression_l2, regression
# rmse, root square loss, root_mean_squared_error, l2_root
# mape, MAPE loss, mean_absolute_percentage_error
lgbtrain = lgb.Dataset(data=X_train, label=Y_train, feature_name=cols)
lgbval = lgb.Dataset(data=X_val, label=Y_val, reference=lgbtrain, feature_name=cols)

model = lgb.train(lgb_params, lgbtrain,
                  valid_sets=[lgbtrain, lgbval],
                  num_boost_round=lgb_params['num_boost_round'],
                  early_stopping_rounds=lgb_params['early_stopping_rounds'],
                  feval=lgbm_smape,
                  verbose_eval=100)

[LightGBM] [Warning] Auto-choosing col-wise multi-threading, the overhead of testing was 0.371499 seconds.
You can set `force_col_wise=true` to remove the overhead.
Training until validation scores don't improve for 200 rounds
[100]	training's l2: 0.0512913	training's SMAPE: 17.5863	valid_1's l2: 0.0529452	valid_1's SMAPE: 17.4429
[200]	training's l2: 0.0352035	training's SMAPE: 14.5662	valid_1's l2: 0.0369583	valid_1's SMAPE: 14.8761
[300]	training's l2: 0.0324151	training's SMAPE: 14.0096	valid_1's l2: 0.0341202	valid_1's SMAPE: 14.4213
[400]	training's l2: 0.0313431	training's SMAPE: 13.7968	valid_1's l2: 0.0332429	valid_1's SMAPE: 14.2604
[500]	training's l2: 0.0306813	training's SMAPE: 13.6622	valid_1's l2: 0.0326193	valid_1's SMAPE: 14.1237
[600]	training's l2: 0.0302368	training's SMAPE: 13.5684	valid_1's l2: 0.032161	valid_1's SMAPE: 14.0164
[700]	training's l2: 0.0298809	training's SMAPE: 13.4971	valid_1's l2: 0.031855	valid_1's SMAPE: 13.9461
[800]	training's l2: 0.0295966	training's SMAPE: 13.4403	valid_1's l2: 0.0316245	valid_1's SMAPE: 13.8914
[900]	training's l2: 0.0293529	training's SMAPE: 13.3908	valid_1's l2: 0.0314269	valid_1's SMAPE: 13.8454
[1000]	training's l2: 0.0291449	training's SMAPE: 13.3487	valid_1's l2: 0.031298	valid_1's SMAPE: 13.8149
Did not meet early stopping. Best iteration is:
[1000]	training's l2: 0.0291449	training's SMAPE: 13.3487	valid_1's l2: 0.031298	valid_1's SMAPE: 13.8149

y_pred_val = model.predict(X_val, num_iteration=model.best_iteration)
smape(np.expm1(y_pred_val), np.expm1(Y_val))

13.814874270278832

Feature Importance Levels

def plot_lgb_importances(model, plot=False, num=10):
    gain = model.feature_importance('gain')
    feat_imp = pd.DataFrame({'feature': model.feature_name(),
                             'split': model.feature_importance('split'),
                             'gain': 100 * gain / gain.sum()}).sort_values('gain', ascending=False)
    if plot:
        plt.figure(figsize=(10, 10))
        sns.set(font_scale=1)
        sns.barplot(x="gain", y="feature", data=feat_imp[0:25])
        plt.title('feature')
        plt.tight_layout()
        plt.show()
    else:
        print(feat_imp.head(num))
    return feat_imp

plot_lgb_importances(model, num=200)

                feature          split       gain
17          sales_roll_mean_546    936  54.465207
13                sales_lag_364   1271  13.219649
16          sales_roll_mean_365    606   9.754873
60   sales_ewm_alpha_05_lag_365    384   4.890732
18   sales_ewm_alpha_095_lag_91     81   2.518718
..                          ...    ...        ...
66                      store_4      0   0.000000
83                      item_11      0   0.000000
111                     item_39      0   0.000000
74                       item_2      0   0.000000
71                      store_9      0   0.000000

[142 rows x 3 columns]

plot_lgb_importances(model, num=30, plot=True)

importance_zero = feat_imp[feat_imp["gain"] == 0]["feature"].values
imp_feats = [col for col in cols if col not in importance_zero]
len(imp_feats)

110

Final Model

train = df.loc[~df.sales.isna()]
Y_train = train['sales']
X_train = train[cols]


test = df.loc[df.sales.isna()]
X_test = test[cols]

lgb_params = {'num_leaves': 10,
              'learning_rate': 0.02,
              'feature_fraction': 0.8,
              'max_depth': 5,
              'verbose': 0,
              'nthread': -1,
              "num_boost_round": model.best_iteration}

lgbtrain_all = lgb.Dataset(data=X_train, label=Y_train, feature_name=cols)

final_model = lgb.train(lgb_params, lgbtrain_all, num_boost_round=model.best_iteration)


test_preds = final_model.predict(X_test, num_iteration=model.best_iteration)

Submission File

Submissions are evaluated on SMAPE between forecasts and actual values. We define SMAPE = 0 when the actual and predicted values are both 0.

submission_df = test.loc[:, ["id", "sales"]]
submission_df['sales'] = np.expm1(test_preds)

submission_df['id'] = submission_df.id.astype(int)

submission_df.to_csv("submission.csv", index=False)

🔔 If you’d like to get in touch, you can find me on Linkedin, drop me an email or follow me 👇

• Follow me: LinkedIn | Kaggle | Github

References

1. https://www.miuul.com/zaman-serileri
2. https://github.com/AutoViML/AutoViz
3. https://www.kaggle.com/code/emineyetm/using-machine-learning-for-time-series-forecasting

Subscribe to DDIntel Here.

Visit our website here: https://www.datadriveninvestor.com

Join our network here: https://datadriveninvestor.com/collaborate

Using Machine Learning for Time Series Forecasting was originally published in DataDrivenInvestor on Medium, where people are continuing the conversation by highlighting and responding to this story.

Statistical models are one of the most common approaches used in time series forecasting. These models are mathematical representations of the underlying relationships between the variables in a time series and are used to make predictions about future values. There are several types of statistical models that can be used for time series forecasting, including:

1.Autoregressive (AR) Model: The autoregressive model is another common method for time series forecasting. It is based on the assumption that the future value of a time series is a linear combination of past values and a white noise error term. The number of past values used in the calculation is referred to as the “order” of the model. The larger the order, the more information is used in the prediction, but the more complicated the model becomes.

2.Moving Average (MA) Model: The moving average model is a simple and straightforward method for forecasting time series data. It involves calculating the average of a set of past observations and using that average as the prediction for future values. The number of past observations used in the calculation is referred to as the “window size”. The larger the window size, the smoother the forecast, but the less responsive it will be to changes in the underlying trends.

3.ARMA (Auto-Regressive Moving Average): The ARMA is one of the most commonly used methods to model univariate time series. ARMA(p,q) combines two components: AR(p), and MA(q).

4.Autoregressive Integrated Moving Average (ARIMA) Model: The autoregressive integrated moving average (ARIMA) model combines the moving average and autoregressive models to produce a more sophisticated time series forecasting method. It can handle both linear and non-linear relationships between the variables in a time series, making it a powerful tool for time series analysis. The order of the autoregression and moving average components of the model, as well as the order of integration, are parameters that must be specified in order to fit the model to the data

# autoregressive: p
# integrated: d
# moving average: q
# order: (p, d, q)  # p: AR, d: I, q: MA

5. Seasonal Autoregressive Integrated Moving Average (SARIMA) Model: The seasonal autoregressive integrated moving average (SARIMA) model is an extension of the ARIMA model that takes into account seasonality in the data. Seasonality refers to patterns in the data that repeat at regular intervals, such as daily, weekly, or yearly patterns. The SARIMA model adds an additional set of parameters that capture the seasonal patterns in the data, allowing for more accurate forecasting of time series with strong seasonality.

Let’s practice in python.

Import Libraries & Load Data

import numpy as np
import pandas as pd
import itertools
import matplotlib.pyplot as plt
import statsmodels.api as sm
from sklearn.metrics import mean_absolute_error

from statsmodels.tsa.api import SARIMAX             # SARIMA: Seasonal Autoregressive Integrated Moving Average
from statsmodels.tsa.arima.model import ARIMA       # ARIMA: Autoregressive Integrated Moving Average

import warnings          # warnings: to avoid seeing error messages
warnings.filterwarnings('ignore')

• Atmospheric CO2 from Continuous Air Samples at Mauna Loa Observatory, Hawaii, U.S.A.
• Period of Record: March 1958 — December 2001

data = sm.datasets.co2.load_pandas()    # loads the CO2 dataset into a pandas DataFrame object, which is stored in the data variable.
y = data.data           # accesses the data in the loaded CO2 dataset
y = y['co2'].resample('MS').mean()
y = y.fillna(y.bfill())
train = y[:'1997-12-01']
test = y['1998-01-01':]

• y = y[‘co2’].resample(‘MS’).mean() resamples the CO2 data to obtain the mean value for each month. The argument ‘MS’ stands for “month start”, and indicates that the data should be resampled to obtain the mean value for each month, starting at the beginning of the month.
• y = y.fillna(y.bfill()) fills any missing values in the data with the values from the previous observation. The method bfill() stands for “backward fill”, which means that the missing values are filled with the values from the previous observation.
• train = y[:’1997–12–01'] and test = y[‘1998–01–01’:] split the resampled and cleaned CO2 data into two separate datasets: the training dataset and the test dataset. The training dataset consists of all the data up to and including December 1st, 1997, while the test dataset consists of all the data starting from January 1st, 1998. The training dataset is used to fit the time series model, while the test dataset is used to evaluate the performance of the model.

train.head()

1958-03-01    316.100000
1958-04-01    317.200000
1958-05-01    317.433333
1958-06-01    315.625000
1958-07-01    315.625000
Freq: MS, Name: co2, dtype: float64

ARIMA(Autoregressive Integrated Moving Average)

The order of the ARIMA model is specified as (p, d, q), and the model is fit to the time series data to make predictions about future values based on the historical patterns in the data.

• autoregressive: p
• integrated: d
• moving average: q
• order: (p, d, q) # p: AR, d: I, q: MA

arima_model = ARIMA(train, order=(1, 1, 1)).fit()
arima_model.summary()

• Dep. Variable: The dependent variable that the model is trying to predict. In this case, it is "co2".
• No. Observations: The number of observations in the data used to fit the model.
• Model: The type of ARIMA model fit, including the parameters used. In this case, the model is ARIMA(1, 1, 1), meaning it has 1 autoregressive term, 1 degree of differencing, and 1 moving average component.
• Log Likelihood: The log-likelihood of the model, which is a measure of how well the model fits the data.
• Date and Time: The date and time when the model was fit.
• Sample: The range of the sample data used to fit the model.
• AIC (Akaike Information Criteria): A measure of the relative quality of the model fit, taking into account both the goodness of fit and the number of parameters used. A lower AIC value indicates a better model fit.
• BIC (Bayesian Information Criteria): Similar to AIC, but with a stronger penalty for the number of parameters used.
• HQIC (Hannan-Quinn Information Criteria): Similar to AIC and BIC, but with a different weighting between goodness of fit and complexity.
• Covariance Type: The type of covariance matrix used to estimate the parameters in the model. In this case, it is "opg", which stands for "outer product of the gradients".
• coef: The coefficient of each term in the model.
• std err: The standard error of the coefficient estimate.
• z: The z-score for each coefficient, which measures the number of standard deviations the coefficient is from zero.
• P>|z|: The p-value for each coefficient, which measures the probability of observing a coefficient as extreme as the one estimated, given that the true coefficient is zero.
• [0.025 0.975]: The 95% confidence interval for each coefficient, which provides an interval estimate of the true value of the coefficient.

From these results, we can see that both the autoregressive term ar.L1 and the moving average term ma.L1 are significant, as indicated by the low p-values and the z-scores that are far from zero. The confidence intervals for these coefficients also do not contain zero, which further supports their significance. The sigma2 term is the variance of the residuals (the differences between the observed values and the values predicted by the model) and its low p-value and significant z-score also support the validity of the model fit.

• Ljung-Box (L1) (Q): The value of the Ljung-Box test for residual autocorrelation, with L1 indicating that the test is performed at lag 1. The test statistic measures whether the residuals are independently distributed, with a low value indicating that the residuals are independent.
• Prob(Q): The p-value of the Ljung-Box test, which measures the probability of observing the test statistic as extreme as the one estimated, given that the residuals are independently distributed.
• Jarque-Bera (JB): The value of the Jarque-Bera test, which tests the residuals for normality. A low value indicates that the residuals are close to normally distributed.
• Prob(JB): The p-value of the Jarque-Bera test, which measures the probability of observing the test statistic as extreme as the one estimated, given that the residuals are normally distributed.
• Heteroskedasticity (H): The value of the Breusch-Pagan test, which tests the residuals for heteroskedasticity (unequal variance). A low value indicates that the residuals have equal variance.
• Prob(H) (two-sided): The p-value of the Breusch-Pagan test, which measures the probability of observing the test statistic as extreme as the one estimated, given that the residuals have equal variance.
• Skew: The skewness of the residuals, which measures the degree of asymmetry in the residual distribution. A value close to zero indicates a symmetric distribution.
• Kurtosis: The kurtosis of the residuals, which measures the "peakedness" of the residual distribution. A value close to 3 indicates a normal distribution.

In summary, these results suggest that the ARIMA(1, 1, 1) model is a good fit for the “co2” time series data, with both the autoregressive and moving average terms having significant impact on the predictions.

y_pred = arima_model.forecast(48)[0]       # 48: timeframe to estimate (size of test set) #0: predicted value 1.index
y_pred

365.5586389302753

# visualization function expects input in form of pandas series or dataframe. For that we need to change it.
y_pred = pd.Series(y_pred, index=test.index)

def plot_co2(train, test, y_pred, title):
    mae = mean_absolute_error(test, y_pred)
    train["1985":].plot(legend=True, label="TRAIN", title=f"{title}, MAE: {round(mae,2)}")
    test.plot(legend=True, label="TEST", figsize=(6, 4))
    y_pred.plot(legend=True, label="PREDICTION")
    plt.show()

plot_co2(train, test, y_pred, "ARIMA")

SARIMA(Seasonal Autoregressive Integrated Moving-Average):

SARIMA stands for Seasonal-ARIMA and it includes seasonality contribution to the forecast. The importance of seasonality is quite evident and ARIMA fails to encapsulate that information implicitly.

The Autoregressive (AR), Integrated (I), and Moving Average (MA) parts of the model remain as that of ARIMA. The addition of Seasonality adds robustness to the SARIMA model. It’s represented as:

• order = (p, d, q) They are trend components from ARIMA.
• seasonal_order = (P, D, Q, m) They are seasonal components of SARIMA.

model = SARIMAX(train, order=(1, 0, 1), seasonal_order=(0, 0, 0, 12))

sarima_model = model.fit(disp=0)
sarima_model.summary()

y_pred_test = sarima_model.get_forecast(steps=48)   # get_forecast: timeframe to be estimated
y_pred_test.predicted_mean

y_pred = y_pred_test.predicted_mean
y_pred = pd.Series(y_pred, index=test.index)
plot_co2(train, test, y_pred, "SARIMA")

Hyperparameter Optimization (Determining Model Degrees)

The function sarima_optimizer_mae takes in the training time series data (train), along with grid search parameters for SARIMA models, pdq and seasonal_pdq. The function iterates over all possible combinations of pdq and seasonal_pdq parameters to fit SARIMA models. For each model, it generates predictions using the get_forecast method and calculates the mean absolute error (MAE) between the predictions and the test set. The function keeps track of the SARIMA model with the lowest MAE and returns the corresponding parameters as the best SARIMA model. Finally, the best SARIMA model is fit on the training set and predictions are generated for the test set. The predictions and actual test set values are then plotted using the plot_co2 function.

p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]


def sarima_optimizer_aic(train, pdq, seasonal_pdq):
    best_aic, best_order, best_seasonal_order = float("inf"), None, None
    for param in pdq:
        for param_seasonal in seasonal_pdq:
            try:
                sarimax_model = SARIMAX(train, order=param, seasonal_order=param_seasonal)
                results = sarimax_model.fit(disp=0)
                aic = results.aic
                if aic < best_aic:
                    best_aic, best_order, best_seasonal_order = aic, param, param_seasonal
                print('SARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, aic))
            except:
                continue
    print('SARIMA{}x{}12 - AIC:{}'.format(best_order, best_seasonal_order, best_aic))
    return best_order, best_seasonal_order

best_order, best_seasonal_order = sarima_optimizer_aic(train, pdq, seasonal_pdq)

SARIMA(0, 0, 0)x(0, 0, 0, 12)12 - AIC:6922.606550808509
SARIMA(0, 0, 0)x(0, 0, 1, 12)12 - AIC:6318.305219679214
SARIMA(0, 0, 0)x(0, 1, 0, 12)12 - AIC:1643.9996657880379
SARIMA(0, 0, 0)x(0, 1, 1, 12)12 - AIC:1428.2938925348508
SARIMA(0, 0, 0)x(1, 0, 0, 12)12 - AIC:1820.7236583825966
SARIMA(0, 0, 0)x(1, 0, 1, 12)12 - AIC:1653.729484027914
SARIMA(0, 0, 0)x(1, 1, 0, 12)12 - AIC:1205.025400971046
SARIMA(0, 0, 0)x(1, 1, 1, 12)12 - AIC:973.207727675453
SARIMA(0, 0, 1)x(0, 0, 0, 12)12 - AIC:6269.190364756832
SARIMA(0, 0, 1)x(0, 0, 1, 12)12 - AIC:5663.567928540047
SARIMA(0, 0, 1)x(0, 1, 0, 12)12 - AIC:1229.8447516025176
SARIMA(0, 0, 1)x(0, 1, 1, 12)12 - AIC:1126.5843176793833
SARIMA(0, 0, 1)x(1, 0, 0, 12)12 - AIC:1473.0898881947135
SARIMA(0, 0, 1)x(1, 0, 1, 12)12 - AIC:1414.49020595465
SARIMA(0, 0, 1)x(1, 1, 0, 12)12 - AIC:1008.160268516428
SARIMA(0, 0, 1)x(1, 1, 1, 12)12 - AIC:767.6332844686492
SARIMA(0, 1, 0)x(0, 0, 0, 12)12 - AIC:1523.9134546907358
SARIMA(0, 1, 0)x(0, 0, 1, 12)12 - AIC:1154.743619904219
SARIMA(0, 1, 0)x(0, 1, 0, 12)12 - AIC:585.8949540598103
SARIMA(0, 1, 0)x(0, 1, 1, 12)12 - AIC:371.57825689147813
SARIMA(0, 1, 0)x(1, 0, 0, 12)12 - AIC:608.9685497853327
SARIMA(0, 1, 0)x(1, 0, 1, 12)12 - AIC:410.71468501006075
SARIMA(0, 1, 0)x(1, 1, 0, 12)12 - AIC:495.3771447981145
SARIMA(0, 1, 0)x(1, 1, 1, 12)12 - AIC:373.1539006474194
SARIMA(0, 1, 1)x(0, 0, 0, 12)12 - AIC:1254.7914169695318
SARIMA(0, 1, 1)x(0, 0, 1, 12)12 - AIC:1039.3952573257652
SARIMA(0, 1, 1)x(0, 1, 0, 12)12 - AIC:543.1399633696033
SARIMA(0, 1, 1)x(0, 1, 1, 12)12 - AIC:393.8943760680236
SARIMA(0, 1, 1)x(1, 0, 0, 12)12 - AIC:579.922934531628
SARIMA(0, 1, 1)x(1, 0, 1, 12)12 - AIC:383.09871119803915
SARIMA(0, 1, 1)x(1, 1, 0, 12)12 - AIC:452.11811365175356
SARIMA(0, 1, 1)x(1, 1, 1, 12)12 - AIC:339.1641111208867
SARIMA(1, 0, 0)x(0, 0, 0, 12)12 - AIC:1540.4140252670818
SARIMA(1, 0, 0)x(0, 0, 1, 12)12 - AIC:1171.2466461793142
SARIMA(1, 0, 0)x(0, 1, 0, 12)12 - AIC:579.208370937031
SARIMA(1, 0, 0)x(0, 1, 1, 12)12 - AIC:376.830651073989
SARIMA(1, 0, 0)x(1, 0, 0, 12)12 - AIC:625.4609723210033
SARIMA(1, 0, 0)x(1, 0, 1, 12)12 - AIC:426.6332821688898
SARIMA(1, 0, 0)x(1, 1, 0, 12)12 - AIC:495.7546159387019
SARIMA(1, 0, 0)x(1, 1, 1, 12)12 - AIC:378.4090442279553
SARIMA(1, 0, 1)x(0, 0, 0, 12)12 - AIC:1271.2934125926395
SARIMA(1, 0, 1)x(0, 0, 1, 12)12 - AIC:1055.8968573894322
SARIMA(1, 0, 1)x(0, 1, 0, 12)12 - AIC:544.036979839994
SARIMA(1, 0, 1)x(0, 1, 1, 12)12 - AIC:342.59441296560226
SARIMA(1, 0, 1)x(1, 0, 0, 12)12 - AIC:596.6333264561456
SARIMA(1, 0, 1)x(1, 0, 1, 12)12 - AIC:400.17578732713
SARIMA(1, 0, 1)x(1, 1, 0, 12)12 - AIC:455.6579852467138
SARIMA(1, 0, 1)x(1, 1, 1, 12)12 - AIC:344.40946539807214
SARIMA(1, 1, 0)x(0, 0, 0, 12)12 - AIC:1207.7849901797972
SARIMA(1, 1, 0)x(0, 0, 1, 12)12 - AIC:1004.4872257540137
SARIMA(1, 1, 0)x(0, 1, 0, 12)12 - AIC:556.2117153544921
SARIMA(1, 1, 0)x(0, 1, 1, 12)12 - AIC:347.6673054338916
SARIMA(1, 1, 0)x(1, 0, 0, 12)12 - AIC:588.1986879134922
SARIMA(1, 1, 0)x(1, 0, 1, 12)12 - AIC:390.8372874271904
SARIMA(1, 1, 0)x(1, 1, 0, 12)12 - AIC:463.3371366789647
SARIMA(1, 1, 0)x(1, 1, 1, 12)12 - AIC:349.441132479151
SARIMA(1, 1, 1)x(0, 0, 0, 12)12 - AIC:1157.5365027891548
SARIMA(1, 1, 1)x(0, 0, 1, 12)12 - AIC:999.6119471633762
SARIMA(1, 1, 1)x(0, 1, 0, 12)12 - AIC:531.0882141113823
SARIMA(1, 1, 1)x(0, 1, 1, 12)12 - AIC:323.3128816560523
SARIMA(1, 1, 1)x(1, 0, 0, 12)12 - AIC:560.3099483839314
SARIMA(1, 1, 1)x(1, 0, 1, 12)12 - AIC:367.09243865603963
SARIMA(1, 1, 1)x(1, 1, 0, 12)12 - AIC:442.4874503875827
SARIMA(1, 1, 1)x(1, 1, 1, 12)12 - AIC:324.9964687919045
SARIMA(1, 1, 1)x(0, 1, 1, 12)12 - AIC:323.3128816560523

Final Model

model = SARIMAX(train, order=best_order, seasonal_order=best_seasonal_order)
sarima_final_model = model.fit(disp=0)
y_pred_test = sarima_final_model.get_forecast(steps=48)

y_pred = y_pred_test.predicted_mean
y_pred = pd.Series(y_pred, index=test.index)
plot_co2(train, test, y_pred, "SARIMA")

SARIMA Optimization by MAE

p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]


def sarima_optimizer_mae(train, pdq, seasonal_pdq):
    best_mae, best_order, best_seasonal_order = float("inf"), None, None
    for param in pdq:
        for param_seasonal in seasonal_pdq:
            try:
                model = SARIMAX(train, order=param, seasonal_order=param_seasonal)
                sarima_model = model.fit(disp=0)
                y_pred_test = sarima_model.get_forecast(steps=48)
                y_pred = y_pred_test.predicted_mean
                mae = mean_absolute_error(test, y_pred)
                if mae < best_mae:
                    best_mae, best_order, best_seasonal_order = mae, param, param_seasonal
                print('SARIMA{}x{}12 - MAE:{}'.format(param, param_seasonal, mae))
            except:
                continue
    print('SARIMA{}x{}12 - MAE:{}'.format(best_order, best_seasonal_order, best_mae))
    return best_order, best_seasonal_order

best_order, best_seasonal_order = sarima_optimizer_mae(train, pdq, seasonal_pdq)

model = SARIMAX(train, order=best_order, seasonal_order=best_seasonal_order)
sarima_final_model = model.fit(disp=0)

y_pred_test = sarima_final_model.get_forecast(steps=48)
y_pred = y_pred_test.predicted_mean
y_pred = pd.Series(y_pred, index=test.index)

plot_co2(train, test, y_pred, "SARIMA")

model = SARIMAX(train, order=best_order, seasonal_order=best_seasonal_order)
sarima_final_model = model.fit(disp=0)

y_pred_test = sarima_final_model.get_forecast(steps=48)

# MAE
y_pred = y_pred_test.predicted_mean
y_pred = pd.Series(y_pred, index=test.index)

plot_co2(train, test, y_pred, "SARIMA")
plt.style.use('ggplot')
sarima_final_model.plot_diagnostics(figsize=(15, 10))
plt.show()

• Histogram plus estimated density plot: The red KDE line follows closely with the N(0,1) line. This is a good indication that the residuals are normally distributed.
• The Q-Q-plot: Shows that the ordered distribution of residuals (blue dots) follows the linear trend of the samples taken from a standard normal distribution with N(0, 1). This is an indication that the residuals are normally distributed.
• The standardize residual plot: The residuals over time don’t display any obvious seasonality and appear to be white noise.
• The Correlogram plot: Shows that the time series residuals have low correlation with lagged versions of itself.

Conclusion

In conclusion, time series forecasting with statistical models can be a powerful tool for making predictions about future events based on historical data. This notebook has shown how different statistical models such as ARIMA and SARIMA can be applied to time series data to make accurate forecasts. The results of the analysis demonstrate that while these models have their limitations, they can provide meaningful insights and predictions when used correctly. The appropriate choice of model will depend on the specific requirements of the problem, such as the nature of the underlying data, the size of the time series, and the desired level of accuracy. It is important to carefully consider these factors and experiment with different models to find the best solution for each specific use case.

🔔 If you’d like to get in touch, you can find me on Linkedin, drop me an email or follow me 👇

• Follow me: LinkedIn | Kaggle | Github

References

1. https://www.kaggle.com/code/emineyetm/time-series-forecasting-with-statistical-models
2. https://www.miuul.com/zaman-serileri

Subscribe to DDIntel Here.

Visit our website here: https://www.datadriveninvestor.com

Join our network here: https://datadriveninvestor.com/collaborate

Time Series Forecasting with Statistical Models was originally published in DataDrivenInvestor on Medium, where people are continuing the conversation by highlighting and responding to this story.

Many disciplines, such as finance, public administration, energy, retail, and healthcare, are dominated by time-series data. Large areas of micro- and macroeconomics rely on applied statistics with an emphasis on time-series analyses and modeling. The following are examples of time-series data:

• Daily closing values of a stock index
• Number of weekly infections of a disease
• Weekly series of train accidents
• Rainfall per day
• Sensor data such as temperature measurements per hour
• Population growth per year
• Quarterly earnings of a company over a number of years

Import Libraries

import itertools    # for grid search
import matplotlib.pyplot as plt  # for plotting graphs
import numpy as np         # for mathematical calculations
import pandas as pd     # for data manipulation
import statsmodels.api as sm        # for statistical models
from sklearn.metrics import mean_absolute_error # for evaluating models
from statsmodels.tsa.holtwinters import ExponentialSmoothing        # for Holt-Winters Exponential Smoothing
from statsmodels.tsa.holtwinters import SimpleExpSmoothing      # for Simple Exponential Smoothing
from statsmodels.tsa.seasonal import seasonal_decompose     # for seasonal decomposition
import statsmodels.tsa.api as smt       # for statistical tests
import warnings     # for ignoring warnings
warnings.filterwarnings('ignore')

data = sm.datasets.co2.load_pandas()        # 1958-03-01 - 2001-12-01
y = data.data                         # 1958-03-01 - 2001-12-01
y = y['co2'].resample('MS').mean()  # Monthly Mean

This code loads the CO2 data from the statsmodels library, and then selects the "co2" column of the data, which represents the CO2 concentration measurements. The resample method is then used to resample the data to a monthly frequency, taking the mean value for each month. This results in a new time series with the average CO2 concentration for each month.

y.isnull().sum()  #5
y = y.fillna(y.bfill())     # Backward Fill

y.describe()

count    526.000000
mean     339.624826
std       17.110954
min      313.400000
25%      324.025000
50%      337.912500
75%      354.537500
max      373.800000
Name: co2, dtype: float64

y.plot(figsize=(15, 6))
plt.show()

Holdout

“Holdout” is a commonly used method in machine learning for evaluating the performance of a model. It involves dividing the data into two parts: a training set and a testing set. The model is trained on the training set, and then the performance is evaluated on the testing set. The holdout method is used to estimate the performance of the model on unseen data, as the testing set is held out from the training process and used only for evaluation.

train = y[:'1997-12-01']        # 1958-03-01 - 1997-12-01
len(train)    #478

test = y['1998-01-01':]     # 1998-01-01 - 2001-12-01
len(test)    #48

Time Series Structural Analysis

Dickey-Fuller Test

def is_stationary(y):

    # "HO: Non-stationary"
    # "H1: Stationary"

    p_value = sm.tsa.stattools.adfuller(y)[1]
    if p_value < 0.05:
        print(F"Result: Stationary (H0: non-stationary, p-value: {round(p_value, 3)})")
    else:
        print(F"Result: Non-Stationary (H0: non-stationary, p-value: {round(p_value, 3)})")

is_stationary(y)

This code defines a function is_stationary which takes in a time series y as input. The function performs the Augmented Dickey-Fuller test to test the null hypothesis that the time series is non-stationary.

The p-value of the test is obtained using the adfuller function from the stattools module of statsmodels. If the p-value is less than 0.05, it suggests that the null hypothesis can be rejected, and the time series is stationary. In this case, the function prints "Result: Stationary". If the p-value is greater than 0.05, it suggests that the null hypothesis cannot be rejected, and the time series is non-stationary. In this case, the function prints "Result: Non-Stationary".

Result: Non-Stationary (H0: non-stationary, p-value: 0.999)

Time Series Components and Stationarity Test

def ts_decompose(y, model="additive", stationary=False):
    result = seasonal_decompose(y, model=model)
    fig, axes = plt.subplots(4, 1, sharex=True, sharey=False)
    fig.set_figheight(10)
    fig.set_figwidth(15)

    axes[0].set_title("Decomposition for " + model + " model")
    axes[0].plot(y, 'k', label='Original ' + model)
    axes[0].legend(loc='upper left')

    axes[1].plot(result.trend, label='Trend')
    axes[1].legend(loc='upper left')

    axes[2].plot(result.seasonal, 'g', label='Seasonality & Mean: ' + str(round(result.seasonal.mean(), 4)))
    axes[2].legend(loc='upper left')

    axes[3].plot(result.resid, 'r', label='Residuals & Mean: ' + str(round(result.resid.mean(), 4)))
    axes[3].legend(loc='upper left')
    plt.show(block=True)

    if stationary:
        is_stationary(y)

ts_decompose(y, stationary=True)

This code defines a function ts_decompose that takes in a time series y and optional parameters model and stationary. The function performs a time series decomposition on y using the seasonal_decompose function from the statsmodels library, with the specified decomposition model (either "additive" or "multiplicative").

The function then creates a plot with 4 subplots to visualize the original time series, the trend, the seasonal component, and the residuals. The plot title and labels are set to reflect the specified decomposition model.

If the stationary parameter is True, the function calls the is_stationary function to test the stationarity of the original time series.

Result: Non-Stationary (H0: non-stationary, p-value: 0.999)

1. Single Exponential Smoothing

ses_model = SimpleExpSmoothing(train).fit(smoothing_level=0.5)

y_pred = ses_model.forecast(48)

mean_absolute_error(test, y_pred)

This code fits a simple exponential smoothing (SES) model to the time series data in train. The SimpleExpSmoothing class from the statsmodels library is used to fit the model. The fit method is used to fit the model to the data, with a smoothing level of 0.5.

The model is then used to make 48-step ahead forecasts for the time series data in test. The forecasts are stored in the y_pred variable.

Finally, the function mean_absolute_error from the statsmodels library is used to calculate the mean absolute error (MAE) between the actual test data and the predicted values y_pred. The MAE is a measure of the average magnitude of the errors between the actual and predicted values. Lower values of MAE indicate better fit of the model to the data.

5.706393369643809

train.plot(title="Single Exponential Smoothing")
test.plot()
y_pred.plot()
plt.show()

train["1985":].plot(title="Single Exponential Smoothing")
test.plot()
y_pred.plot()
plt.show()

def plot_co2(train, test, y_pred, title):
    mae = mean_absolute_error(test, y_pred)
    train["1985":].plot(legend=True, label="TRAIN", title=f"{title}, MAE: {round(mae,2)}")
    test.plot(legend=True, label="TEST", figsize=(6, 4))
    y_pred.plot(legend=True, label="PREDICTION")
    plt.show()

plot_co2(train, test, y_pred, "Single Exponential Smoothing")

This code defines a function plot_co2 which plots the original time series data in train and test, and the predicted values in y_pred. The function also calculates the mean absolute error (MAE) between the actual test values and the predicted values y_pred, using the mean_absolute_error function from the statsmodels library. The MAE is then displayed in the plot title along with the title argument.

The function uses the plot method from the pandas library to plot the time series data. The legend parameter is set to True to display a legend in the plot, and the label parameter is used to specify the labels for each time series data. The figsize parameter is used to specify the size of the plot.

The function can be used to plot the results of any time series forecasting model, by passing in the relevant train, test, and y_pred arguments.

ses_model.params

{'smoothing_level': 0.5,
 'smoothing_trend': nan,
 'smoothing_seasonal': nan,
 'damping_trend': nan,
 'initial_level': 316.4419309772571,
 'initial_trend': nan,
 'initial_seasons': array([], dtype=float64),
 'use_boxcox': False,
 'lamda': None,
 'remove_bias': False}

2. Double Exponential Smoothing (DES)

Double Exponential Smoothing, also known as Holt’s method, is a time series forecasting technique that extends simple exponential smoothing by adding a trend component to the forecast model. It is used to handle time series data with a trend, such as sales data with an upward or downward trend over time. The method models the data as a combination of a level and a trend component. The level component captures the average value of the data, and the trend component captures changes in the level over time.

The method uses a weighted average of the past values and the trend to make predictions for future time steps. The weighting factors are determined by two smoothing parameters, alpha and beta, which control the weight given to the past values and the trend, respectively. The method can be implemented in a variety of ways, including simple exponential smoothing with a trend, and Holt-Winters with a seasonal period of 1.

ts_decompose(y)

des_model = ExponentialSmoothing(train, trend="add").fit(smoothing_level=0.5,
                                                         smoothing_trend=0.5)

y_pred = des_model.forecast(48)

plot_co2(train, test, y_pred, "Double Exponential Smoothing")

des_model.params

{'smoothing_level': 0.5,
 'smoothing_trend': 0.5,
 'smoothing_seasonal': nan,
 'damping_trend': nan,
 'initial_level': 317.5299919184764,
 'initial_trend': -0.3037511835408627,
 'initial_seasons': array([], dtype=float64),
 'use_boxcox': False,
 'lamda': None,
 'remove_bias': False}

3. Triple Exponential Smoothing (Holt-Winters)

Triple Exponential Smoothing, also known as the Holt-Winters method, is a time series forecasting technique that extends exponential smoothing by adding a seasonal component to the forecast model. It is used to handle time series data with a trend and seasonality, such as monthly sales data with a consistent upward trend and a repeating seasonal pattern from year to year.

The method models the data as a combination of a level, trend, and seasonality component. The level component captures the average value of the data, the trend component captures changes in the level over time, and the seasonality component captures repeating patterns in the data. The method uses a weighted average of the past values, level, and trend, to make predictions for future time steps.

tes_model = ExponentialSmoothing(train,
                                 trend="add",
                                 seasonal="add",
                                 seasonal_periods=12).fit(smoothing_level=0.5,
                                                          smoothing_slope=0.5,
                                                          smoothing_seasonal=0.5)

y_pred = tes_model.forecast(48)
plot_co2(train, test, y_pred, "Triple Exponential Smoothing")

Model Comparison

model_names = ["SES", "DES", "TES"]
model_mse = [round(np.mean((ses_model.resid)**2), 4),             
             round(np.mean((des_model.resid)**2), 4),             
             round(np.mean((tes_model.resid)**2), 4)]

plt.figure(figsize=(7, 5))
plt.bar(model_names, model_mse)
plt.xlabel('Models')
plt.ylabel('Mean Squared Error (MSE)')
plt.title('MSE of Different Time Series Models')
plt.show()

Conclusion

Time series forecasting is a powerful tool for organizations in a variety of industries. By using historical data to make predictions about future trends, organizations can make informed decisions and allocate resources effectively. There are many different types of time series forecasting models, ranging from simple methods to more complex models, and the choice of method will depend on the characteristics of the data and the goals of the problem. With the right tools and techniques, time series forecasting can provide valuable insights into future trends and help organizations plan for the future.

There are many resources available for learning about time series forecasting and how to perform it using Python. Here are a few recommended references:

• Time Series Forecasting in Python-by Marco Peixeiro
• Time Series Analysis with Python Cookbook: Practical Recipes for Exploratory Data Analysis, Data Preparation, Forecasting, and Model Evaluation-by Tarek A. Atwan
• Python for Time Series Data Analysis by Travis Mahoney
• Time Series Analysis and Forecasting with Python by Dr. PKS Prakash
• Hands-On Time Series Analysis with R by Vishwanathan N
• Time Series Analysis and Its Applications: With R Examples by Robert H. Shumway, David S. Stoffer

You can find the detailed part of the project on my Kaggle notebook.

References

1. https://www.miuul.com/zaman-serileri

🔔 If you’d like to get in touch, you can find me on Linkedin, drop me an email or follow me 👇

• Follow me: LinkedIn | Kaggle | Github

Subscribe to DDIntel Here.

Visit our website here: https://www.datadriveninvestor.com

Join our network here: https://datadriveninvestor.com/collaborate

Exponential Smoothing Techniques for Time Series Forecasting in Python: A Guide was originally published in DataDrivenInvestor on Medium, where people are continuing the conversation by highlighting and responding to this story.

Introduction

Unbalanced datasets refer to instances where the number of examples in one class is significantly different from the number of examples in another class. For example, in a binary classification problem where one class represents a small percentage of the total examples, this is known as an imbalanced dataset. In such scenarios, the class imbalance can lead to several issues while building a machine learning model.

One of the major issues with unbalanced datasets is that the model may be biased towards the majority class, leading to poor performance in predicting the minority class. This is because the model is trained to minimize the error rate, and when the majority class is over-represented, the model tends to predict the majority class more often. This leads to a higher accuracy score, but a poor recall and precision score for the minority class.

Another issue is that the model may not generalize well when exposed to new, unseen data. This is because the model is trained on a skewed dataset and may not be able to handle the imbalance in the test data.

In this article, we will discuss various tips and strategies to handle unbalanced datasets and improve the performance of machine learning models. Some of the techniques that will be covered include resampling techniques, cost-sensitive learning, using appropriate performance metrics, ensemble methods, and other strategies. By following these guidelines, we can build effective models for unbalanced datasets.

Tips for handling unbalanced datasets

Resampling techniques are one of the most popular ways to handle unbalanced datasets. These techniques involve either reducing the number of examples in the majority class or increasing the number of examples in the minority class.

Undersampling is a technique where we randomly remove examples from the majority class to reduce its size and balance the dataset. This technique is simple and easy to implement, but it can lead to information loss as it discards some of the majority class examples.

Oversampling is the opposite of undersampling, where we randomly replicate examples from the minority class to increase its size. This technique can lead to overfitting as the model is trained on repeated examples of the minority class.

SMOTE (Synthetic Minority Over-sampling Technique) is a more advanced technique that creates synthetic examples of the minority class rather than replicating existing examples. This technique helps to balance the dataset without introducing duplicates.

Cost-sensitive learning is another technique that can be used to handle unbalanced datasets. In this approach, different misclassification costs are assigned to different classes. This means that the model is penalized more heavily for misclassifying examples from the minority class than for misclassifying examples from the majority class.

Using appropriate performance metrics is also important when working with unbalanced datasets. Accuracy is not always the best metric as it can be misleading when dealing with unbalanced datasets. Instead, using a metric such as AUC-ROC (Area Under the Receiver Operating Characteristic Curve) can provide a better indication of model performance.

Ensemble methods, such as bagging and boosting, can also be effective for modeling unbalanced datasets. These methods combine the predictions of multiple models to improve the overall performance. Bagging involves training multiple models independently and averaging their predictions, while boosting involves training multiple models sequentially, where each model attempts to correct the errors of the previous model.

In summary, resampling techniques, cost-sensitive learning, using appropriate performance metrics, and ensemble methods are some of the tips and strategies that can help to handle unbalanced datasets and improve the performance of machine learning models.

Strategies for improving model performance on unbalanced datasets

Collecting more data is one of the most straightforward strategies for improving model performance on unbalanced datasets. By increasing the number of examples in the minority class, the model will have more information to learn from and will be less likely to be biased towards the majority class. This strategy is especially useful when the number of examples in the minority class is very small.

Generating synthetic samples is another strategy that can be used to improve model performance. Synthetic samples are artificially created examples that are similar to the real examples in the minority class. These samples can be generated using techniques such as SMOTE (Synthetic Minority Over-sampling Technique), which creates synthetic examples by interpolating between existing examples. Generating synthetic samples can help to balance the dataset and provide more examples for the model to learn from.

Using domain knowledge to focus on important samples is a strategy that can be used to improve model performance by identifying the most informative examples in the dataset. For example, if we are working on a medical dataset, we may know that certain symptoms or lab results are more indicative of a certain disease. By focusing on these examples, we can improve the model’s ability to accurately predict the minority class.

Finally, advanced techniques such as anomaly detection can be used to identify and focus on the minority class examples. These techniques can be used to identify examples that are different from the majority class and are likely to be minority class examples. This can help to improve model performance by identifying the most informative examples in the dataset.

In summary, collecting more data, generating synthetic samples, using domain knowledge to focus on important samples and using advanced techniques such as anomaly detection are some of the strategies that can be used to improve model performance on unbalanced datasets. These strategies can help to balance the dataset, provide more examples for the model to learn from, and identify the most informative examples in the dataset.

Let’s practise in python.

Import Libraries

import pandas as pd
import numpy as np
from sklearn.preprocessing import RobustScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix, classification_report,f1_score,recall_score,roc_auc_score, roc_curve
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib import rc,rcParams
import itertools

import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning) 
warnings.filterwarnings("ignore", category=FutureWarning) 
warnings.filterwarnings("ignore", category=UserWarning)

Read data

# Reading the dataset
df = pd.read_csv("creditcard.csv")
df.head()

# Number of variables and observations in the data set
print("Number of observations : " ,len(df))
print("Number of variables : ", len(df.columns))

Number of observations :  284807
Number of variables :  31

# We want to observe the types of variables in the dataset and whether they contain nulls
df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 284807 entries, 0 to 284806
Data columns (total 31 columns):
 #   Column  Non-Null Count   Dtype  
---  ------  --------------   -----  
 0   Time    284807 non-null  float64
 1   V1      284807 non-null  float64
 2   V2      284807 non-null  float64
 3   V3      284807 non-null  float64
 4   V4      284807 non-null  float64
 5   V5      284807 non-null  float64
 6   V6      284807 non-null  float64
 7   V7      284807 non-null  float64
 8   V8      284807 non-null  float64
 9   V9      284807 non-null  float64
 10  V10     284807 non-null  float64
 11  V11     284807 non-null  float64
 12  V12     284807 non-null  float64
 13  V13     284807 non-null  float64
 14  V14     284807 non-null  float64
 15  V15     284807 non-null  float64
 16  V16     284807 non-null  float64
 17  V17     284807 non-null  float64
 18  V18     284807 non-null  float64
 19  V19     284807 non-null  float64
 20  V20     284807 non-null  float64
 21  V21     284807 non-null  float64
 22  V22     284807 non-null  float64
 23  V23     284807 non-null  float64
 24  V24     284807 non-null  float64
 25  V25     284807 non-null  float64
 26  V26     284807 non-null  float64
 27  V27     284807 non-null  float64
 28  V28     284807 non-null  float64
 29  Amount  284807 non-null  float64
 30  Class   284807 non-null  int64  
dtypes: float64(30), int64(1)
memory usage: 67.4 MB

# The rate of presence of class 1 in the data set is 0.2%, and class 0 is 99.8%.

f,ax=plt.subplots(1,2,figsize=(18,8))
df['Class'].value_counts().plot.pie(explode=[0,0.1],autopct='%1.1f%%',ax=ax[0],shadow=True)
ax[0].set_title('dağılım')
ax[0].set_ylabel('')
sns.countplot('Class',data=df,ax=ax[1])
ax[1].set_title('Class')
plt.show()

# Standardizing Time and Amount variables
rob_scaler = RobustScaler()
df['Amount'] = rob_scaler.fit_transform(df['Amount'].values.reshape(-1,1))
df['Time'] = rob_scaler.fit_transform(df['Time'].values.reshape(-1,1))
df.head()

X = df.drop("Class", axis=1)
y = df["Class"]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.20, random_state=123456)

model = LogisticRegression(random_state=123456)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)
print("Accuracy: %.3f"%(accuracy))

Accuracy is the ratio of correct predictions to all predictions in the system. The accuracy score of the model we created is 0.999. We can say that our model works perfectly, right?

Let’s take a look at the Confusion Matrix to examine its performance.

Confusion Matrix

The Confusion Matrix is a table used to describe the performance of a classification model’s true values on test data. It contains 4 different combinations of estimated and actual values.

Terminology and derivations from a confusion matrix

• condition positive (P): the number of real positive cases in the data
• condition negative (N): the number of real negative cases in the data
• true positive (TP): A test result that correctly indicates the presence of a condition or characteristic
• true negative (TN): A test result that correctly indicates the absence of a condition or characteristic
• false positive (FP): A test result which wrongly indicates that a particular condition or attribute is present
• false negative (FN): A test result which wrongly indicates that a particular condition or attribute is absent

def plot_confusion_matrix(cm, classes,
                          title='Confusion matrix',
                          cmap=plt.cm.Blues):

    plt.rcParams.update({'font.size': 19})
    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    plt.title(title,fontdict={'size':'16'})
    plt.colorbar()
    tick_marks = np.arange(len(classes))
    plt.xticks(tick_marks, classes, rotation=45,fontsize=12,color="blue")
    plt.yticks(tick_marks, classes,fontsize=12,color="blue")
    rc('font', weight='bold')
    fmt = '.1f'
    thresh = cm.max()
    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
        plt.text(j, i, format(cm[i, j], fmt),
                 horizontalalignment="center",
                 color="red")

    plt.ylabel('True label',fontdict={'size':'16'})
    plt.xlabel('Predicted label',fontdict={'size':'16'})
    plt.tight_layout()

plot_confusion_matrix(confusion_matrix(y_test, y_pred=y_pred), classes=['Non Fraud','Fraud'],
                      title='Confusion matrix')

• A total of 56875 predictions were made for the Non Fraud class, with 56870 (TP) correct and 5 (FP) incorrect.

• A total of 87 predictions were made for the Fraud class, with 31 (FN) incorrect and 56 (TN) correct.

The model tells us that it can predict the Fraud state with 0.99 accuracies. But when we examine the confusion Matrix, the rate of false predictions in the Fraud class is quite high. While it is good at predicting the majority class, it is not good at predicting the minority class. In other words, the model correctly predicts the non-fraud class with a rate of 0.99. The fact that the number of observations belonging to the Non-Fraud class is higher than the number of observations of the fraud class causes the model to be successful in estimating the Non-Fraud class. With this observation we have made, we can say that the accuracy score is not a good performance measure for classification models, especially if it contains imbalances like the data set we have. We’ve looked at the dataset, we can look at how we can deal with the imbalance and build a model, what methods can be applied, and with what metrics we can measure its performance.

Choosing the Right Metric

When dealing with imbalanced datasets, it is important to choose the right metric to evaluate the performance of the model. Traditional metrics, such as accuracy, precision, and recall, may not be appropriate for imbalanced datasets because they do not take into account the distribution of the classes in the data.

One metric that is often used for imbalanced datasets is the F1 score. The F1 score is the harmonic mean of precision and recall, and it provides a balance between the two metrics. It is calculated as follows:

F1 = 2 * (precision * recall) / (precision + recall)

Another metric that is often used for imbalanced datasets is the area under the receiver operating characteristic curve (AUC-ROC). The AUC-ROC is a measure of the model’s ability to distinguish between positive and negative classes. It is calculated by plotting the true positive rate (TPR) against the false positive rate (FPR) at different classification thresholds. AUC-ROC values range from 0.5 (random guessing) to 1.0 (perfect classification).

print(classification_report(y_test, y_pred))

                precision   recall   f1-score   support

           0       1.00      1.00      1.00     56875
           1       0.92      0.64      0.76        87

    accuracy                           1.00     56962
   macro avg       0.96      0.82      0.88     56962
weighted avg       1.00      1.00      1.00     56962

Let’s examine the Precision measure for each class.

Returns how many of the predictions made for the 0 (non-fraud) class are correct. Looking at the confusion matrix, 56870 + 31 = 56901 non-fraud class predictions were made and 56870 of them were estimated correctly. The precision value for class 0 is 1 (56870 / 56901)

Returns how many of the predictions for the 1 (fraud) class are correct. Looking at the confusion matrix, 5 + 56 = 61 fraud class predictions were made and 56 of them were correctly estimated. The precision for class 0 is 0.92 (56 / 61)

Let’s examine the Recall measure for each class.

It shows how many of the values we need to predict for the 0(non-fraud) class are correctly predicted. We have 56870 + 5 = 56875 observations belonging to the non-fraud class and 56870 of them have been estimated correctly. The recall value for class 0 is 56870 / 56875 = 1.

It shows how many of the values we need to predict for the 1(fraud) class correctly. 31 + 56 = We have observations of 87 fraud classes and 56 of them are correctly estimated. The recall value for class 1 is 56 / 87 = 0.64.

When we look at the Recall values, we can easily see the failure of estimating 1 class. F1-score also expresses the harmonic mean of the recall and precision values.

Support refers to the number of actual values of the classes. We can show the structural weaknesses of the measurements, that is, we can say that the imbalance in the number of observations between the classes affects the measurements.

ROC curve

A receiver operating characteristic (ROC) curve is a graphical plot that shows the performance of a binary classification model at all classification thresholds. The ROC curve is created by plotting the true positive rate (TPR) against the false positive rate (FPR) at different classification thresholds. The true positive rate is the proportion of positive samples that are correctly classified as positive. It is calculated as the number of true positives divided by the number of actual positives. The false positive rate is the proportion of negative samples that are incorrectly classified as positive. It is calculated as the number of false positives divided by the number of actual negatives.

AUC (Area under the ROC curve)

AUC stands for “Area under the ROC Curve.” That is, AUC measures the entire two-dimensional area underneath the entire ROC curve (think integral calculus) from (0,0) to (1,1).

# Auc Roc Curve
def generate_auc_roc_curve(clf, X_test):
    y_pred_proba = clf.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test,  y_pred_proba)
    auc = roc_auc_score(y_test, y_pred_proba)
    plt.plot(fpr,tpr)
    plt.show()
    pass

generate_auc_roc_curve(model, X_test)

y_pred_proba = model.predict_proba(X_test)[:, 1]
auc = roc_auc_score(y_test, y_pred_proba)
print("AUC ROC Curve with Area Under the curve = %.3f"%auc)

AUC ROC Curve with Area Under the curve = 0.985

Resampling

Oversampling

Stabilizes the dataset by duplicating minority-class samples.

Random Oversampling:

It is the balancing of the data set by adding randomly selected samples from the minority class. This technique can be used if your dataset is small. It may cause overfitting. The RandomOverSampler method takes the sampling_strategy argument, and when sampling_stratefy=’minority’ is called, it increments the number of the minority class to equal the number of the majority class.

We can enter a float value in this argument. For example, let the number of our minority class be 1000 and the number of the majority class be 100. If we say sampling_stratefy = 0.5, the minority class will be added to be 500.

# Number of classes in the training set before random oversampling
y_train.value_counts()

0    227440
1       405
Name: Class, dtype: int64

 # Implementing RandomOver Sampling (Applying to the training set)
from imblearn.over_sampling import RandomOverSampler
oversample = RandomOverSampler(sampling_strategy='minority')
X_randomover, y_randomover = oversample.fit_resample(X_train, y_train)

# Number of classes of the training set after random oversampling
y_randomover.value_counts()

#training the model and its success rate
model.fit(X_randomover, y_randomover)
y_pred = model.predict(X_test)

accuracy = accuracy_score(y_test, y_pred)
print("Accuracy: %.3f%%" % (accuracy))

plot_confusion_matrix(confusion_matrix(y_test, y_pred=y_pred), classes=['Non Fraud','Fraud'],
                      title='Confusion matrix')

#classification report
print(classification_report(y_test, y_pred))

                 precision    recall  f1-score   support

           0       1.00      0.98      0.99     56875
           1       0.05      0.86      0.10        87

    accuracy                           0.98     56962
   macro avg       0.53      0.92      0.55     56962
weighted avg       1.00      0.98      0.99     56962

After applying Random Oversampling, the accuracy value of the trained model is 0.97, a decrease is observed. Looking at the Confusion Matrix and Classification report, the predicted fraud classes seem to have a high false rate, which lowered the precision of 1 class. But there is also an increase in the recall value of the 1 class, the rate of correctly estimating the fraud class of the model has increased. According to the first model, the prediction success of the Non fraud class has decreased, but the increase in the correct estimation of the fraud class is a big factor in our preference for the model created after randomoversampling.

SMOTE Oversampling:

Generating minority class synthetic samples to prevent overfitting.

First, a random sample from the minority class is selected. Then, k nearest neighbors are found for this sample. One of the k nearest neighbors is randomly selected and the synthetic sample is formed by combining it with the randomly selected sample from the minority class and forming a line segment in the feature space.

y_train.value_counts()

0    227440
1       405
Name: Class, dtype: int64

# Smote
from imblearn.over_sampling import SMOTE
oversample = SMOTE()
X_smote, y_smote = oversample.fit_resample(X_train, y_train)

# class number of training set after smote
y_smote.value_counts()

1    227440
0    227440
Name: Class, dtype: int64

model.fit(X_smote, y_smote)
y_pred = model.predict(X_test)

accuracy = accuracy_score(y_test, y_pred)
print("Accuracy: %.3f%%" % (accuracy))

Accuracy: 0.975%

plot_confusion_matrix(confusion_matrix(y_test, y_pred=y_pred), classes=['Non Fraud','Fraud'],
                      title='Confusion matrix')

#classification report
print(classification_report(y_test, y_pred))

precision    recall  f1-score   support

           0       1.00      0.98      0.99     56875
           1       0.05      0.86      0.10        87

    accuracy                           0.98     56962
   macro avg       0.53      0.92      0.54     56962
weighted avg       1.00      0.98      0.99     56962

Undersampling

It is a technique of balancing the data set by removing samples belonging to the majority class. It can be used when you have a large data set. Since the data set we have is not large, efficient results will not be obtained. But I will explain the methods and show how some of them can be applied.

Random Undersampling:

• Extracted samples are randomly selected.
• You can use this technique if you have a large data set.
• Information may be lost due to random selection.

# random undersampling
y_train.value_counts()

0    227440
1       405
Name: Class, dtype: int64

from imblearn.under_sampling import RandomUnderSampler
# transform the dataset
ranUnSample = RandomUnderSampler()
X_ranUnSample, y_ranUnSample = ranUnSample.fit_resample(X_train, y_train)

# After Random undersampling 
y_ranUnSample.value_counts()

Change the algorithm

from sklearn.ensemble import RandomForestClassifier
# train model
rfc = RandomForestClassifier(n_estimators=10).fit(X_train, y_train)
# predict on test set
rfc_pred = rfc.predict(X_test)
accuracy_score(y_test, rfc_pred)

0.9996313331694814

# f1 score
f1_score(y_test, rfc_pred)

0.8645161290322582

# confusion matrix
pd.DataFrame(confusion_matrix(y_test, rfc_pred))

# recall score
recall_score(y_test, rfc_pred)

0.7976190476190477

Conclusion

In this article, we discussed various tips and strategies for handling unbalanced datasets and improving the performance of machine learning models. Unbalanced datasets can be a common problem in machine learning and can lead to poor performance in predicting the minority class.

We discussed various resampling techniques such as undersampling, oversampling, and SMOTE that can be used to balance the dataset. We also discussed cost-sensitive learning and using appropriate performance metrics such as AUC-ROC, which can provide a better indication of model performance. Additionally, we discussed ensemble methods such as bagging and boosting that can also be effective for modeling unbalanced datasets.

We also discussed strategies for improving model performance on unbalanced datasets, such as collecting more data, generating synthetic samples, using domain knowledge to focus on important samples, and using advanced techniques such as anomaly detection.

In summary, unbalanced datasets can be challenging to work with, but by following the tips and strategies discussed in this article, we can build effective models that can accurately predict the minority class. It’s important to remember that the best approach will depend on the specific dataset and problem, and that a combination of techniques may be necessary to achieve the best results. Therefore, it’s important to experiment with different techniques and evaluate their performance using appropriate metrics.

References

1. https://imbalanced-learn.org/stable/introduction.html
2. https://courses.miuul.com/p/machine-learning
3. https://imbalanced-learn.org/dev/references/over_sampling.html

🔔 If you’d like to get in touch, you can find me on Linkedin, drop me an email or follow me 👇

• Follow me: LinkedIn | Kaggle | Github

Subscribe to DDIntel Here.

Visit our website here: https://www.datadriveninvestor.com

Join our network here: https://datadriveninvestor.com/collaborate

Modeling for Unbalanced Datasets: Tips and Strategies was originally published in DataDrivenInvestor on Medium, where people are continuing the conversation by highlighting and responding to this story.

Many disciplines, such as finance, public administration, energy, retail, and healthcare, are dominated by time-series data. Large areas of micro- and macroeconomics rely on applied statistics with an emphasis on time-series analyses and modeling.

Introduction

Time series data refers to a set of observations collected over a period of time. It is widely used in various fields such as finance, economics, meteorology, and more.

The following are examples of time-series data:
• Daily closing values of a stock index
• Number of weekly infections of a disease
• Weekly series of train accidents
• Rainfall per day
• Sensor data such as temperature measurements per hour
• Population growth per year
• Quarterly earnings of a company over a number of years

Time series data can be either continuous or discrete, and it can be univariate (one variable) or multivariate (multiple variables)

Characteristics:

• Time-dependent: the observations are collected at specific time intervals and the order of observations matters.
• Stationarity: the statistical properties of a time series do not change over time.
• Seasonality: a repeating pattern within a year or another fixed period of time.
• Trend: a long-term increase or decrease in the data.
• Noise: random fluctuations in the data.

Understanding these characteristics is important for selecting the appropriate time series modeling technique and for interpreting the results of the analysis.

Types of Time Series Data:

1. Additive time series: the individual components (trend, seasonality, and noise) can be separated and added together to obtain the total value.
2. Multiplicative time series: the individual components are multiplied together to obtain the total value.

Applications:

1. Financial data: stock prices, exchange rates, interest rates, etc.
2. Economic data: GDP, inflation, unemployment, etc.
3. Meteorological data: temperature, precipitation, wind speed, etc.
4. Healthcare data: hospital visits, prescription drug sales, etc.
5. Industrial data: production levels, energy consumption, etc.
6. Social media data: number of likes, shares, comments, etc.
7. Transportation data: traffic volume, fuel consumption, etc.

These are some examples of the types of data that can be analyzed using time series methods, and the specific type of data will affect the choice of model used in the analysis.

Time series modeling techniques including trend analysis, decomposition, and forecasting methods

Time series modeling techniques:

1. Trend analysis: identifies and models the long-term trend in the data. Techniques include moving averages, linear regression, and polynomial regression.
2. Decomposition: separates the time series into its individual components (trend, seasonality, and noise). Techniques include additive and multiplicative decomposition, and the classical decomposition of X-11.
3. Forecasting methods: use historical data to make predictions about future values. Techniques include simple exponential smoothing, Holt’s linear trend method, and ARIMA (Auto-Regressive Integrated Moving Average) model.

It’s important to note that many time series datasets are non-stationary and may require a technique like differencing to make it stationary before applying any of the above methods.

Importance of data preprocessing and cleaning in time series analysis

Data preprocessing and cleaning is an essential step in time series analysis as it helps to ensure the quality and reliability of the data. This is particularly important when working with real-world data, as it may contain missing values, outliers, and other errors.

Here are a few key aspects of data preprocessing and cleaning for time series analysis:

1. Missing values: Time series data often contains missing values, which can affect the accuracy of the analysis. These missing values must be handled appropriately, either by interpolation or by deletion.
2. Outliers: Outliers can have a significant impact on the results of time series analysis. They should be identified and handled by replacing them with a more appropriate value or by deletion.
3. Data transformation: Data may need to be transformed, such as taking the logarithm of the data, to stabilize the variance and make the data more suitable for analysis.
4. Data normalization: normalizing the data to a common scale is important for comparing different time series or for using some methods like ARIMA.
5. Data validation: it’s important to validate the data before analyzing it to ensure that it meets the assumptions of the chosen model.

By performing these steps, you can ensure that the data is of high quality, which will lead to more accurate and reliable results in the analysis.

Commonly used tools and software for time series analysis

There are several tools and software available for time series analysis, depending on the specific needs of the project and the level of complexity of the analysis. Some of the most commonly used tools and software include:

1. R: A programming language and software environment for statistical computing and graphics. It has a wide variety of packages for time series analysis, including “forecast”, “tseries”, “astsa”, “prophet” and more.
2. Python: A popular programming language that also has a wide variety of libraries for time series analysis such as “pandas”, “statsmodels”, “prophet”, “sklearn”, “numpy” and more.
3. Excel: a spreadsheet software that is widely used for simple time series analysis. It has built-in functions for handling time series data, such as moving averages, trendlines, and forecasting.
4. Tableau: a business intelligence and data visualization software, it can be used to create interactive dashboards and visualizations for time series data.
5. SAS: A widely used software for statistical analysis and data visualization, it has many capabilities for time series analysis like ARIMA, ETS, and more.
6. MATLAB: A multi-paradigm numerical computing environment and programming language, it has a wide range of toolboxes for time series analysis, signal processing, and control systems.

It’s also worth noting that some of these tools and software have online platforms such as RStudio, Jupyter, and Colab that can be used to perform time series analysis.

Real-world examples of time series analysis in various industries

Time series analysis is widely used in various industries to gain insights from historical data and make predictions about future events. Here are a few examples of how time series analysis is used in different industries:

1. Finance: Time series analysis is used to analyze stock prices, exchange rates, and interest rates. Financial institutions use it to identify trends and make predictions about future market conditions.
2. Economics: Time series analysis is used to analyze economic data such as GDP, inflation, and unemployment. This helps policymakers to understand the current state of the economy and make informed decisions.
3. Meteorology: Time series analysis is used to analyze weather data such as temperature, precipitation, and wind speed. This helps meteorologists to make weather predictions and issue alerts for severe weather conditions.
4. Healthcare: Time series analysis is used to analyze healthcare data such as hospital visits and prescription drug sales. This helps healthcare providers to understand patient needs and allocate resources more effectively.
5. Manufacturing: Time series analysis is used to analyze production levels, energy consumption, and other operational data. This helps manufacturers to optimize their production processes and reduce costs.
6. Retail: Retail companies use time series analysis on historical sales data to understand patterns in customer behavior, predict future sales and optimize inventory management.
7. Energy: Time series analysis is used to analyze energy consumption data, it helps energy companies to predict energy demand and optimize the operation of power plants, reduce costs and plan for the future.

These are just a few examples of how time series analysis is used in various industries, but the possibilities are endless. With the help of time series analysis, organizations can make data-driven decisions, optimize processes and improve their overall performance.

Conclusion

Time series analysis is a powerful tool for understanding historical data and making predictions about future events. By understanding the characteristics of time series data, selecting the appropriate modeling techniques, and properly preprocessing and cleaning the data, it is possible to gain valuable insights and make data-driven decisions.

🔔 If you’d like to get in touch, you can find me on Linkedin, drop me an email or follow me 👇

• Follow me: LinkedIn | Kaggle | Github

Data visualization is a crucial step in the data analysis process, as it allows us to quickly identify patterns and relationships within the data. However, creating visualizations can be a time-consuming and tedious task, especially when dealing with large and complex datasets. This is where Autoviz comes in. Autoviz is a powerful tool that allows users to easily visualize any dataset without the need for manual data manipulation or coding.

One of the key features of Autoviz is its ability to handle large and complex datasets. It can quickly identify patterns and relationships within the data, making it an ideal tool for data exploration and discovery. Additionally, Autoviz is highly customizable, allowing users to specify the types of visualizations they want to generate, as well as the color schemes and layouts of the charts.

The process of using Autoviz is simple and straightforward. Users can import their data into the tool using a variety of file formats, including CSV, Excel, and JSON. Once the data is imported, Autoviz will automatically generate a variety of visualizations, which can be further customized and refined by the user.

Another benefit of Autoviz is its user-friendliness. It is built with a graphical user interface that makes it easy for users of all skill levels to navigate and understand. Users can also create and share visualizations with others by exporting them to various file formats, such as PDF, PNG, and JPEG.

Import Libraries

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)

import warnings
warnings.filterwarnings('ignore')

The process of using Autoviz is simple and straightforward. First, you need to install Autoviz by running the following command in your terminal:

pip install autoviz

Once Autoviz is installed, you can import your data into the tool using a variety of file formats, including CSV, Excel, and JSON. For example, if your data is in a CSV file called “data.csv”, you can import it using the following code:

data = pd.read_csv('Iris.csv')
data.head()

Generate visualizations

Now that your data is imported, you can use Autoviz to generate visualizations. Autoviz have a function called “AutoViz_Class” that takes in a pandas DataFrame and generates visualizations. Here is an example of how to use it:

from autoviz.AutoViz_Class import AutoViz_Class
AV = AutoViz_Class()

%matplotlib inline # AutoViz no displays plots automatically. You must perform %matplotlib inline just before you run AutoViz on your data.

# Generate visualizations
AV.AutoViz("/kaggle/input/iris/Iris.csv")

This will generate a variety of visualizations, including bar charts, line graphs, and scatter plots. These visualizations will be saved in the current working directory as HTML files that can be opened in a web browser for viewing.

6 Predictors classified...
        1 variables removed since they were ID or low-information variables
        List of variables removed: ['Id']
Number of All Scatter Plots = 10

All Plots done
Time to run AutoViz = 6 seconds 

 ###################### AUTO VISUALIZATION Completed ########################

Note: Autoviz only works for continuous and categorical variables. You can use dep_var to define continuous variables in your dataset.

This will generate a variety of visualizations, including bar charts, line graphs, and scatter plots. These visualizations will be saved in the current working directory as HTML files that can be opened in a web browser for viewing. In addition to generating visualizations, Autoviz also allows users to customize and refine their visualizations. For example, you can specify the types of visualizations you want to generate, as well as the color schemes and layouts of the charts.

Parameters

• filename: The name of the file containing the dataset.
• sep: The delimiter used in the file.
• depVar: The name of the target variable.
• dfte: The dataframe containing the dataset.
• header: The row number of the header in the file.
• verbose: The level of verbosity. If it is 0 or 1, it only shows the plots on the notebook. If we give 2, it also saves the plots.
• lowess: A boolean indicating whether to use the LOWESS smoother for scatter plots.
• chart_format: The format of the generated charts. In this case, it is 'png'.
• max_rows_analyzed: The maximum number of rows to be analyzed from the dataset.
• max_cols_analyzed: The maximum number of columns to be analyzed from the dataset.

Setting a Target Variable

While visualizing the data, we can give a parameter called depVar. In this way, we can determine a target variable and create plots accordingly. You can also specify the columns of the dataframe to be plotted and excluded in the visualization.

Here is an example of how to do this:

AV = AutoViz_Class()
dftc = AV.AutoViz(
    filename='', 
    sep='' , 
    depVar='Species', 
    dfte=data, 
    header=0, 
    verbose=1, 
    lowess=False, 
    chart_format='png', 
    max_rows_analyzed=300000, 
    max_cols_analyzed=30
)

Shape of your Data Set loaded: (150, 6)
#######################################################################################
######################## C L A S S I F Y I N G  V A R I A B L E S  ####################
#######################################################################################
Classifying variables in data set...
Data cleaning improvement suggestions. Complete them before proceeding to ML modeling.

Multi_Classification problem

AutoViz provides data cleaning suggestions

AutoViz automatically analyzes your dataset and provides suggestions for how to clean your data set. It detects missing values, identifies rare categories, finds infinite values, detects mixed data types, and so much more. This will help you tremendously speed up your data cleaning activities.

from autoviz import data_cleaning_suggestions
data_cleaning_suggestions(data)

Besides these AutoViz can create Wordclouds automatically for your NLP variables in data. It detects NLP variables automatically and creates wordclouds for them. And AutoViz runs on Docker containers as part of MLOps pipelines. Check out Orchest.io

Conclusion

In conclusion, Autoviz is a powerful and user-friendly tool that makes data visualization easy and effortless. Its ability to handle large and complex datasets, highly customizable and ability to specify columns to be plotted makes it an ideal tool for data exploration and discovery. By using Autoviz, data scientists and analysts can save time and effort when creating visualizations, allowing them to focus on interpreting the data and uncovering insights.

Thanks for reading this article. You can access the detailed codes of the project and other projects on my Github account or Kaggle account. Happy coding!

Please feel free to contact me if you need any further information.

References

1. https://pypi.org/project/autoviz/
2. https://github.com/AutoViML/AutoViz
3. https://github.com/rsesha/autoviz_pipeline

Subscribe to DDIntel Here.

Visit our website here: https://www.datadriveninvestor.com

Join our network here: https://datadriveninvestor.com/collaborate

Autoviz: The Key to Effortless Data Visualization was originally published in DataDrivenInvestor on Medium, where people are continuing the conversation by highlighting and responding to this story.

Hello everyone!

In this post, we will explore the use of Python to analyze and predict World Happiness. The World Happiness Report, which is published by the United Nations Sustainable Development Solutions Network, ranks 156 countries by their overall happiness levels.

We will begin by importing the necessary libraries, including pandas for data manipulation and matplotlib for visualization. Next, we will load in the World Happiness Report dataset and take a look at the data to understand the different factors that contribute to a country’s overall happiness score.

Using statistical analysis techniques, such as correlation and regression analysis, we will then identify the strongest predictors of happiness. We can use machine learning algorithms like linear regression, Random Forest and Decision Tree.

Finally, we will use our findings to create a predictive model that can estimate a country’s happiness score based on its other characteristics.

This post will serve as a comprehensive guide for those looking to learn about data analysis and prediction using Python in the context of World Happiness. By the end of this post, readers will have a better understanding of the factors that contribute to a country’s overall happiness, and will be able to create their own predictive models using Python.

The World Happiness Report is a survey that ranks countries by their happiness levels, using data from the Gallup World Poll. The report uses the Cantril ladder question, which asks respondents to rate their current lives on a scale of 0 to 10, with 10 being the best possible life. The report also includes estimates of how six factors — economic production, social support, life expectancy, freedom, absence of corruption, and generosity — contribute to making life evaluations higher in each country. The report also includes a hypothetical country called Dystopia, which has the lowest national averages for each of the six factors, and serves as a benchmark against which all countries can be compared.

Variables

• Overall rank: Rank of countries based on the Happiness Score
• Score: A metric measured in 2019 by asking the sampled people the question: “How would you rate your happiness on a scale of 0 to 10 where 10 is the happiest.”
• GDP per capita: The extent to which GDP contributes to the calculation of the Happiness Score.
• Social support: The extent to which Social support contributes to the calculation of the Happiness Score
• Healthy life expectancy: The extent to which Life expectancy contributed to the calculation of the Happiness Score
• Freedom to make life choices: The extent to which Freedom contributed to the calculation of the Happiness Score
• Generosity: The extent to which Generosity contributed to the calculation of the Happiness Score
• Perceptions of corruption: The extent to which Perception of Corruption contributes to Happiness Score.

Import libraries

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
import warnings
warnings.filterwarnings('ignore')

Read data

df = pd.read_csv('/kaggle/input/world-happiness/2019.csv')
df.head()

Summary Statistics

To get a summary of the numerical columns in the data, you can use the some function.

This code defines a function called check_df that takes a DataFrame as input and prints various information about the DataFrame, including its shape, data types, head, tail, missing values, and quantiles.

The function has two parameters: dataframe, which is the DataFrame to be analyzed, and head, which is the number of rows to print for the head and tail of the DataFrame. The default value for head is 5.

def check_df(dataframe, head=5):
    print("##################### Shape #####################")
    print(dataframe.shape)
    print("##################### Types #####################")
    print(dataframe.dtypes)
    print("##################### Head #####################")
    print(dataframe.head(head))
    print("##################### Tail #####################")
    print(dataframe.tail(head))
    print("##################### NA #####################")
    print(dataframe.isnull().sum())
    print("##################### Quantiles #####################")
    print(dataframe.quantile([0, 0.05, 0.50, 0.95, 0.99, 1]).T)

check_df(df)

##################### Shape #####################
(156, 9)
##################### Types #####################
Overall rank                      int64
Country or region                object
Score                           float64
GDP per capita                  float64
Social support                  float64
Healthy life expectancy         float64
Freedom to make life choices    float64
Generosity                      float64
Perceptions of corruption       float64
dtype: object
##################### Head #####################
   Overall rank Country or region  Score  GDP per capita  Social support  \
0             1           Finland  7.769           1.340           1.587   
1             2           Denmark  7.600           1.383           1.573   
2             3            Norway  7.554           1.488           1.582   
3             4           Iceland  7.494           1.380           1.624   
4             5       Netherlands  7.488           1.396           1.522   

   Healthy life expectancy  Freedom to make life choices  Generosity  \
0                    0.986                         0.596       0.153   
1                    0.996                         0.592       0.252   
2                    1.028                         0.603       0.271   
3                    1.026                         0.591       0.354   
4                    0.999                         0.557       0.322   

   Perceptions of corruption  
0                      0.393  
1                      0.410  
2                      0.341  
3                      0.118  
4                      0.298  
##################### Tail #####################
     Overall rank         Country or region  Score  GDP per capita  \
151           152                    Rwanda  3.334           0.359   
152           153                  Tanzania  3.231           0.476   
153           154               Afghanistan  3.203           0.350   
154           155  Central African Republic  3.083           0.026   
155           156               South Sudan  2.853           0.306   

     Social support  Healthy life expectancy  Freedom to make life choices  \
151           0.711                    0.614                         0.555   
152           0.885                    0.499                         0.417   
153           0.517                    0.361                         0.000   
154           0.000                    0.105                         0.225   
155           0.575                    0.295                         0.010   

     Generosity  Perceptions of corruption  
151       0.217                      0.411  
152       0.276                      0.147  
153       0.158                      0.025  
154       0.235                      0.035  
155       0.202                      0.091  
##################### NA #####################
Overall rank                    0
Country or region               0
Score                           0
GDP per capita                  0
Social support                  0
Healthy life expectancy         0
Freedom to make life choices    0
Generosity                      0
Perceptions of corruption       0
dtype: int64
##################### Quantiles #####################
                               0.00     0.05     0.50       0.95       0.99  \
Overall rank                  1.000  8.75000  78.5000  148.25000  154.45000   
Score                         2.853  3.48150   5.3795    7.28525    7.57470   
GDP per capita                0.000  0.25200   0.9600    1.44150    1.58865   
Social support                0.000  0.64325   1.2715    1.53800    1.58425   
Healthy life expectancy       0.000  0.28825   0.7890    1.03975    1.10330   
Freedom to make life choices  0.000  0.14275   0.4170    0.58425    0.60570   
Generosity                    0.000  0.04675   0.1775    0.34950    0.45455   
Perceptions of corruption     0.000  0.01900   0.0855    0.32225    0.41045   

                                 1.00  
Overall rank                  156.000  
Score                           7.769  
GDP per capita                  1.684  
Social support                  1.624  
Healthy life expectancy         1.141  
Freedom to make life choices    0.631  
Generosity                      0.566  
Perceptions of corruption       0.453

The describe method is a function of the Pandas DataFrame class that generates descriptive statistics of a DataFrame.

When you call the describe method on a DataFrame, it returns a new DataFrame that contains various statistical summary measures of the data in the original DataFrame. These summary measures include the count, mean, standard deviation, minimum, maximum, 25th percentile, 50th percentile (also known as the median), and 75th percentile of the data.

# Generate statistics
df.describe()

Remove outliers

This code defines a function called remove_outliers that takes a DataFrame and a column name as input and returns a new DataFrame with the outliers in the specified column removed.

The function works by first sorting the DataFrame by the specified column. It then calculates the first quartile (Q1) and third quartile (Q3) of the column, as well as the interquartile range (IQR) which is the difference between Q3 and Q1. The function then calculates the lower and upper range of the column data by adding and subtracting 1.5 * IQR from Q1 and Q3, respectively.

Finally, the function filters the original DataFrame to include only rows where the value in the specified column is within the lower and upper range, and returns the resulting DataFrame.

def remove_outliers(df, column):
    sorted(df)
    Q1, Q3 = np.percentile(df[column], [25, 75])
    IQR = Q3 - Q1
    lower_range = Q1 - (1.5 * IQR)
    upper_range = Q3 + (1.5 * IQR)
    df = df[(df[column] > lower_range) & (df[column] < upper_range)]
    return df

df = remove_outliers(df, 'GDP per capita')
df = remove_outliers(df, 'Social support')
df = remove_outliers(df, 'Healthy life expectancy')
df = remove_outliers(df, 'Freedom to make life choices')
df = remove_outliers(df, 'Generosity')
df = remove_outliers(df, 'Perceptions of corruption')

# Check if there is any outliers

df.boxplot(figsize=(20, 8))
plt.show()

# Correlation matrix

corr_matrix = df.corr()
sns.heatmap(corr_matrix, annot=True, cmap='coolwarm')
plt.show()

# Scatter plot

sns.scatterplot(x = 'GDP per capita', y = 'Happiness Score', data = df)
plt.show()

sns.scatterplot(x = 'Social support', y = 'Happiness Score', data = df)
plt.show()

sns.scatterplot(x = 'Healthy life expectancy', y = 'Happiness Score', data = df)
plt.show()

sns.scatterplot(x = 'Freedom to make life choices', y = 'Happiness Score', data = df)
plt.show()

sns.scatterplot(x = 'Generosity', y = 'Happiness Score', data = df)
plt.show()

sns.scatterplot(x = 'Perceptions of corruption', y = 'Happiness Score', data = df)
plt.show()

# Histogram

sns.distplot(df['Happiness Score'])
plt.show()

Barplot

This code uses the Matplotlib and Seaborn libraries to create a bar plot that visualizes the relationship between the Country and Happiness Score columns in the DataFrame df.

The figure function from Matplotlib is used to set the size of the plot, and the barplot function from Seaborn is used to create the bar plot. The x and y parameters specify the variables to use for the x-axis and y-axis, respectively, and the data parameter specifies the DataFrame to use for the plot.

The xticks function from Matplotlib is used to rotate the labels on the x-axis, and the show function is used to display the plot.

plt.figure(figsize=(25, 10))
sns.barplot(x='Country', y='Happiness Score', data=df)
plt.xticks(rotation=90)
plt.show()

# Top 10 countries with the highest happiness score

top10 = df.sort_values(by = 'Happiness Score', ascending = False).head(10)
plt.figure(figsize=(20, 10))
sns.barplot(x = 'Country', y = 'Happiness Score', data = top10)
plt.xticks(rotation = 90)
plt.show()

# Top 10 countries with the highest GDP per capita

top10_gdp = df.sort_values(by = 'GDP per capita', ascending = False).head(10)
plt.figure(figsize=(25, 10))
sns.barplot(x = 'Country', y = 'GDP per capita', data = top10_gdp)
plt.xticks(rotation = 90)
plt.show()

# Top 10 countries with the highest social support

top10_social = df.sort_values(by = 'Social support', ascending = False).head(10)
plt.figure(figsize=(25, 10))
sns.barplot(x = 'Country', y = 'Social support', data = top10_social)
plt.xticks(rotation = 90)
plt.show()

# Top 10 countries with the highest healthy life expectancy

top10_health = df.sort_values(by = 'Healthy life expectancy', ascending = False).head(10)
plt.figure(figsize=(25, 10))
sns.barplot(x = 'Country', y = 'Healthy life expectancy', data = top10_health)
plt.xticks(rotation = 90)
plt.show()

from sklearn.model_selection import train_test_split

# Splitting the data into features and target
X = df.drop(['Country', 'Happiness Score'], axis=1)
y = df['Happiness Score']

# Splitting the dataset into the Training set and Test set
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state = 0)

# Feature Scaling
sc = StandardScaler()
X_train = sc.fit_transform(X_train)
X_test = sc.transform(X_test)

This code uses the train_test_split function from the sklearn.model_selection module to split a dataset into a training set and a test set. The X variable contains the features of the dataset, and the y variable contains the target values. The train_test_split function takes the features and target values as input, along with the test size (in this case, 20% of the data is used for the test set) and a random seed (to ensure reproducibility).

The StandardScaler class from the sklearn.preprocessing module is then used to perform feature scaling on the training and test sets. Feature scaling is a technique used to standardize the range of independent variables or features of data. It is important to scale the features in a machine learning model to ensure that no single feature dominates the others and to avoid any bias in the model.

To use the StandardScaler, an instance of the class is first created and then fit to the training set using the fit_transform method. The resulting transformed data is stored in X_train. The transform method is then used to apply the same transformation to the test set, and the resulting transformed data is stored in X_test.

Model Building

This code imports several classes from the xgboost, sklearn.linear_model, sklearn.tree, sklearn.ensemble, and sklearn.ensemble modules, which are commonly used for building regression models in Python.

• The XGBRegressor class from the xgboost library is a gradient boosting model that uses decision trees as its base learners. It is a popular choice for building regression models because of its good performance and ability to handle large datasets efficiently.
• The LinearRegression class from the sklearn.linear_model module is a simple linear regression model that fits a linear equation to the data. It is a good choice for modeling linear relationships between variables.
• The DecisionTreeRegressor class from the sklearn.tree module is a decision tree model that splits the data into groups based on the value of a feature. It is a good choice for modeling non-linear relationships between variables.
• The RandomForestRegressor class from the sklearn.ensemble module is an ensemble model that combines multiple decision trees and averages their predictions. It is a good choice for improving the performance of a decision tree model and reducing overfitting.
• The GradientBoostingRegressor class from the sklearn.ensemble module is another gradient boosting model that uses decision trees as its base learners. It is similar to the XGBRegressor class, but has more options for model tuning and is slower to train.
• The AdaBoostRegressor class from the sklearn.ensemble module is an ensemble model that combines multiple weak models and weights their predictions to make a final prediction. It is a good choice for improving the performance of a simple model.

from xgboost import XGBRegressor
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.ensemble import AdaBoostRegressor
from sklearn.neighbors import KNeighborsRegressor

# Linear Regression
lr = LinearRegression()
lr.fit(X_train, y_train)
lr_pred = lr.predict(X_test)
lr_score = lr.score(X_test, y_test)
print('Linear Regression Score: ', lr_score)

# Decision Tree
dt = DecisionTreeRegressor()
dt.fit(X_train, y_train)
dt_pred = dt.predict(X_test)
dt_score = dt.score(X_test, y_test)
print('Decision Tree Score: ', dt_score)

# Random Forest
rf = RandomForestRegressor()
rf.fit(X_train, y_train)
rf_pred = rf.predict(X_test)
rf_score = rf.score(X_test, y_test)
print('Random Forest Score: ', rf_score)

# XGBoost
xgb = XGBRegressor()
xgb.fit(X_train, y_train)
xgb_pred = xgb.predict(X_test)
xgb_score = xgb.score(X_test, y_test)
print('XGBoost Score: ', xgb_score)

# RandomForestRegressor
rf = RandomForestRegressor()
rf.fit(X_train, y_train)
rf_pred = rf.predict(X_test)
rf_score = rf.score(X_test, y_test)
print('RandomForestRegressor Score: ', rf_score)

# GradientBoostingRegressor
gb = GradientBoostingRegressor()
gb.fit(X_train, y_train)
gb_pred = gb.predict(X_test)
gb_score = gb.score(X_test, y_test)
print('GradientBoostingRegressor Score: ', gb_score)

# AdaBoostRegressor
ada = AdaBoostRegressor()
ada.fit(X_train, y_train)
ada_pred = ada.predict(X_test)
ada_score = ada.score(X_test, y_test)
print('AdaBoostRegressor Score: ', ada_score)

# KNeighborsRegressor
knn = KNeighborsRegressor()
knn.fit(X_train, y_train)
knn_pred = knn.predict(X_test)
knn_score = knn.score(X_test, y_test)
print('KNeighborsRegressor Score: ', knn_score)

Linear Regression Score:  0.6841216592604057
Decision Tree Score:  0.28616030271278814
Random Forest Score:  0.69443250417991
XGBoost Score:  0.5981049248655365
RandomForestRegressor Score:  0.7092322318852047
GradientBoostingRegressor Score:  0.6715963944352763
AdaBoostRegressor Score:  0.6841751067060653
KNeighborsRegressor Score:  0.6816380800543443

Model Evaluation

Evaluating a model’s performance is an important step in any data analysis or machine learning project. There are many different ways to evaluate a model’s performance, depending on the specific problem you are trying to solve and the type of model you are using.

For example, if you are working on a classification problem (e.g. predicting whether someone is happy or not), you might use evaluation metrics such as accuracy, precision, recall, and F1 score to measure how well your model is performing.

On the other hand, if you are working on a regression problem (e.g. predicting happiness on a continuous scale), you might use evaluation metrics such as mean absolute error (MAE), mean squared error (MSE), and root mean squared error (RMSE) to measure the model’s performance.

To evaluate a model’s performance, you will typically split your data into a training set and a test set. You will train your model on the training set and then use the test set to evaluate the model’s performance. This helps ensure that the model is not overfitting to the training data and can generalize well to new data.

It is also important to consider the context of your analysis and the goals of your model. For example, if you are trying to predict happiness around the world, you might want to consider factors such as cultural differences and the impact of specific events on happiness levels. You might also want to consider using multiple evaluation metrics to get a more comprehensive view of the model’s performance.

from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score
from sklearn.metrics import mean_absolute_error

The code you provided will evaluate the performance of a linear regression model using four different evaluation metrics: mean squared error (MSE), root mean squared error (RMSE), mean absolute error (MAE), and R2 score.

Here is a brief explanation of each metric:

• Mean squared error (MSE) is a measure of the average squared difference between the predicted values and the true values. It is calculated by taking the sum of the squared differences and dividing by the number of samples. Lower values of MSE indicate a better fit.
• Root mean squared error (RMSE) is the square root of the MSE. It is a measure of the average magnitude of the error. Lower values of RMSE indicate a better fit.
• Mean absolute error (MAE) is a measure of the average magnitude of the error. It is calculated by taking the absolute value of the difference between the predicted values and the true values and averaging over all samples. Lower values of MAE indicate a better fit.
• R2 score (also known as the coefficient of determination) is a measure of the proportion of the variance in the dependent variable that is explained by the independent variables. It ranges from 0 to 1, with higher values indicating a better fit. An R2 score of 1 indicates that the model perfectly explains the variance in the dependent variable.

# Linear Regression
print('Linear Regression')
print('Mean Squared Error:', mean_squared_error(y_test, lr_pred))
print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, lr_pred)))
print('Mean Absolute Error:', mean_absolute_error(y_test, lr_pred))
print('R2 Score:', r2_score(y_test, lr_pred))

print('*************************')

# Decision Tree
print('Decision Tree')
print('Mean Squared Error:', mean_squared_error(y_test, dt_pred))
print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, dt_pred)))
print('Mean Absolute Error:', mean_absolute_error(y_test, dt_pred))
print('R2 Score:', r2_score(y_test, dt_pred))

print('*************************')

# Random Forest
print('Random Forest')
print('Mean Squared Error:', mean_squared_error(y_test, rf_pred))
print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, rf_pred)))
print('Mean Absolute Error:', mean_absolute_error(y_test, rf_pred))
print('R2 Score:', r2_score(y_test, rf_pred))

print('*************************')

# XGBoost
print('XGBoost')
print('Mean Squared Error:', mean_squared_error(y_test, xgb_pred))
print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, xgb_pred)))
print('Mean Absolute Error:', mean_absolute_error(y_test, xgb_pred))
print('R2 Score:', r2_score(y_test, xgb_pred))

print('*************************')  # Add a newline between the output for the two models

# RandomForestRegressor
print('RandomForestRegressor')
print('Mean Squared Error:', mean_squared_error(y_test, rf_pred))
print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, rf_pred)))
print('Mean Absolute Error:', mean_absolute_error(y_test, rf_pred))
print('R2 Score:', r2_score(y_test, rf_pred))

print('*************************')

# GradientBoostingRegressor
print('GradientBoostingRegressor')
print('Mean Squared Error:', mean_squared_error(y_test, gb_pred))
print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, gb_pred)))
print('Mean Absolute Error:', mean_absolute_error(y_test, gb_pred))
print('R2 Score:', r2_score(y_test, gb_pred))

print('*************************')

# AdaBoostRegressor
print('AdaBoostRegressor')
print('Mean Squared Error:', mean_squared_error(y_test, ada_pred))
print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, ada_pred)))
print('Mean Absolute Error:', mean_absolute_error(y_test, ada_pred))
print('R2 Score:', r2_score(y_test, ada_pred))

print('*************************')

# KNeighborsRegressor
print('KNeighborsRegressor')
print('Mean Squared Error:', mean_squared_error(y_test, knn_pred))
print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, knn_pred)))
print('Mean Absolute Error:', mean_absolute_error(y_test, knn_pred))
print('R2 Score:', r2_score(y_test, knn_pred))

Linear Regression
Mean Squared Error: 0.3029439735952702
Root Mean Squared Error: 0.5504034643743353
Mean Absolute Error: 0.42688322026284425
R2 Score: 0.6841216592604057
*************************
Decision Tree
Mean Squared Error: 0.5624305769230769
Root Mean Squared Error: 0.7499537165206109
Mean Absolute Error: 0.5437307692307694
R2 Score: 0.41355612619967286
*************************
Random Forest
Mean Squared Error: 0.2824164245730768
Root Mean Squared Error: 0.5314286636728177
Mean Absolute Error: 0.38225346153846146
R2 Score: 0.7055256437913671
*************************
XGBoost
Mean Squared Error: 0.38543855442742897
Root Mean Squared Error: 0.6208369789465097
Mean Absolute Error: 0.5258417670910175
R2 Score: 0.5981049248655365
*************************
RandomForestRegressor
Mean Squared Error: 0.2824164245730768
Root Mean Squared Error: 0.5314286636728177
Mean Absolute Error: 0.38225346153846146
R2 Score: 0.7055256437913671
*************************
GradientBoostingRegressor
Mean Squared Error: 0.31603251429639845
Root Mean Squared Error: 0.562167692327119
Mean Absolute Error: 0.44836291139259415
R2 Score: 0.67047429578818
*************************
AdaBoostRegressor
Mean Squared Error: 0.28672776810554257
Root Mean Squared Error: 0.535469670574854
Mean Absolute Error: 0.40385947538218353
R2 Score: 0.701030224967776
*************************
KNeighborsRegressor
Mean Squared Error: 0.30532585692307673
Root Mean Squared Error: 0.5525629890999547
Mean Absolute Error: 0.3932692307692307
R2 Score: 0.6816380800543443

Model Tuning

Model tuning, also known as hyperparameter optimization or hyperparameter tuning, refers to the process of finding the best values for the hyperparameters of a machine learning model. Hyperparameters are parameters that are set before the model is trained and are not learned from the data during training. They control the model’s behavior and can significantly impact its performance.

from sklearn.model_selection import GridSearchCV

# Random Forest
rf_params = {"max_depth": [3, 5, 8, None],
                "max_features": [3, 5, 10],
                "min_samples_split": [2, 5, 10],
                "min_samples_leaf": [1, 3, 5, 10],
                "bootstrap": [True, False],
                "n_estimators": [100, 500, 1000],
                "random_state": [42]}
rf_cv_model = GridSearchCV(rf, rf_params, cv=5, n_jobs=-1, verbose=2)
rf_cv_model.fit(X_train, y_train)
print("Best Parameters: " + str(rf_cv_model.best_params_))

# GradientBoostingRegressor
gb_params = {"learning_rate": [0.01, 0.1, 0.05],
                "max_depth": [3, 5, 8],
                "max_features": [3, 5, 10],
                "min_samples_split": [2, 5, 10],
                "min_samples_leaf": [1, 3, 5, 10],
                "n_estimators": [100, 500, 1000],
                "subsample": [1, 0.5, 0.75],
                "random_state": [42]}
gb_cv_model = GridSearchCV(gb, gb_params, cv=5, n_jobs=-1, verbose=2)
gb_cv_model.fit(X_train, y_train)
print("Best Parameters: " + str(gb_cv_model.best_params_))

# AdaBoostRegressor
ada_params = {"learning_rate": [0.01, 0.1, 0.05],
                "n_estimators": [100, 500, 1000],
                "random_state": [42]}
ada_cv_model = GridSearchCV(ada, ada_params, cv=5, n_jobs=-1, verbose=2)
ada_cv_model.fit(X_train, y_train)
print("Best Parameters: " + str(ada_cv_model.best_params_))

# KNeighborsRegressor
knn_params = {"n_neighbors": np.arange(1, 50)}
knn_cv_model = GridSearchCV(knn, knn_params, cv=5, n_jobs=-1, verbose=2)
knn_cv_model.fit(X_train, y_train)
print("Best Parameters: " + str(knn_cv_model.best_params_))

Fitting 5 folds for each of 9 candidates, totalling 45 fits
Best Parameters: {'learning_rate': 0.05, 'n_estimators': 100, 'random_state': 42}
Fitting 5 folds for each of 49 candidates, totalling 245 fits
Best Parameters: {'n_neighbors': 3}

Model Comparison

models = [
    RandomForestRegressor(),
    GradientBoostingRegressor(),
    AdaBoostRegressor(),
    KNeighborsRegressor()
]

for model in models:
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    print(model)
    print('Mean Squared Error:', mean_squared_error(y_test, y_pred))
    print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, y_pred)))
    print('Mean Absolute Error:', mean_absolute_error(y_test, y_pred))
    print('R2 Score:', r2_score(y_test, y_pred))
    print('---------------------------------------------------')

RandomForestRegressor()
Mean Squared Error: 0.2892094872846152
Root Mean Squared Error: 0.5377820072153914
Mean Absolute Error: 0.39066615384615405
R2 Score: 0.6984425473613733
---------------------------------------------------
GradientBoostingRegressor()
Mean Squared Error: 0.31686866643486644
Root Mean Squared Error: 0.5629108867617204
Mean Absolute Error: 0.44874855982340944
R2 Score: 0.6696024436533756
---------------------------------------------------
AdaBoostRegressor()
Mean Squared Error: 0.29817323290920117
Root Mean Squared Error: 0.5460524085737569
Mean Absolute Error: 0.3969715415454
R2 Score: 0.6890960894632248
---------------------------------------------------
KNeighborsRegressor()
Mean Squared Error: 0.30532585692307673
Root Mean Squared Error: 0.5525629890999547
Mean Absolute Error: 0.3932692307692307
R2 Score: 0.6816380800543443
---------------------------------------------------

# Random Forest
rf_tuned = RandomForestRegressor(bootstrap=False, max_depth=8, max_features=5, min_samples_leaf=1, min_samples_split=2, n_estimators=1000, random_state=42)
rf_tuned.fit(X_train, y_train)
rf_tuned_pred = rf_tuned.predict(X_test)
print('Random Forest Tuned')
print('Mean Squared Error:', mean_squared_error(y_test, rf_tuned_pred))
print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, rf_tuned_pred)))
print('Mean Absolute Error:', mean_absolute_error(y_test, rf_tuned_pred))
print('R2 Score:', r2_score(y_test, rf_tuned_pred))

# GradientBoostingRegressor
gb_tuned = GradientBoostingRegressor(learning_rate=0.01, max_depth=3, max_features=5, min_samples_leaf=1, min_samples_split=2, n_estimators=1000, random_state=42, subsample=1)
gb_tuned.fit(X_train, y_train)
gb_tuned_pred = gb_tuned.predict(X_test)
print('GradientBoostingRegressor Tuned')
print('Mean Squared Error:', mean_squared_error(y_test, gb_tuned_pred))
print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, gb_tuned_pred)))
print('Mean Absolute Error:', mean_absolute_error(y_test, gb_tuned_pred))
print('R2 Score:', r2_score(y_test, gb_tuned_pred))

# AdaBoostRegressor
ada_tuned = AdaBoostRegressor(learning_rate=0.01, n_estimators=1000, random_state=42)
ada_tuned.fit(X_train, y_train)
ada_tuned_pred = ada_tuned.predict(X_test)
print('AdaBoostRegressor Tuned')
print('Mean Squared Error:', mean_squared_error(y_test, ada_tuned_pred))
print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, ada_tuned_pred)))
print('Mean Absolute Error:', mean_absolute_error(y_test, ada_tuned_pred))
print('R2 Score:', r2_score(y_test, ada_tuned_pred))

# KNeighborsRegressor
knn_tuned = KNeighborsRegressor(n_neighbors=3)
knn_tuned.fit(X_train, y_train)
knn_tuned_pred = knn_tuned.predict(X_test)
print('KNeighborsRegressor Tuned')
print('Mean Squared Error:', mean_squared_error(y_test, knn_tuned_pred))
print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, knn_tuned_pred)))
print('Mean Absolute Error:', mean_absolute_error(y_test, knn_tuned_pred))
print('R2 Score:', r2_score(y_test, knn_tuned_pred))

Random Forest Tuned
Mean Squared Error: 0.45508467626236465
Root Mean Squared Error: 0.6745996414632642
Mean Absolute Error: 0.4930034936591129
R2 Score: 0.5254852218125935
GradientBoostingRegressor Tuned
Mean Squared Error: 0.29546406657813795
Root Mean Squared Error: 0.5435660645939351
Mean Absolute Error: 0.4266534944405497
R2 Score: 0.6919209252085536
AdaBoostRegressor Tuned
Mean Squared Error: 0.284297692788161
Root Mean Squared Error: 0.5331957359058314
Mean Absolute Error: 0.385970685031292
R2 Score: 0.7035640537481176
KNeighborsRegressor Tuned
Mean Squared Error: 0.4097848888888888
Root Mean Squared Error: 0.6401444281479679
Mean Absolute Error: 0.4633589743589744
R2 Score: 0.5727191096551916

Model Saving

The pickle.dump() function is used to serialize an object hierarchy and write it to a file. The first argument to pickle.dump() is the object you want to serialize, and the second argument is a file-like object, such as a file opened in write mode using the open() function. In your code, you are creating new files with the names ‘rf_tuned.pkl’, ‘gb_tuned.pkl’, ‘ada_tuned.pkl’, and ‘knn_tuned.pkl’ and opening them in write mode (‘wb’) before passing them to pickle.dump().

By serializing the models and saving them to files, you can later load them using the pickle.load() function and use them to make predictions or perform other operations. This can be useful if you want to use the models in a different script or after your Python interpreter has been closed.

import pickle

pickle.dump(rf_tuned, open('rf_tuned.pkl', 'wb'))
pickle.dump(gb_tuned, open('gb_tuned.pkl', 'wb'))
pickle.dump(ada_tuned, open('ada_tuned.pkl', 'wb'))
pickle.dump(knn_tuned, open('knn_tuned.pkl', 'wb'))

Model Loading

The pickle module is a part of the Python Standard Library and provides a way to serialize and deserialize Python objects. The pickle.load() function is used to deserialize a pickled object hierarchy stored in a file, which you can then use to make predictions or perform other operations with the model.

import pickle

rf_tuned = pickle.load(open('rf_tuned.pkl', 'rb'))
gb_tuned = pickle.load(open('gb_tuned.pkl', 'rb'))
ada_tuned = pickle.load(open('ada_tuned.pkl', 'rb'))
knn_tuned = pickle.load(open('knn_tuned.pkl', 'rb'))

Model Visualization

# Random Forest
rf_tuned_pred = rf_tuned.predict(X_test)
plt.figure(figsize=(10, 5))
plt.scatter(y_test, rf_tuned_pred, color='red')
plt.plot(y_test, y_test, color='blue')
plt.xlabel('Actual')
plt.ylabel('Predicted')
plt.title('Random Forest Tuned')
plt.show()

Conclusion

We have the results of some final model evaluations for four different models: a random forest model (rf_tuned), a gradient boosting regressor model (gb_tuned), an AdaBoost regressor model (ada_tuned), and a K-nearest neighbors regressor model (knn_tuned).

For each model, the results include several evaluation metrics: mean squared error (MSE), root mean squared error (RMSE), mean absolute error (MAE), and R2 score.

The MSE is a measure of the average squared difference between the predicted values and the true values. A lower MSE indicates that the model is making more accurate predictions. The RMSE is the square root of the MSE, which gives the error in the same units as the target variable. The MAE is the average absolute difference between the predicted values and the true values. A lower MAE indicates that the model is making more accurate predictions.

The R2 score, also known as the coefficient of determination, is a measure of how well the model explains the variance in the target variable. It ranges from 0 to 1, with a higher score indicating a better fit.

Looking at the results, it appears that the gradient boosting regressor model (gb_tuned) and the AdaBoost regressor model (ada_tuned) have the lowest MSE and RMSE, indicating that they are making more accurate predictions on average. The gradient boosting regressor model also has the highest R2 score, indicating that it is the best at explaining the variance in the target variable.

Thanks for reading this article. You can access the detailed codes of the project and other projects on my Github account or Kaggle account. Happy coding!

Please feel free to contact me if you need any further information.

References

1. https://www.kaggle.com/datasets/unsdsn/world-happiness
2. https://www.kaggle.com/code/emineyetm/happiness-around-the-world-a-comparative-analysis/notebook
3. https://worldhappiness.report/

Subscribe to DDIntel Here.

Visit our website here: https://www.datadriveninvestor.com

Join our network here: https://datadriveninvestor.com/collaborate

Insights from the World Happiness Report: A Data-Driven Perspective was originally published in DataDrivenInvestor on Medium, where people are continuing the conversation by highlighting and responding to this story.

1. Data Wrangling with Python: Tips and Tools to Make Your Life Easier- by Jacqueline Kazil and David Beazley

This book covers a wide range of data cleaning and preprocessing techniques, including handling missing values, working with strings and dates, and dealing with messy data. It also covers how to use Python’s pandas library for data manipulation.

The book begins with an introduction to Python and the fundamental concepts of data wrangling, including how to read and write data, how to explore and summarize data, and how to visualize data using Python libraries such as Pandas, NumPy, and Matplotlib.

As the book progresses, it covers more advanced topics such as how to work with large datasets, how to perform statistical analysis, and how to use machine learning algorithms to analyze data. It also includes practical examples and exercises to help readers apply what they have learned.

Overall, “Data Wrangling with Python” is a comprehensive guide to working with data in Python, and is suitable for people who are new to both programming and data analysis. It is a useful resource for anyone who wants to learn how to manipulate, analyze, and visualize data using Python

2. Bad Data Handbook: Cleaning Up The Data So You Can Get Back To Work-by Q. Ethan McCallum

The “Bad Data Handbook” is a book that provides guidance on how to identify and fix problems with data. Data quality is an important aspect of any project that involves analysis or decision-making based on data. Poor quality data can lead to incorrect conclusions and ineffective decisions, which can have serious consequences for organizations.

The book covers a variety of topics related to data quality, including how to identify bad data, how to clean and transform data, and how to maintain data quality over time. It also includes practical tips and best practices for working with data, as well as case studies and real-world examples to illustrate the concepts.

Some of the key topics covered in the “Bad Data Handbook” include:

• Identifying common problems with data, such as missing values, incorrect data types, and incorrect formatting
• Techniques for cleaning and transforming data, including using functions and scripts, using data quality tools, and working with data in a spreadsheet or database
• Strategies for maintaining data quality, including implementing data governance policies and processes, and using data quality tools and techniques to monitor and improve data quality
• Case studies and examples of how to apply these concepts in real-world situations

Overall, the “Bad Data Handbook” is a useful resource for anyone who works with data and wants to ensure that they are using high-quality data for their analysis and decision-making.

3. Feature Engineering and Selection: A Practical Approach for Predictive Models- by Kjell Johnson and Max Kuhn

“Feature Engineering and Selection: A Practical Approach for Predictive Models” is a book that teaches readers how to create and select features for use in predictive modeling. Feature engineering and selection are important aspects of the data science process, as they can significantly affect the accuracy and performance of a predictive model.

The book covers a wide range of topics related to feature engineering and selection, including how to create new features from existing data, how to select the most relevant features for a particular problem, and how to evaluate the effectiveness of different features. It also includes practical examples and case studies to illustrate the concepts and show how they can be applied in real-world situations.

Some of the key topics covered in “Feature Engineering and Selection” include:

• The role of feature engineering and selection in the data science process
• Techniques for creating and transforming features, including feature extraction, feature selection, and feature construction
• Strategies for evaluating and comparing the effectiveness of different features
• Case studies and examples of how to apply these concepts in real-world situations

Overall, “Feature Engineering and Selection” is a useful resource for anyone who is interested in predictive modeling and wants to learn more about how to create and select effective features for their models.

4. Feature Engineering for Machine Learning: Principles and Techniques for Data Scientists- by Alice Zheng and Amanda Casari

“Feature Engineering for Machine Learning” by Alice Zheng and Amanda Casari is a comprehensive guide to the principles and techniques of feature engineering for machine learning.

The book covers a wide range of topics, including:

• The role of feature engineering in the machine learning process
• Techniques for preprocessing and cleaning data
• Encoding and transforming features
• Feature selection and dimensionality reduction
• Creating new features through feature construction and feature transformation
• Working with different types of data, including numerical, categorical, and text data

The book also includes practical examples and case studies to illustrate the concepts and techniques discussed. It is suitable for both beginner and experienced data scientists who want to learn more about feature engineering and how to effectively preprocess and transform data for machine learning.

5. Python Feature Engineering Cookbook: Over 70 Recipes for Creating, Engineering, and Transforming Features to Build Machine Learning Models- by Soledad Galli

Feature engineering is an important step in the machine learning process, as it involves creating and transforming features to make them more suitable for building effective models. The “Python Feature Engineering Cookbook” by Soledad Galli offers a collection of over 70 recipes for creating, engineering, and transforming features in Python.

Some of the topics covered in the book include:

• Working with numerical and categorical features
• Handling missing values
• Encoding categorical features
• Extracting features from text
• Feature selection and dimensionality reduction
• Creating new features through feature transformation and feature construction

By following the recipes in the book, you can learn how to effectively preprocess and transform your data to improve the performance of your machine learning models. The book also includes practical tips and best practices for feature engineering, making it a useful resource for both beginner and experienced data scientists.

6. Best Practices in Data Cleaning: A Complete Guide to Everything You Need to Do Before and After Collecting Your Data- by Jason W. Osborne

“Best Practices in Data Cleaning: A Complete Guide to Everything You Need to Do Before and After Collecting Your Data” by Jason W. Osborne is a comprehensive guide to the process of data cleaning, which is an important step in the data preparation process for machine learning and data analysis.

The book covers a wide range of topics, including:

• The importance of data cleaning and the impact of dirty data on the accuracy of machine learning models
• Techniques for identifying and correcting errors and inconsistencies in data
• Strategies for handling missing values and outliers
• Best practices for data cleaning, including planning, documentation, and testing
• Tips for working with different types of data, including numerical, categorical, and text data

By following the best practices and techniques outlined in the book, you can ensure that your data is clean and ready for analysis or machine learning. The book is suitable for both beginner and experienced data scientists who want to learn more about data cleaning and how to effectively prepare data for analysis and machine learning.

🌸Thanks for reading this article. You can access my projects on my Github account or Kaggle account. Happy coding!

Please feel free to contact me if you need any further information.🌸