| import numpy as np | | | import pandas as pd | | | import matplotlib.pyplot as plt | | | import seaborn as sns | | | from sklearn.model_selection import train_test_split, GridSearchCV, cross_validate | | | from sklearn.pipeline import Pipeline | | | from sklearn.compose import ColumnTransformer, TransformedTargetRegressor | | | from sklearn.preprocessing import StandardScaler, OneHotEncoder, TargetEncoder | | | from sklearn.impute import SimpleImputer | | | from sklearn.feature_selection import mutual_info_regression | | | from sklearn.cluster import KMeans | | | from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV | | | from xgboost import XGBRegressor | | | from sklearn.metrics import mean_absolute_error, root_mean_squared_error, r2_score, median_absolute_error | | | from sklearn.inspection import permutation_importance | |
| # ========================================== | |
| # 0. LOAD DATA (Mock Example) | |
| # ========================================== | |
| # Features: 'price', 'surface_area', 'rooms', 'zip_code', 'property_type', 'latitude', 'longitude' | |
| # df = pd.read_csv('data.csv') | |
| # ========================================== | |
| # 1. EXPLORATORY DATA ANALYSIS (EDA) | |
| # ========================================== | |
| sns.set_theme(style="whitegrid") | |
| # Target Distribution | |
| plt.figure(figsize=(10, 5)) | |
| sns.histplot(data=df, x='price', kde=True, bins=50) | |
| plt.title("Target Distribution: Price") | |
| plt.show() | |
| # Geographical Price Heatmap | |
| plt.figure(figsize=(10, 8)) | |
| sns.scatterplot(x='longitude', y='latitude', hue='price', size='surface_area', | |
| sizes=(10, 200), palette='viridis', data=df, alpha=0.7) | |
| plt.title("Geographical Distribution of Prices") | | | plt.show() | | | # Multivariate Relationships (Pairplot on subset to save time) | | | sns.pairplot(df[['price', 'surface_area', 'rooms']], diag_kind='kde') | | | plt.suptitle("Pairplot of Key Numerical Features", y=1.02) | | | plt.show() | | | # Correlation Matrix | |
| plt.figure(figsize=(8, 6)) | |
| sns.heatmap(data=df.select_dtypes(include=np.number).corr(), annot=True, cmap='coolwarm', fmt=".2f") | |
| plt.title("Feature Correlation Matrix") | |
| plt.show() | |
| # ========================================== | |
| # 2. SPLIT & MANUAL OUTLIER FILTERING | | | # ========================================== | | | # We can also add a step here to remove rows with NAs for a subset of columns. Do before split to avoid crashing... | | | # If we don't have a price, we can't train or test. | | | # If we don't have a surface area, we refuse to predict. | |
| columns_to_drop_if_na = ['price', 'surface_area'] | |
| df_clean = df.dropna(subset=columns_to_drop_if_na).copy() | |
| X = df.drop('price', axis=1) | |
| y = df['price'] | |
| # Split BEFORE filtering to prevent leakage | | | X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42) | | | # Manual Outlier Filtering (Logic-based from observing scatterplots) | | | # Example: Remove properties > 500 sqm or prices that are obvious data entry errors | |
| valid_idx = (X_train['surface_area'] < 500) & (y_train > 10000) | |
| X_train = X_train[valid_idx].copy() | |
| y_train = y_train[valid_idx].copy() | |
| # ========================================== | |
| # 3. FEATURE ENGINEERING | | | # ========================================== | | | # Spatial Clustering (Fit on train, transform both) | |
| kmeans = KMeans(n_clusters=8, random_state=42, n_init=10) | |
| X_train['geo_cluster'] = kmeans.fit_predict(X_train[['latitude', 'longitude']]).astype(str) | |
| X_test['geo_cluster'] = kmeans.predict(X_test[['latitude', 'longitude']]).astype(str) | |
| # ========================================== | |
| # 4. PREPROCESSING | | | # ========================================== | | | # Group features by the exact mathematical transformation they require | | | features_to_scale = ['surface_area', 'rooms', 'latitude', 'longitude'] | |
| features_to_ohe = ['property_type', 'geo_cluster'] | |
| features_to_target_encode = ['zip_code'] | |
| features_to_passthrough = ['has_balcony'] # Boolean flags that need no transformation | | | # Pipeline for features that need scaling (with fallback median imputation) | |
| scale_pipeline = Pipeline([ | |
| ('imputer', SimpleImputer(strategy='median')), | |
| ('scaler', StandardScaler()) | |
| ]) | |
| # Pipeline for OHE (with fallback mode imputation) | |
| ohe_pipeline = Pipeline([ | |
| ('imputer', SimpleImputer(strategy='most_frequent')), | |
| ('ohe', OneHotEncoder(handle_unknown='ignore', sparse_output=False)) | |
| ]) | |
| # Pipeline for Target Encoding | |
| target_enc_pipeline = Pipeline([ | |
| ('imputer', SimpleImputer(strategy='most_frequent')), | |
| ('target_enc', TargetEncoder(smooth="auto")) | |
| ]) | |
| preprocessor = ColumnTransformer( | |
| transformers=[ | |
| ('scale', scale_pipeline, features_to_scale), | | | ('ohe', ohe_pipeline, features_to_ohe), | | | ('target_enc', target_enc_pipeline, features_to_target_encode), | | | ('passthrough', 'passthrough', features_to_passthrough) | | | ], | | | remainder='drop' # Ensures any unassigned columns are discarded | | | ) | | | # ========================================== | | | # 5. MODELING (LINEAR & XGBOOST VARIANTS) | |
| # ========================================== | |
| # --- 5.1 Standard Linear Regression --- | |
| # Model 1: No Transformation | |
| pipeline_lr = Pipeline([ | |
| ('preprocessor', preprocessor), | |
| ('model', LinearRegression()) | |
| ]) | |
| # Model 2: Transformed Target (log1p) | |
| pipeline_lr_log = Pipeline([ | |
| ('preprocessor', preprocessor), | |
| ('model', TransformedTargetRegressor( | |
| regressor=LinearRegression(), func=np.log1p, inverse_func=np.expm1 | |
| )) | |
| ]) | |
| # --- 5.2 Regularized Linear Models (Transformed Target) --- | |
| # Model 3: Ridge Regression (log1p) | |
| pipeline_ridge_log = Pipeline([ | |
| ('preprocessor', preprocessor), | |
| ('model', TransformedTargetRegressor( | |
| regressor=RidgeCV(alphas=np.logspace(-3, 3, 10)), func=np.log1p, inverse_func=np.expm1 | |
| )) | |
| ]) | |
| # Model 4: Lasso Regression (log1p) | |
| # Note: Increased max_iter for Lasso as it can struggle to converge on unscaled dummy variables | |
| pipeline_lasso_log = Pipeline([ | |
| ('preprocessor', preprocessor), | |
| ('model', TransformedTargetRegressor( | |
| regressor=LassoCV(alphas=np.logspace(-5, 2, 10), max_iter=10000), func=np.log1p, inverse_func=np.expm1 | |
| )) | |
| ]) | |
| # --- 5.3 XGBoost Models --- | |
| # Model 5: XGBoost (No Transformation) | |
| pipeline_xgb = Pipeline([ | |
| ('preprocessor', preprocessor), | |
| ('model', XGBRegressor(random_state=42)) | |
| ]) | |
| # Model 6: XGBoost (Transformed Target log1p) | |
| pipeline_xgb_log = Pipeline([ | |
| ('preprocessor', preprocessor), | |
| ('model', TransformedTargetRegressor( | |
| regressor=XGBRegressor(random_state=42), func=np.log1p, inverse_func=np.expm1 | |
| )) | |
| ]) | |
| # ========================================== | |
| # 6. HYPERPARAMETER TUNING & FITTING | | | # ========================================== | | | # 6.1 Fit all Linear Models directly (Ridge/Lasso handle tuning internally via CV) | | | pipeline_lr.fit(X_train, y_train) | | | pipeline_lr_log.fit(X_train, y_train) | | | pipeline_ridge_log.fit(X_train, y_train) | | | pipeline_lasso_log.fit(X_train, y_train) | | | # 6.2 Tune XGBoost (No Transformation) | | | # Notice the param grid keys start with 'model__' because the regressor is directly in the pipeline | |
| param_grid_standard = { | |
| 'model__n_estimators': [100, 250], | |
| 'model__learning_rate': [0.05, 0.1], | |
| 'model__max_depth': [4, 6] | |
| } | | | grid_search_xgb = GridSearchCV(pipeline_xgb, param_grid_standard, cv=5, scoring='neg_root_mean_squared_error', n_jobs=-1) | | | grid_search_xgb.fit(X_train, y_train) | | | best_xgb = grid_search_xgb.best_estimator_ | | | # 6.3 Tune XGBoost (Transformed Target) | | | # Notice the keys start with 'model__regressor__' because of the TransformedTargetRegressor wrapper | |
| param_grid_log = { | |
| 'model__regressor__n_estimators': [100, 250], | |
| 'model__regressor__learning_rate': [0.05, 0.1], | |
| 'model__regressor__max_depth': [4, 6] | |
| } | | | grid_search_xgb_log = GridSearchCV(pipeline_xgb_log, param_grid_log, cv=5, scoring='neg_root_mean_squared_error', n_jobs=-1) | | | grid_search_xgb_log.fit(X_train, y_train) | | | best_xgb_log = grid_search_xgb_log.best_estimator_ | | | # ========================================== | | | # 7. RESULTS ANALYSIS | |
| # ========================================== | |
| def evaluate_model(name, model, X_test, y_test): | |
| y_pred = model.predict(X_test) | |
| r2 = r2_score(y_test, y_pred) | |
| n, p = X_test.shape[0], X_test.shape[1] | |
| adj_r2 = 1 - (1 - r2) * (n - 1) / (n - p - 1) | |
| print(f"--- {name} Performance ---") | |
| print(f"RMSE: {root_mean_squared_error(y_test, y_pred):,.2f}") | |
| print(f"MAE: {mean_absolute_error(y_test, y_pred):,.2f}") | |
| print(f"MedAE: {median_absolute_error(y_test, y_pred):,.2f}") | |
| print(f"R2: {r2:.3f}") | |
| print(f"Adj R2: {adj_r2:.3f}\n") | |
| evaluate_model("Linear Regression (Standard)", pipeline_lr, X_test, y_test) | | | evaluate_model("Linear Regression (Log Target)", pipeline_lr_log, X_test, y_test) | | | evaluate_model("Ridge Regression (Log Target)", pipeline_ridge_log, X_test, y_test) | | | evaluate_model("Lasso Regression (Log Target)", pipeline_lasso_log, X_test, y_test) | | | evaluate_model("XGBoost (Standard, Tuned)", best_xgb, X_test, y_test) | | | evaluate_model("XGBoost (Log Target, Tuned)", best_xgb_log, X_test, y_test) | | | # ========================================== | | | # 8. FEATURE IMPORTANCE & EXPLAINABILITY | | | # ========================================== | | | # Extract feature names from ColumnTransformer | | | # Note: Adjust 'ohe' and 'categorical_features' based on the transformation-centric lists defined earlier | | | cat_features_out = preprocessor.named_transformers_['ohe'].named_steps['ohe'].get_feature_names_out(features_to_ohe).tolist() | | | all_feature_names = features_to_scale + cat_features_out + features_to_target_encode + features_to_passthrough | | | # --- 8.1 Linear Model Coefficients (Example using Ridge Log) --- | | | # Because Ridge is wrapped in TransformedTargetRegressor, we must access .regressor_ | | | ridge_estimator = pipeline_ridge_log.named_steps['model'].regressor_ | |
| coefs = pd.Series(ridge_estimator.coef_, index=all_feature_names).sort_values(key=abs, ascending=False) | |
| plt.figure(figsize=(10, 6)) | |
| sns.barplot(x=coefs.head(15).values, y=coefs.head(15).index, palette="coolwarm") | |
| plt.title("Top 15 Feature Coefficients (Ridge Regression - Log Target)") | |
| plt.xlabel("Coefficient Value (Log Scale Impact)") | |
| plt.show() | |
| # --- 8.2 Tree Model Permutation Importance (Example using XGBoost Log) --- | |
| # Transform test data to feed directly into the base regressor | | | X_test_transformed = preprocessor.transform(X_test) | | | # Since we are evaluating the base estimator of the wrapped model, we must transform the target manually | | | y_test_log = np.log1p(y_test) | | | xgb_estimator_log = best_xgb_log.named_steps['model'].regressor_ | | | # If we were evaluating the non-transformed XGBoost, it would be: | | | # xgb_estimator_standard = best_xgb.named_steps['model'] | | | # and we would pass y_test without transforming it. | | | result = permutation_importance(xgb_estimator_log, X_test_transformed, y_test_log, | |
| n_repeats=5, random_state=42, scoring='neg_root_mean_squared_error') | |
| perm_importance = pd.Series(result.importances_mean, index=all_feature_names).sort_values(ascending=False) | |
| plt.figure(figsize=(10, 6)) | |
| sns.barplot(x=perm_importance.head(15).values, y=perm_importance.head(15).index, palette="viridis") | |
| plt.title("Top 15 Feature Permutation Importances (XGBoost - Log Target)") | |
| plt.xlabel("Increase in RMSE when feature is shuffled") | | | plt.show() |