Beyond the Straight Line: Choosing Between OLS, Interaction Terms, and Tweedie Regression

Beyond the Straight Line: Choosing Between OLS, Interaction Terms, and Tweedie Regression

Linear regression is often the first tool data scientists reach for when modeling a continuous outcome. But real-world data rarely follows a straight line. Relationships can be curved, heterogeneous, or fraught with zeros and skewness. In this article, we explore three increasingly sophisticated approaches: ordinary least squares (OLS), OLS with interaction terms, and Tweedie regression. We'll compare their assumptions, strengths, and pitfalls, then show you how to implement each in Python with real commands.

The Limits of Ordinary Least Squares

Ordinary least squares regression assumes a linear relationship between predictors and the response, normally distributed errors with constant variance, and independence of observations. When these assumptions hold, OLS is the best linear unbiased estimator. But in practice, many datasets violate them.

Consider predicting insurance claim costs. Claim amounts are strictly positive, heavily right-skewed, and often contain many zeros (no claim). OLS would predict negative values, fail to capture the zero mass, and produce heteroscedastic residuals. Similarly, in marketing mix modeling, the effect of advertising spend may be nonlinear—diminishing returns at high spend levels. OLS would miss this curvature.

The first step beyond the straight line is to add interaction terms, which allow the effect of one predictor to depend on another.

Interaction Terms: Capturing Conditional Relationships

Interaction terms model how the relationship between a predictor \(X_1\) and the outcome \(Y\) changes with the value of another predictor \(X_2\). For example, the effect of education on income might be stronger for younger workers than older workers. You can model this by including the product \(X_1 \times X_2\) in the regression.

Interaction terms are intuitive and easy to implement. They allow slopes to vary across groups or levels of a continuous variable. However, they still assume a linear relationship between the predictors and the outcome after accounting for interactions. They also require careful interpretation: the main effects change meaning when an interaction is present.

For data with severe skewness, zero inflation, or non-constant variance, interactions alone are insufficient. You need a distributional model that respects the data's nature. This is where Tweedie regression shines.

Tweedie Regression: Handling Skewed, Zero-Inflated, and Continuous Data

Tweedie regression belongs to the family of generalized linear models (GLMs). It is designed for data that are continuous but have a non-zero probability of being exactly zero, and are positively skewed. This combination makes it ideal for insurance claims, healthcare costs, retail sales, and other "semi-continuous" outcomes.

The Tweedie distribution has a power parameter \(p\) (between 1 and 2) that determines the distribution's shape:

• \(p = 1\): Poisson (count data, discrete)
• \(p = 2\): Gamma (continuous, positive, skewed)
• \(1 < p < 2\): Compound Poisson–Gamma (continuous with zeros)

By choosing the appropriate \(p\), you can model data that have a spike at zero and a long right tail. Tweedie regression also handles heteroscedasticity naturally because the variance is a function of the mean.

Requirements

Before we dive into code, ensure your environment has the following Python packages. We'll use `statsmodels` for OLS and interaction terms, and `scikit-learn` for Tweedie regression (via `TweedieRegressor`).

• Python 3.8 or later
• `numpy`
• `pandas`
• `statsmodels`
• `scikit-learn`
• `matplotlib` (for optional plotting)

Step-by-Step Installation

We'll create a dedicated virtual environment and install the required packages.

First, create and activate a virtual environment. This isolates dependencies for this project.

python3 -m venv regression_env
source regression_env/bin/activate   # On Windows: regression_env\Scripts\activate

Next, upgrade `pip` to the latest version to avoid dependency conflicts.

pip install --upgrade pip

Now install the core packages. We'll install `numpy`, `pandas`, `statsmodels`, `scikit-learn`, and `matplotlib`.

pip install numpy pandas statsmodels scikit-learn matplotlib

To verify the installation, run a quick Python check. This should print the versions without errors.

python -c "import numpy; import pandas; import statsmodels; import sklearn; print('All packages installed successfully')"

Usage Examples

We'll simulate a realistic dataset: insurance claim costs. It will have a continuous predictor `age` and a categorical predictor `vehicle_type` (0 = sedan, 1 = truck). The response `claim_amount` is zero-inflated and right-skewed.

1. Simulate the Data

First, generate the data. We'll create 1000 observations with a known structure.

import numpy as np
import pandas as pd

np.random.seed(42)
n = 1000
age = np.random.uniform(18, 70, n)
vehicle_type = np.random.binomial(1, 0.4, n)

# Base claim probability: 30% for sedan, 50% for truck
prob_claim = 0.3 + 0.2 * vehicle_type
has_claim = np.random.binomial(1, prob_claim)

# Claim amount: gamma distributed, mean depends on age and vehicle type
mean_claim = 200 + 5 * age + 50 * vehicle_type
claim_amount = has_claim * np.random.gamma(shape=2, scale=mean_claim/2)

df = pd.DataFrame({
    'age': age,
    'vehicle_type': vehicle_type,
    'claim_amount': claim_amount
})
print(df.head())

2. Fit OLS with Interaction Terms

We use `statsmodels` to fit OLS with an interaction between `age` and `vehicle_type`. Note that OLS is not appropriate for this data (zeros and skewness), but we include it for comparison.

import statsmodels.api as sm

# Add interaction term
df['age_vehicle'] = df['age'] * df['vehicle_type']
X = df[['age', 'vehicle_type', 'age_vehicle']]
X = sm.add_constant(X)
y = df['claim_amount']

ols_model = sm.OLS(y, X).fit()
print(ols_model.summary())

The output shows coefficients, standard errors, and diagnostics. Notice the R-squared might be low because OLS fails to capture the zero mass.

3. Fit Tweedie Regression

Now we fit a Tweedie regression using `scikit-learn`'s `TweedieRegressor`. We'll use `power=1.5` (typical for insurance data). We also scale the features for better convergence.

from sklearn.linear_model import TweedieRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split

# Prepare features (no interaction needed—Tweedie handles nonlinearity via link function)
X = df[['age', 'vehicle_type']]
y = df['claim_amount']

# Scale features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Split data
X_train, X_test, y_train, y_test = train_test_split(
    X_scaled, y, test_size=0.2, random_state=42
)

# Fit Tweedie regression
tweedie_model = TweedieRegressor(
    power=1.5,   # compound Poisson–Gamma
    alpha=0.0,   # no regularization
    max_iter=1000,
    tol=1e-4
)
tweedie_model.fit(X_train, y_train)

# Predict and evaluate
y_pred = tweedie_model.predict(X_test)
print("Coefficients:", tweedie_model.coef_)
print("Intercept:", tweedie_model.intercept_)

4. Compare Model Performance

We'll compare the mean absolute error (MAE) and mean squared error (MSE) between OLS and Tweedie. Note that MSE is more sensitive to large errors, which Tweedie handles better.

from sklearn.metrics import mean_absolute_error, mean_squared_error

# OLS predictions on test set
X_test_ols = sm.add_constant(pd.DataFrame({
    'age': df.loc[X_test_indices, 'age'],
    'vehicle_type': df.loc[X_test_indices, 'vehicle_type'],
    'age_vehicle': df.loc[X_test_indices, 'age'] * df.loc[X_test_indices, 'vehicle_type']
}))
# (Assume X_test_indices are the indices from train_test_split)
# For simplicity, we re-predict on full data (this is illustrative)
ols_pred = ols_model.predict(X)

# Tweedie predictions (already computed)
print("OLS MAE:", mean_absolute_error(y, ols_pred))
print("Tweedie MAE:", mean_absolute_error(y_test, y_pred))
print("OLS MSE:", mean_squared_error(y, ols_pred))
print("Tweedie MSE:", mean_squared_error(y_test, y_pred))

In practice, Tweedie will yield lower MSE because it accounts for the zero mass and skewness.

5. Visualizing the Fit

Plot the predicted vs. actual values for both models to see the difference.

import matplotlib.pyplot as plt

plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.scatter(y, ols_pred, alpha=0.3)
plt.plot([0, y.max()], [0, y.max()], 'r--')
plt.xlabel('Actual')
plt.ylabel('OLS Predicted')
plt.title('OLS: Predicted vs Actual')

plt.subplot(1, 2, 2)
plt.scatter(y_test, y_pred, alpha=0.3)
plt.plot([0, y_test.max()], [0, y_test.max()], 'r--')
plt.xlabel('Actual')
plt.ylabel('Tweedie Predicted')
plt.title('Tweedie: Predicted vs Actual')

plt.tight_layout()
plt.show()

The Tweedie plot should show less scatter and better alignment along the diagonal.

When to Use Each Approach

• **OLS**: Use only when your data is approximately normal, homoscedastic, and has no zeros or bounded range. Good for controlled experiments or standardized test scores.
• **OLS with interactions**: Use when you suspect the effect of one variable depends on another, but the data still satisfy OLS assumptions. For example, marketing ROI by channel.
• **Tweedie regression**: Use when your outcome is continuous, positive (or zero), and right-skewed. Ideal for insurance, healthcare, and sales data. It naturally handles zero inflation and heteroscedasticity.

Conclusion

Linear regression is a powerful baseline, but real-world data often demands more. Interaction terms add flexibility for conditional relationships, but they don't address distributional issues like skewness and zeros. Tweedie regression bridges the gap—it respects the data's true nature while remaining interpretable and computationally efficient. By choosing the right tool for the job, you move beyond the straight line into models that reflect the complexity of the world. Start with OLS, add interactions when needed, and switch to Tweedie when your data calls for it. Your models—and your stakeholders—will thank you.