Regression Prediction

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import warnings
warnings.filterwarnings('ignore')
 
# Set seed for reproducible results
np.random.seed(42)
 
# Generate customer feature data
n_customers = 800
 
customer_features = pd.DataFrame({
    'user_id': range(1, n_customers + 1),
    'total_orders': np.random.poisson(5, n_customers) + 1,
    'total_items': np.random.poisson(15, n_customers) + 1,
    'avg_order_value': np.random.exponential(80, n_customers) + 20,
    'order_std': np.random.exponential(30, n_customers),
    'tenure_days': np.random.randint(30, 730, n_customers),
    'avg_order_gap': np.random.exponential(30, n_customers) + 5,
    'unique_categories': np.random.randint(1, 10, n_customers),
    'unique_brands': np.random.randint(1, 15, n_customers)
})
 
# Generate CLV (target) - with relationship to features
customer_features['total_spent'] = (
    customer_features['total_orders'] * customer_features['avg_order_value'] +
    np.random.normal(0, 100, n_customers)
).clip(50, None)
 
# Handle missing values
customer_features = customer_features.fillna(0)
 
print(f"Number of customers: {len(customer_features)}")
print(f"Average CLV: ${customer_features['total_spent'].mean():,.2f}")
print(f"Median CLV: ${customer_features['total_spent'].median():,.2f}")
print(f"CLV range: ${customer_features['total_spent'].min():,.2f} ~ ${customer_features['total_spent'].max():,.2f}")

# Separate features and target
feature_cols = ['total_orders', 'total_items', 'avg_order_value', 'order_std',
                'tenure_days', 'avg_order_gap', 'unique_categories', 'unique_brands']
X = customer_features[feature_cols]
y = customer_features['total_spent']
 
# Train/test split
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)
 
# Scaling
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
 
print(f"Training set: {len(X_train)} samples")
print(f"Test set: {len(X_test)} samples")
print(f"Training CLV mean: ${y_train.mean():,.2f}")
print(f"Test CLV mean: ${y_test.mean():,.2f}")

y = β₀ + β₁x₁ + β₂x₂ + ... + βₙxₙ + ε

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
 
# Model training
lr_model = LinearRegression()
lr_model.fit(X_train_scaled, y_train)
 
# Prediction
y_pred_lr = lr_model.predict(X_test_scaled)
 
# Evaluation
print("=== Linear Regression Results ===")
print(f"MAE: ${mean_absolute_error(y_test, y_pred_lr):,.2f}")
print(f"RMSE: ${np.sqrt(mean_squared_error(y_test, y_pred_lr)):,.2f}")
print(f"R²: {r2_score(y_test, y_pred_lr):.3f}")

import matplotlib.pyplot as plt
 
# Coefficients by feature
coef_df = pd.DataFrame({
    'feature': feature_cols,
    'coefficient': lr_model.coef_
}).sort_values('coefficient', key=abs, ascending=False)
 
print("\nCoefficients by Feature (Influence):")
print(coef_df.to_string(index=False))
 
# Visualization
plt.figure(figsize=(10, 6))
colors = ['green' if c > 0 else 'red' for c in coef_df['coefficient']]
plt.barh(coef_df['feature'], coef_df['coefficient'], color=colors)
plt.xlabel('Coefficient')
plt.title('Linear Regression Feature Coefficients', fontsize=14, fontweight='bold')
plt.axvline(x=0, color='black', linestyle='-', linewidth=0.5)
plt.tight_layout()
plt.show()
 
# Interpretation example
top_feature = coef_df.iloc[0]['feature']
top_coef = coef_df.iloc[0]['coefficient']
print(f"\nInterpretation: When {top_feature} increases by 1 standard deviation, CLV changes by ${top_coef:,.2f}")

from sklearn.linear_model import Ridge
 
# Test multiple alpha values
alphas = [0.01, 0.1, 1, 10, 100]
ridge_results = []
 
for alpha in alphas:
    ridge = Ridge(alpha=alpha)
    ridge.fit(X_train_scaled, y_train)
    y_pred = ridge.predict(X_test_scaled)
    r2 = r2_score(y_test, y_pred)
    ridge_results.append({'alpha': alpha, 'r2': r2})
 
ridge_df = pd.DataFrame(ridge_results)
print("Ridge Regression R² by alpha:")
print(ridge_df.to_string(index=False))
 
# Train model with optimal alpha
best_alpha = ridge_df.loc[ridge_df['r2'].idxmax(), 'alpha']
ridge_model = Ridge(alpha=best_alpha)
ridge_model.fit(X_train_scaled, y_train)
y_pred_ridge = ridge_model.predict(X_test_scaled)
 
print(f"\nOptimal alpha: {best_alpha}")
print(f"Ridge R²: {r2_score(y_test, y_pred_ridge):.3f}")

from sklearn.linear_model import Lasso
 
# Lasso regression
lasso_model = Lasso(alpha=0.1, max_iter=10000)
lasso_model.fit(X_train_scaled, y_train)
y_pred_lasso = lasso_model.predict(X_test_scaled)
 
# Selected features (non-zero coefficients)
selected_features = pd.DataFrame({
    'feature': feature_cols,
    'coefficient': lasso_model.coef_
})
selected_features = selected_features[selected_features['coefficient'] != 0]
print(f"Lasso Selected Features ({len(selected_features)}):")
print(selected_features.to_string(index=False))
 
print(f"\nLasso R²: {r2_score(y_test, y_pred_lasso):.3f}")

from sklearn.ensemble import RandomForestRegressor
 
# Model training
rf_model = RandomForestRegressor(
    n_estimators=100,
    max_depth=10,
    min_samples_split=10,
    random_state=42,
    n_jobs=-1
)
rf_model.fit(X_train, y_train)
 
# Prediction (no scaling needed)
y_pred_rf = rf_model.predict(X_test)
 
# Evaluation
print("=== Random Forest Regression Results ===")
print(f"MAE: ${mean_absolute_error(y_test, y_pred_rf):,.2f}")
print(f"RMSE: ${np.sqrt(mean_squared_error(y_test, y_pred_rf)):,.2f}")
print(f"R²: {r2_score(y_test, y_pred_rf):.3f}")

# Feature importance
importance_df = pd.DataFrame({
    'feature': feature_cols,
    'importance': rf_model.feature_importances_
}).sort_values('importance', ascending=False)
 
print("Feature Importance:")
print(importance_df.to_string(index=False))
 
# Visualization
plt.figure(figsize=(10, 6))
plt.barh(importance_df['feature'], importance_df['importance'], color='forestgreen')
plt.xlabel('Importance')
plt.title('Random Forest Feature Importance', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

from xgboost import XGBRegressor
 
# Model training
xgb_model = XGBRegressor(
    n_estimators=100,
    max_depth=6,
    learning_rate=0.1,
    subsample=0.8,
    colsample_bytree=0.8,
    random_state=42
)
xgb_model.fit(X_train, y_train)
 
# Prediction
y_pred_xgb = xgb_model.predict(X_test)
 
# Evaluation
print("=== XGBoost Regression Results ===")
print(f"MAE: ${mean_absolute_error(y_test, y_pred_xgb):,.2f}")
print(f"RMSE: ${np.sqrt(mean_squared_error(y_test, y_pred_xgb)):,.2f}")
print(f"R²: {r2_score(y_test, y_pred_xgb):.3f}")

# Model performance comparison
models = {
    'Linear Regression': y_pred_lr,
    'Ridge': y_pred_ridge,
    'Lasso': y_pred_lasso,
    'Random Forest': y_pred_rf,
    'XGBoost': y_pred_xgb
}
 
results = []
for name, y_pred in models.items():
    results.append({
        'Model': name,
        'MAE': mean_absolute_error(y_test, y_pred),
        'RMSE': np.sqrt(mean_squared_error(y_test, y_pred)),
        'R²': r2_score(y_test, y_pred)
    })
 
results_df = pd.DataFrame(results).round(2)
print("=== Model Performance Comparison ===")
print(results_df.to_string(index=False))

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
best_models = [('Linear Regression', y_pred_lr),
               ('Random Forest', y_pred_rf),
               ('XGBoost', y_pred_xgb)]
 
for ax, (name, y_pred) in zip(axes, best_models):
    ax.scatter(y_test, y_pred, alpha=0.5, s=20)
    ax.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()],
            'r--', linewidth=2, label='Perfect Prediction')
    ax.set_xlabel('Actual CLV ($)')
    ax.set_ylabel('Predicted CLV ($)')
    ax.set_title(f'{name}\nR² = {r2_score(y_test, y_pred):.3f}')
    ax.legend()
 
plt.tight_layout()
plt.show()

# Residual analysis (based on best model)
residuals = y_test - y_pred_xgb
 
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
 
# Residual distribution
axes[0].hist(residuals, bins=30, edgecolor='black', alpha=0.7)
axes[0].axvline(x=0, color='red', linestyle='--')
axes[0].set_xlabel('Residual ($)')
axes[0].set_ylabel('Frequency')
axes[0].set_title('Residual Distribution')
 
# Residual vs Predicted
axes[1].scatter(y_pred_xgb, residuals, alpha=0.5, s=20)
axes[1].axhline(y=0, color='red', linestyle='--')
axes[1].set_xlabel('Predicted CLV ($)')
axes[1].set_ylabel('Residual ($)')
axes[1].set_title('Residual vs Predicted')
 
plt.tight_layout()
plt.show()
 
# Residual statistics
print(f"Residual Mean: ${residuals.mean():,.2f}")
print(f"Residual Std: ${residuals.std():,.2f}")

from sklearn.model_selection import cross_val_score
 
# 5-Fold cross-validation
cv_scores = cross_val_score(
    xgb_model, X, y,
    cv=5,
    scoring='r2'
)
 
print("=== 5-Fold Cross-Validation ===")
print(f"R² Scores: {cv_scores.round(3)}")
print(f"Mean R²: {cv_scores.mean():.3f} (+/- {cv_scores.std():.3f})")

from sklearn.model_selection import GridSearchCV
 
# Grid search
param_grid = {
    'n_estimators': [50, 100, 200],
    'max_depth': [4, 6, 8],
    'learning_rate': [0.05, 0.1, 0.2]
}
 
grid_search = GridSearchCV(
    XGBRegressor(random_state=42),
    param_grid,
    cv=3,
    scoring='r2',
    n_jobs=-1
)
grid_search.fit(X_train, y_train)
 
print("Best Parameters:", grid_search.best_params_)
print(f"Best R²: {grid_search.best_score_:.3f}")