regression model outperform every models

Question

I followed from this question.

---

• Case1: I have the following task: Train for consecutive 3 days to predict each fourth day. Each day's data represents one CSV file, which has dimensions 24x25. Each datapoint of each CSV file is pixels. In this case, I want to use Ridge() regression, LinearRegression(), ARIMA(), and the LSTM() model. We want to have 90 mse for each model(more detail in my linked previous question).

• Case2: Here I am using, what is known as "the" naive forecast, or a "random walk" forecast. It is often hard to beat. The naive approach is: The guess for any day is simply the map of the previous day. Here also We want to have 90 mse for this naive approach(more detail in my linked previous question).

My code is:

import os
import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
from keras.models import Sequential
from keras.layers import LSTM, Dense
from statsmodels.tsa.arima.model import ARIMA
import warnings
from time import time

# Suppress warnings from ARIMA
warnings.filterwarnings("ignore")

# Paths
data_folder = r'C:\Users\alokj\OneDrive\Desktop\jupyter_proj\All_data'
output_folder = r'C:\Users\alokj\OneDrive\Desktop\jupyter_proj\90_days_merged'

# Ensure the output folder exists
os.makedirs(output_folder, exist_ok=True)

# List all CSV files in the folder
csv_files = [f for f in os.listdir(data_folder) if f.endswith('.csv')]
csv_files = sorted(csv_files, key=lambda x: int(x.split('_Day')[1].split('_')[0]))

# Prepare data
data_list = [pd.read_csv(os.path.join(data_folder, file), header=None).values for file in csv_files]
data_array = np.array(data_list)  # Shape: (num_days, 24, 25)

# Flatten the data for easier handling
num_days, rows, cols = data_array.shape
data_flattened = data_array.reshape(num_days, -1)  # Shape: (num_days, 600)

# Normalize the data using MinMaxScaler
scaler = MinMaxScaler()
data_flattened_scaled = scaler.fit_transform(data_flattened)

# Prepare data for model inputs
X = np.array([data_flattened_scaled[i-3:i].flatten() for i in range(3, 93)])  # Shape: (90, 1800)
y = data_flattened_scaled[3:93]  # Target is the 4th day in each sequence, shape: (90, 600)

# ---------- Linear Regression ----------
lr_model = LinearRegression()
lr_model.fit(X, y)
y_pred_lr = lr_model.predict(X)
residuals_lr = [np.mean((y[i] - y_pred_lr[i]) ** 2) for i in range(len(y))]

# ---------- Ridge Regression ----------
ridge_model = Ridge()
ridge_model.fit(X, y)
y_pred_ridge = ridge_model.predict(X)
residuals_ridge = [np.mean((y[i] - y_pred_ridge[i]) ** 2) for i in range(len(y))]

# ---------- LSTM Model ----------
X_lstm = np.array([data_flattened_scaled[i-3:i] for i in range(3, 93)])  # Shape: (90, 3, 600)
y_lstm = data_flattened_scaled[3:93].reshape(len(y), 1, 600)  # Shape: (90, 1, 600)

lstm_model = Sequential()
lstm_model.add(LSTM(64, activation='relu', input_shape=(X_lstm.shape[1], X_lstm.shape[2])))
lstm_model.add(Dense(600))
lstm_model.compile(optimizer='adam', loss='mse')
lstm_model.fit(X_lstm, y_lstm, epochs=50, batch_size=8, verbose=1)

y_pred_lstm = lstm_model.predict(X_lstm)
residuals_lstm = [np.mean((y_lstm[i][0] - y_pred_lstm[i]) ** 2) for i in range(len(y_lstm))]

# ---------- Naive Prediction ----------
residuals_naive = [np.mean((data_flattened_scaled[i+1] - data_flattened_scaled[i]) ** 2) for i in range(90)]

# ---------- ARIMA Model ----------
residuals_arima = []
X_arima = np.array([data_flattened_scaled[i-3:i] for i in range(3, 93)])  # Shape: (90, 3, 600)

# Timing ARIMA model
start_time = time()

for day in range(X_arima.shape[0]):  # Loop over 90 sequences
    train_data = X_arima[day]        # Shape: (3, 600)
    test_data = y[day]              # Shape: (600)
    day_predictions = []

    for feature in range(train_data.shape[1]):  # Loop over 600 features
        feature_train = train_data[:, feature]  # Shape: (3,)
        model = ARIMA(feature_train, order=(1, 0, 0))  # ARIMA(1, 0, 0)
        arima_model = model.fit()

        # Forecast for the next time step
        forecast = arima_model.forecast(steps=1)[0]
        day_predictions.append(forecast)

    # Calculate MSE for the day's predictions
    residuals_arima.append(np.mean((test_data - np.array(day_predictions)) ** 2))
    print(f"Day {day + 1}/{X_arima.shape[0]} completed in {time() - start_time:.2f}s")

# ---------- Plotting ----------
days = [f'Day {i+1}' for i in range(90)]  # Labels for Days 4 to 93

plt.figure(figsize=(12, 6))
plt.plot(days, residuals_lr, label='Linear Regression Residuals', marker='o')
plt.plot(days, residuals_ridge, label='Ridge Regression Residuals', marker='o', color='orange')
plt.plot(days, residuals_lstm, label='LSTM Residuals', marker='o', color='green')
plt.plot(days, residuals_naive, label='Naive Prediction Residuals', marker='o', color='red')
plt.plot(days, residuals_arima, label='ARIMA Residuals', marker='o', color='purple')

# Configure plot
plt.xticks(rotation=90)
plt.xlabel('Predicted Day')
plt.ylabel('Residuals (MSE)')
plt.title('Residuals for Predicted Days (Linear, Ridge, LSTM, Naive, ARIMA)')
plt.legend()
plt.grid(True)

# Save and show plot
plt.savefig(os.path.join(output_folder, 'residuals_plot_combined.png'))
plt.show()

My result is:

Here we see that the regression model outperforms every model, would anybody check all the models inside the code?

My all 90 days data folder link that I used for code.

Answer 1

OK. It's taken me longer than I expected, and I wish I had more time to put finishing touches on, and check for errors. Feel free to point out any errors :) I hope it is of some help :)

Code Review and Suggestions

Strengths

• Comprehensive Approach:

  • The code compares multiple methods (Linear Regression, Ridge, LSTM, ARIMA, and Naive Prediction), providing a well-rounded view of model performance.
  • Includes both traditional statistical models (ARIMA) and modern machine learning techniques (LSTM), demonstrating versatility.

• Clear Organisation:

  • Logical sections for preprocessing, model training, evaluation, and plotting make the code easy to follow and maintain.
  • Each model is implemented separately, enabling modular testing and debugging.

• Scalability Considerations:

  • The data is preprocessed using scaling, which ensures compatibility with machine learning models and avoids numerical instability.
  • Flattening the input arrays into vectors prepares the dataset for models requiring tabular input, demonstrating adaptability to different architectures.

• Error Handling:

  • Suppresses ARIMA warnings, keeping the output cleaner and more readable for the user.

• Reusable Framework:

  • The code is structured to handle time-series predictions systematically, making it reusable for similar tasks or datasets.
  • Methods for extracting features (eg., consecutive days for training) and targets are general enough to be adapted to other sequences or window sizes.

• Data Handling:

  • Works with multi-dimensional time-series data (24x25 grids), demonstrating the ability to manage and process spatial-temporal datasets effectively.
  • Scales data appropriately for machine learning models, ensuring consistent preprocessing across approaches.

---

Areas for Improvement

The following suggestions are listed (roughly) in decreasing order of seriousness.

1. Data Splitting and Scaling

• Issue: The MinMaxScaler is applied to the entire dataset before splitting into training and test sets. This introduces data leakage because information from the test set (eg., its minimum and maximum values) is used to normalise the training data. This artificially inflates model performance, as the test set is no longer a truly unseen dataset.

• Why This Is a Problem:

  • Data leakage occurs when information from outside the training data influences the model during training. In this case, knowing the test set's range (minimum and maximum values) allows the model to "cheat" by benefiting from test set information during training.
  • For example, if the test set contains extreme values absent from the training set, scaling the entire dataset compresses the range of the training data unnaturally, leading to overly optimistic evaluation.
  • For time series data, preserving temporal order is critical. The test set must simulate future, unseen data, and any preprocessing must rely solely on past (training) data.

• Fix:

    scaler.fit(X_train)  # Apply to training data only
    X_train_scaled = scaler.transform(X_train)
    X_test_scaled = scaler.transform(X_test)

• Why This Fix Works:
  • By fitting the MinMaxScaler only on the training data, the test set remains unseen during scaling. This helps to ensure that model evaluation simulates predicting future data, preserving the test set's integrity as a truly unseen dataset.

---

2. Train-Test Split and Validation

• Issue: Models are evaluated on the same dataset they are trained on, invalidating residuals as a performance metric.

• Why This Is a Problem:

  • Evaluating models on training data leads to overfitting and an overly optimistic view of performance. The residuals fail to reflect how the model generalises to unseen data.

• Fix:

    train_size = int(0.8 * len(data_flattened_scaled))
    X_train = X[:train_size]
    y_train = y[:train_size]
    X_test = X[train_size:]
    y_test = y[train_size:]

• Why This Fix Works:
  • Separating the data into training and validation sets ensures that model evaluation reflects performance on unseen data, providing an unbiased assessment.

---

3. Time Series Cross-Validation

• Issue: The code does not use appropriate cross-validation for time series, such as TimeSeriesSplit. This risks information leakage and overfitting by not respecting the temporal order of data.

• Why This Is a Problem:

  • Time series data has inherent temporal dependencies; future observations must not influence model training on past data.
  • Standard cross-validation methods (eg., k-fold) can allow future data to leak into the training process, invalidating model evaluation.

• Fix:

    from sklearn.model_selection import TimeSeriesSplit
    tscv = TimeSeriesSplit(n_splits=5)
    for train_index, test_index in tscv.split(X):
        X_train, X_test = X[train_index], X[test_index]
        y_train, y_test = y[train_index], y[test_index]

• Why This Fix Works:
  • TimeSeriesSplit helps to ensure that training data always precedes validation data, preserving the temporal order. This prevents future data from influencing model training, leading to a more realistic evaluation.

---

4. LSTM Input and Output

• Issue: This is a curious one! The LSTM target (y_lstm) is reshaped to (90, 1, 600), which looks to be inconsistent with the intended regression task.

• Why This Is a Problem:

  • LSTM output should match the dimensionality of the target variable ((90, 600) in this case) to ensure compatibility with regression.
  • LSTM models expect 3D input tensors (samples, timesteps, features), so incorrect reshaping can misalign input-output dimensions, causing training or prediction errors.

• Fix:

    y_lstm = data_flattened_scaled[3:93]  # Shape: (90, 600)

• Why This Fix Works:
  • Ensuring that the target shape matches the output dimension allows the model to perform regression correctly, aligning predictions with the target. Although the tensor shape for y_lstm may be incorrect and we would typically expect a runtime error due to mismatched dimensions, the code likely runs (I am assuming that it does, I haven't actually run it myself) because Keras can internally reshape or broadcast the target tensor to fit the model's output during training. However, this silent handling of shape mismatches can lead to unintended behaviour: the model may not optimise correctly for the intended regression task, and its predictions might not align directly with the flattened day-wise targets. Adjusting the target shape ensures consistency with the regression, allowing the model to produce outputs that are directly comparable to the ground truth and easier to interpret. I won't rule out the possibility that I have got this one wrong!

---

5. Naive Prediction

• Issue: The naive method calculates residuals based on consecutive days without accounting for test-set separation.

• Why This Is a Problem:

  • Without separating training and test sets, naive predictions fail to simulate real-world forecasting scenarios, leading to inaccurate comparisons.

• Fix: Ensure naive predictions are evaluated only on validation data:

    residuals_naive = [np.mean((X_test[i] - X_test[i - 1]) ** 2) for i in range(1, len(X_test))]

• Why This Fix Works:
  • Evaluating naive predictions on the validation set helps to ensure a fair comparison with other models, reflecting real-world forecasting challenges.

---

6. ARIMA Loop Efficiency

• Issue: The ARIMA model is fit independently for each feature of each day, leading to inefficiency and excessive computation time.

• Why This Is a Problem:

• Fitting ARIMA independently for 600 features per day is computationally expensive and redundant, especially when many features may be correlated or irrelevant.

• This approach will likely not scale well for larger datasets.

• If there are no current performance issues and you don't intend to run it a bigger scale, then this is a non-issue.

• Fix:

    from sklearn.decomposition import PCA
    pca = PCA(n_components=10)  # Choose appropriate number of components
    X_pca = pca.fit_transform(data_flattened_scaled)

• Why This Fix Works:
• Reducing dimensionality via PCA aggregates feature information, significantly decreasing computational overhead while retaining most of the variance.

---

7. Visualisation - Ambiguous Labels

• Issue: The x-axis labels (eg., Day 1, Day 2) do not clarify the relationship between predicted days and residuals.

• Why This Is a Problem:

  • Ambiguous labels make it harder for the audience to understand the results. Without clarity, viewers might misinterpret the residuals' relationship to the data (eg., thinking residuals for Day 1 reflect errors on training data when they actually reflect predictions for test data).

• Fix:

    plt.title('Residuals for Models (Validation Set)')
    plt.xlabel('Days (Validation Set)')

• Why This Fix Works:
  • Including a phrase like "Validation Set" explicitly ties the residuals to a specific data split (validation) rather than leaving it ambiguous. This tells the audience that these are errors from predictions on unseen data, not training data.

---

8. Visualisation - Overlapping Labels

• Issue: The x-axis labels (eg., Day 1, Day 2, etc.) are printed too close together, causing them to overlap and making them illegible.

• Why This Is a Problem:

  • Overlapping labels make it difficult for the audience to identify individual days or interpret the residuals associated with each day.
  • This diminishes the clarity of the visualisation and can lead to frustration or misinterpretation of the data.

• Fix:

    plt.xticks(ticks=range(0, len(days), 2), labels=[f'Day {i}' for i in range(0, len(days), 2)], rotation=45, ha='right')

• Why This Fix Works:
  • Reducing the number of labels by skipping every 2 or more days, combined with rotating and aligning the remaining labels to the right, ensures clarity and readability while preserving the overall layout and scale of the plot. You may need to experiment with the number of skipped days to achieve an optimal plot for your results.

---

Are Residuals the Best Way to Visually Compare Models?

Plotting just the residuals is a useful starting point for comparing models, as they highlight the magnitude and distribution of errors. However, residual plots alone will likely not provide a comprehensive comparison:

• Limitations: Residuals focus only on individual errors and do not summarise overall performance (eg., MSE or MAE). They also lack insight into how well models capture underlying patterns or trends in the data.

• Alternatives: Consider complementing residual plots with the following visualisations:

1. Prediction vs. Actual plot:

Overlay model predictions with actual values to assess how well the model captures trends and patterns:

import matplotlib.pyplot as plt

plt.plot(y_test, label='Ground Truth')
plt.plot(y_pred, label='Predictions', alpha=0.7)
plt.title('Predictions vs. Ground Truth')
plt.legend()
plt.show()

2. Residual Histograms:

Examine the distribution of residuals to detect systematic biases or outliers:

import numpy as np

residuals = y_test - y_pred
plt.hist(residuals, bins=20, edgecolor='black')
plt.title('Residual Distribution')
plt.xlabel('Residuals')
plt.ylabel('Frequency')
plt.show()

3. Heatmaps:

Use heatmaps to visualise residuals across time and features, revealing patterns that might otherwise go unnoticed:

import seaborn as sns
import numpy as np

# Example residual matrix
residual_matrix = np.random.rand(10, 25)  # Replace with actual residuals
sns.heatmap(residual_matrix, annot=True, fmt=".2f", cmap="coolwarm")
plt.title("Residuals Heatmap")
plt.xlabel("Features")
plt.ylabel("Days")
plt.show()

Statistical Methods for Model Comparison

While visualisations provide qualitative insights into model performance, statistical methods offer a formal framework to compare models quantitatively. These methods ensure that observed differences in performance are meaningful and not due to random variation. Below are some common statistical approaches for model comparison:

---

1. Aggregate Metrics

Aggregate metrics summarise model performance across all predictions, providing a clear quantitative measure of accuracy. Common metrics include:

• Mean Absolute Error (MAE): Measures the average magnitude of errors in predictions, without considering their direction.
• Mean Squared Error (MSE): Penalises larger errors more heavily, making it sensitive to outliers.

Example implementation:

from sklearn.metrics import mean_squared_error, mean_absolute_error

mse = mean_squared_error(y_test, y_pred)
mae = mean_absolute_error(y_test, y_pred)
print(f"MSE: {mse:.2f}, MAE: {mae:.2f}")

2. Likelihood-Based Criteria (AIC/BIC):

For models based on likelihood estimation (eg., ARIMA), metrics such as:

Akaike Information Criterion (AIC): Balances goodness-of-fit with model complexity by penalising the number of parameters. Bayesian Information Criterion (BIC): Similar to AIC but imposes a stricter penalty on model complexity. These criteria allow for the comparison of models by evaluating their relative likelihoods while discouraging overfitting.

3. Diebold-Mariano Test:

In time-series contexts, the Diebold-Mariano test is often used to compare the predictive accuracy of two models. It tests whether the forecast errors from two models differ significantly, accounting for potential autocorrelation in the residuals.

This test is especially valuable when comparing models that produce sequential predictions over time.

4. Residual-Based Tests:

For general regression tasks, residual-based tests can help assess performance differences:

Paired t-tests: Compare the absolute, or squared, residuals from two models to determine if one produces smaller errors than the other.

ANOVA-based tests: Analyse residual variance across multiple models to identify differences.

Summing up

This review highlights both the strengths of the current implementation and the areas where refinements can improve the robustness and interpretability of results. The code demonstrates a comprehensive approach by combining traditional statistical models like ARIMA with modern machine learning methods such as LSTM, showcasing versatility in addressing complex time-series problems. However, several improvements, such as addressing data leakage, adopting proper time-series cross-validation, and refining visualisations techniques can improve the reliability and clarity of the analysis.

In particular, ensuring that the data pipeline preserves temporal integrity (eg., through proper train-test splits and scaling) should avoid inflated performance metrics. Exploring alternative visualisation methods, like heatmaps and prediction vs. actual plots, can provide additional insights into model behaviour. Complementing these with statistical methods, such as aggregate metrics, AIC/BIC, and the Diebold-Mariano test, allows for rigorous quantitative comparisons between models.

The unusual results reported in the OP likely stem from some of the identified issues, such as improper data splitting or mismatched model evaluation techniques. Implementing the suggested fixes and enhancements will not only address these concerns but also establish a robust framework for fair and meaningful model comparison. This foundation can be extended to handle larger datasets, more sophisticated models, and deeper statistical insights, making the analysis both scalable and adaptable for future applications.