Farm Productivity Prediction using Stepwise Regression in ML

Accurately forecasting farm productivity—measured as crop yield (hectograms per hectare)—is essential for farmers, agribusinesses, and policy‑makers to optimise input use, manage food supply, and ensure economic sustainability. In this project, we will predict farm yield based on environmental (rainfall, temperature), agronomic (pesticide usage), and temporal (year) factors. By applying stepwise regression, we aim to select the most significant predictors and build a concise, interpretable linear model that balances simplicity with predictive performance—enabling data‑driven decisions to boost productivity.

Libraries Required

import pandas as pd               # Data manipulation  
import numpy as np                # Numerical operations  
import statsmodels.api as sm      # Statistical modeling  
from sklearn.model_selection import train_test_split   # Data splitting  
from sklearn.metrics import r2_score, mean_squared_error  # Evaluation  
import matplotlib.pyplot as plt   # Visualization  

Dataset

Crop Yield Prediction Dataset

Step-by-Step Code Implementation

Data Loading & Initial Inspection

We read ~28k records of farm yield data, inspecting variable types and basic statistics to understand distributions (rainfall, temperature, pesticide usage).

# Block 1: Load dataset
url = "https://www.kaggle.com/datasets/patelris/crop-yield-prediction-dataset/download"
df = pd.read_csv(url)

# Inspect data
print(df.head())
print(df.info())
print(df.describe())

Data Preprocessing

• We drop any missing rows to simplify modelling. Features (X) include agro‑climatic inputs and Year; target (y) is Yield per Hectare. We split the data into 80% for training and 20% for testing.
• The dataset contains ~28,000 records of farm-level yields, along with annual rainfall, average temperature, and pesticide use.

# Block 2: Clean & prepare
# Handle missing values (if any)
df = df.dropna()

# Separate features and target
X = df[["Annual Rainfall (mm)", "Average Temperature (°C)", "Pesticide Usage (kg/ha)", "Year"]]
y = df["Yield per Hectare (hg/ha)"]

# Train–test split
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

Stepwise Regression Function

Our stepwise_selection function alternates forward inclusion (adding the excluded predictor with p < 0.01) and backward elimination (removing included predictors with p > 0.05) until no changes remain—yielding a parsimonious set of significant features.

# Block 3: Forward–backward stepwise selection
def stepwise_selection(X, y, 
                       initial_list=[], 
                       threshold_in=0.01, 
                       threshold_out=0.05, 
                       verbose=True):
    included = list(initial_list)
    while True:
        changed = False
        # Forward step
        excluded = list(set(X.columns) - set(included))
        new_pvals = pd.Series(index=excluded, dtype=float)
        for col in excluded:
            model = sm.OLS(y, sm.add_constant(X[included + [col]])).fit()
            new_pvals[col] = model.pvalues[col]
        best_pval = new_pvals.min()
        if best_pval < threshold_in:
            best_var = new_pvals.idxmin()
            included.append(best_var)
            changed = True
            if verbose:
                print(f"Add  {best_var:30} p-value {best_pval:.6f}")

        # Backward step
        model = sm.OLS(y, sm.add_constant(X[included])).fit()
        pvals = model.pvalues.iloc[1:]  # exclude intercept
        worst_pval = pvals.max()
        if worst_pval > threshold_out:
            worst_var = pvals.idxmax()
            included.remove(worst_var)
            changed = True
            if verbose:
                print(f"Drop {worst_var:30} p-value {worst_pval:.6f}")
        if not changed:
            break
    return included

Model Building & Evaluation

We fit an Ordinary Least Squares regression using statsmodels on the selected predictors. The .summary() provides coefficient estimates, p-values, R², and diagnostic metrics, revealing each predictor’s impact.

Predictions on the test set yield R² (explained variance) and RMSE (error scale), quantifying generalisation performance.

# Block 4: Select features
selected = stepwise_selection(X_train, y_train)

# Fit final OLS model
X_train_sel = sm.add_constant(X_train[selected])
model = sm.OLS(y_train, X_train_sel).fit()
print(model.summary())

# Predict on test set
X_test_sel = sm.add_constant(X_test[selected])
y_pred = model.predict(X_test_sel)

# Metrics
print("Test R²:", r2_score(y_test, y_pred))
print("Test RMSE:", np.sqrt(mean_squared_error(y_test, y_pred)))

Residual Diagnostics

Plotting residuals against predicted yields checks for non‑random patterns or heteroscedasticity, validating linear model assumptions.

# Block 5: Residual plot
residuals = y_test - y_pred
plt.scatter(y_pred, residuals)
plt.axhline(0, linestyle="--")
plt.xlabel("Predicted Yield (hg/ha)")
plt.ylabel("Residuals")
plt.title("Residuals vs. Predicted Yield")
plt.show()

Summary

Applying stepwise regression to the farm productivity dataset isolates the most powerful drivers—such as rainfall, temperature, and pesticide application—while ignoring less informative variables.

The final linear model achieves a robust balance of interpretability and accuracy (high R², low RMSE), offering farmers and policymakers a transparent tool to forecast yields, optimise resource use, and make informed agronomic decisions.