Capstone: Comparative Machine Learning Approaches for Predicting Scale Score Averages¶
In [ ]:
In [6]:
# Import packages
import pandas as pd
import seaborn as sns
import joblib
import os
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
import optuna
from optuna import visualization as opt_viz
from xgboost import XGBRegressor
from scipy.stats import pearsonr, ttest_rel, f_oneway, ttest_ind
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error, make_scorer
from statsmodels.formula.api import ols
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split, cross_val_score, KFold
from sklearn.model_selection import RandomizedSearchCV
# Suppress all warnings
import warnings
warnings.filterwarnings('ignore')
# magic word for producing visualizations in notebook
%matplotlib inline
# Ensure all columns are displayed
pd.set_option('display.max_columns', None) # Show all columns
pd.set_option('display.width', None) # Adjust display width to avoid wrapping
pd.set_option('display.max_colwidth', None) # Show full content of each column (if needed)
# Define a function to print bold headers
def print_bold(text):
print(f"\033[1m{text}\033[0m")
In [ ]:
Data¶
The datasets used in this study is taken from Delaware Open Data Council. Two datasets are used in this study, ‘student assessment performance’ dataset and ‘student discipline’ dataset.
Student Assessment Performance dataset (assessment.csv): https://data.delaware.gov/Education/Student-Assessment-Performance/ms6b-mt82/about_data Student Discipline (discipline.csv): https://data.delaware.gov/Education/Student-Discipline/yr4w-jdi4/about_data
In [ ]:
Research Question¶
How do traditional regression models compare to tree-based models in predicting school-level Scale Score Averages, based on predictive accuracy and model performance metrics?
In [ ]:
Hypothesis¶
Null Hypothesis (H₀): There is no significant difference in predictive accuracy between tree-based models (Decision Tree, XGBoost) and traditional regression models (Ridge Regression, Lasso Regression) in predicting Scale Score Averages.
Alternative Hypothesis (H₁): Tree-based models (Decision Tree, XGBoost) demonstrate significantly higher predictive accuracy than traditional regression models (Ridge Regression, Lasso Regression) in predicting Scale Score Averages.
In [ ]:
In [7]:
def split_data(filepath='4_ready_training.csv', test_size=0.2, random_state=42):
"""
Load the dataset, convert specific columns to categorical, dynamically select predictor columns,
and split the data into training and testing sets.
Parameters:
filepath (str): Path to the CSV file.
test_size (float): Proportion of the dataset to include in the test split.
random_state (int): Random seed for reproducibility.
Returns:
X_train (DataFrame): Training predictors.
X_test (DataFrame): Testing predictors.
y_train (Series): Training target variable.
y_test (Series): Testing target variable.
"""
# Load the dataset
df = pd.read_csv(filepath)
# Convert relevant columns to 'category' data type
df['School Year'] = df['School Year'].astype('category')
df['Grade_Numeric'] = df['Grade_Numeric'].astype('category')
# Dynamically find relevant columns
grade_numeric_cols = [col for col in df.columns if col.startswith('Grade_Numeric_')]
grade_group_cols = [col for col in df.columns if col.startswith('Grade_Group_')]
category_cols = [col for col in df.columns if col.lower().startswith('category_')]
content_area_cols = [col for col in df.columns if col.startswith('ContentArea_')]
special_demo_cols = [col for col in df.columns if col.startswith('SpecialDemo_')]
district_code_cols = [col for col in df.columns if col.startswith('DistrictCode_')]
school_year_cols = [col for col in df.columns if col.startswith('SchoolYear_')]
assessment_name_cols = [col for col in df.columns if col.startswith('AssessmentName_')]
# Define Dependent and Independent Variables
y = df['ScaleScoreAvg']
# Independent Variables (Predictors)
independent_vars = (
['Incidents', 'AvgDuration', 'PctEnrollment'] +
grade_numeric_cols + grade_group_cols + category_cols +
content_area_cols + special_demo_cols + district_code_cols +
school_year_cols + assessment_name_cols
)
X = df[independent_vars]
# --- Split Data into Training and Testing Sets ---
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=test_size, random_state=random_state
)
print("\n\nShape of X_train:", X_train.shape)
print("Shape of X_test:", X_test.shape)
print("Process successfully!\n")
return X_train, X_test, y_train, y_test
In [ ]:
Saved Tuned ML Models¶
In [3]:
## Traditional Regression
# Ridge_Regression_tuned
loaded_ridge_model_tuned = joblib.load("saved_models/ridge_model_tuned.joblib")
# lasso_regression_model_tuned
loaded_lasso_regression_tuned = joblib.load("saved_models/lasso_regression_model_tuned.joblib")
## Tree-based
# XGBoost_Regressor_Tuned
loaded_XGBoost_Regressor_Tuned_model = joblib.load("saved_models/XGBoost_Regressor_Tuned.joblib")
# decision_tree_model_tuned
loaded_decisiontree_model = joblib.load("saved_models/decision_tree_model_tuned.joblib")
In [ ]:
In [8]:
# Load Models
ridge = joblib.load("saved_models/ridge_model_tuned.joblib")
lasso = joblib.load("saved_models/lasso_regression_model_tuned.joblib")
dtree = joblib.load("saved_models/decision_tree_model_tuned.joblib")
xgb = joblib.load("saved_models/XGBoost_Regressor_Tuned.joblib")
# Define Models Dictionary
models = {
"Ridge Regression": ridge,
"Lasso Regression": lasso,
"Decision Tree": dtree,
"XGBoost": xgb
}
# Evaluation Results Storage
results = []
residuals_dict = {}
predictions_dict = {}
# Evaluation Loop
for name, model in models.items():
y_pred = model.predict(X_test)
r2 = r2_score(y_test, y_pred)
mse = mean_squared_error(y_test, y_pred)
mae = mean_absolute_error(y_test, y_pred)
rmse = np.sqrt(mse)
corr, p_val = pearsonr(y_test, y_pred)
residuals = y_test - y_pred
residuals_dict[name] = residuals
predictions_dict[name] = y_pred
residual_mean = np.mean(residuals)
residual_std = np.std(residuals)
cv_scores = cross_val_score(model, X_train, y_train, cv=3, scoring='r2')
cv_mean = cv_scores.mean()
results.append({
"Model": name,
"R² Score": round(r2, 6),
"MSE": round(mse, 6),
"RMSE": round(rmse, 6),
"MAE": round(mae, 6),
"Pearson Correlation": round(corr, 3),
"P-Value": f"{p_val:.3f}",
"Mean Residual": round(residual_mean, 6),
"Residual Std Dev": round(residual_std, 6),
"CV R² Scores": np.round(cv_scores, 6).tolist(),
"Mean CV R² Score": round(cv_mean, 6)
})
performance_table = pd.DataFrame(results)
# Display Performance Table
pd.set_option('display.max_colwidth', None)
display(performance_table)
# Residuals Scatter Plot for All Models
plt.figure(figsize=(12, 7))
for model_name, residuals in residuals_dict.items():
plt.scatter(y_test, residuals, alpha=0.5, label=f"Residuals {model_name}")
plt.axhline(0, color='black', linewidth=1, linestyle='--')
plt.title('Residual Analysis of Tuned Models', fontsize=16)
plt.xlabel('Actual ScaleScoreAvg', fontsize=14)
plt.ylabel('Residuals (Actual - Predicted)', fontsize=14)
plt.legend(title='Model', fontsize=10, title_fontsize=11)
plt.grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()
# Individual Residual Distribution per Model
for model_name, residuals in residuals_dict.items():
plt.figure(figsize=(10, 6))
sns.histplot(residuals, bins=50, kde=True, color='blue')
plt.axvline(0, color='black', linestyle='--', linewidth=1)
plt.title(f'Residual Distribution - {model_name}', fontsize=14)
plt.xlabel("Residual Error", fontsize=12)
plt.ylabel("Frequency", fontsize=12)
plt.grid(True, linestyle="--", alpha=0.6)
plt.tight_layout()
plt.show()
# === Trend Analysis by School Year ===
# Reconstruct 'School Year' from one-hot encoding
X_test_with_year = X_test.copy()
school_year_cols = [col for col in X_test_with_year.columns if 'SchoolYear_' in col]
X_test_with_year['School Year'] = X_test_with_year[school_year_cols].idxmax(axis=1).str.replace('SchoolYear_', '')
# Create trend_df for predictions
trend_df = pd.DataFrame({
'School Year': X_test_with_year['School Year'],
'Actual': y_test
})
for name, preds in predictions_dict.items():
trend_df[name] = preds
trend_df[f'Residuals_{name}'] = trend_df['Actual'] - preds
trend_analysis = trend_df.groupby('School Year').mean(numeric_only=True).reset_index()
print("\n\nTrend Analysis")
display(trend_analysis.head())
# Trend Plot: Predicted vs Actual
plt.figure(figsize=(12, 7))
plt.plot(trend_analysis['School Year'], trend_analysis['Actual'],
marker='o', linestyle='-', linewidth=2, color='black', label='Actual')
for name in models.keys():
plt.plot(trend_analysis['School Year'], trend_analysis[name],
marker='o', linestyle='--', linewidth=2, label=name)
plt.title('\n\nTrend Analysis: Predicted vs Actual by School Year', fontsize=16)
plt.xlabel('School Year', fontsize=14)
plt.ylabel('ScaleScoreAvg', fontsize=14)
plt.xticks(rotation=45)
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend(title='Model', fontsize=10, title_fontsize=11)
plt.tight_layout()
plt.show()
# Residual Plot: by School Year
plt.figure(figsize=(12, 7))
for name in models.keys():
plt.plot(trend_analysis['School Year'], trend_analysis[f'Residuals_{name}'],
marker='o', linestyle='--', label=f'Residuals {name}')
plt.axhline(0, color='black', linewidth=1, linestyle='--')
plt.title('\n\nResiduals by Model and School Year', fontsize=16)
plt.xlabel('School Year', fontsize=14)
plt.ylabel('Residuals (Actual - Predicted)', fontsize=14)
plt.xticks(rotation=45)
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend(title='Model', fontsize=10, title_fontsize=11)
plt.tight_layout()
plt.show()
# Bar Plot: Model Performance Metrics Comparison
metrics_bar = performance_table[['Model', 'MAE', 'MSE', 'RMSE', 'R² Score']].set_index('Model')
metrics_bar.plot(kind='bar', figsize=(12, 7))
plt.title('\n\nTuned Model Performance Metrics Comparison', fontsize=16)
plt.ylabel('Metric Value', fontsize=14)
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()
| Model | R² Score | MSE | RMSE | MAE | Pearson Correlation | P-Value | Mean Residual | Residual Std Dev | CV R² Scores | Mean CV R² Score | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Ridge Regression | 0.292510 | 0.279764 | 0.528927 | 0.401849 | 0.541 | 0.000 | 0.003585 | 0.528915 | [0.295127, 0.292894, 0.295866] | 0.294629 |
| 1 | Lasso Regression | 0.289127 | 0.281102 | 0.530191 | 0.403611 | 0.538 | 0.000 | 0.003544 | 0.530179 | [0.291559, 0.289416, 0.292087] | 0.291020 |
| 2 | Decision Tree | 0.852405 | 0.058364 | 0.241586 | 0.117435 | 0.923 | 0.000 | -0.000297 | 0.241586 | [0.84778, 0.849624, 0.848568] | 0.848657 |
| 3 | XGBoost | 0.864874 | 0.053433 | 0.231156 | 0.118371 | 0.930 | 0.000 | 0.000088 | 0.231156 | [0.8618, 0.861578, 0.862006] | 0.861795 |
Trend Analysis
| School Year | Actual | Ridge Regression | Residuals_Ridge Regression | Lasso Regression | Residuals_Lasso Regression | Decision Tree | Residuals_Decision Tree | XGBoost | Residuals_XGBoost | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2019 | 0.305992 | 0.307298 | -0.001306 | 0.300830 | 0.005162 | 0.299860 | 0.006132 | 0.305479 | 0.000512 |
| 1 | 2020 | 0.161507 | 0.155501 | 0.006006 | 0.155388 | 0.006119 | 0.162667 | -0.001160 | 0.161055 | 0.000451 |
| 2 | 2021 | 0.213464 | 0.208760 | 0.004704 | 0.219351 | -0.005886 | 0.225232 | -0.011768 | 0.218139 | -0.004674 |
| 3 | 2022 | 0.104016 | 0.098043 | 0.005973 | 0.101230 | 0.002786 | 0.102489 | 0.001527 | 0.101364 | 0.002653 |
| 4 | 2023 | 0.165420 | 0.162845 | 0.002575 | 0.159777 | 0.005643 | 0.167931 | -0.002511 | 0.167270 | -0.001849 |
In [ ]:
ANOVA Test (Individual). T-test (groups)¶
In [9]:
# Call function: split_data
X_train, X_test, y_train, y_test = split_data()
# Load Models
ridge = joblib.load("saved_models/ridge_model_tuned.joblib")
lasso = joblib.load("saved_models/lasso_regression_model_tuned.joblib")
dtree = joblib.load("saved_models/decision_tree_model_tuned.joblib")
xgb = joblib.load("saved_models/XGBoost_Regressor_Tuned.joblib")
# Set CV and Scorer
cv = 3
scorer = make_scorer(r2_score)
# Cross-Validated R² Scores
ridge_scores = cross_val_score(ridge, X_train, y_train, cv=cv, scoring=scorer)
lasso_scores = cross_val_score(lasso, X_train, y_train, cv=cv, scoring=scorer)
dtree_scores = cross_val_score(dtree, X_train, y_train, cv=cv, scoring=scorer)
xgb_scores = cross_val_score(xgb, X_train, y_train, cv=cv, scoring=scorer)
# Model Scores Dict
model_scores = {
"Ridge Regression": ridge_scores,
"Lasso Regression": lasso_scores,
"Decision Tree": dtree_scores,
"XGBoost": xgb_scores
}
# ANOVA Across All Models
f_stat, p_val_anova = f_oneway(*model_scores.values())
# Group Comparison: Regression vs Tree-Based
regression_group = np.concatenate([ridge_scores, lasso_scores])
tree_group = np.concatenate([dtree_scores, xgb_scores])
t_stat, p_val_group = ttest_ind(tree_group, regression_group, equal_var=False)
# Determine Winning Group
group_means = {
"Regression Models": regression_group.mean(),
"Tree-Based Models": tree_group.mean()
}
winning_group = max(group_means, key=group_means.get)
# Rank Models
model_means = {name: scores.mean() for name, scores in model_scores.items()}
ranked_models = sorted(model_means.items(), key=lambda x: x[1], reverse=True)
# Final Report
print("\nFinal Model Comparison Report\n")
print(f"ANOVA across all models: F = {f_stat:.4f}, p = {p_val_anova:.4f}")
if p_val_anova < 0.05:
print("Reject H₀: At least one model's predictive accuracy differs significantly from the others.")
else:
print("Fail to reject H₀: No significant difference in predictive accuracy across the models.")
print(f"\nT-test (Tree-Based vs Regression): T = {t_stat:.4f}, p = {p_val_group:.4f}")
if p_val_group < 0.05:
print("Reject H₀: Tree-based models have significantly different predictive accuracy than regression models.")
if np.mean(tree_group) > np.mean(regression_group):
print("Tree-based models outperform regression models (H₁ supported).")
else:
print("Regression models outperform tree-based models.")
else:
print("Fail to reject H₀: No significant difference in predictive accuracy between tree-based and regression models.")
print(f"Winning Model Group: {winning_group}")
print("\nModel Ranking (based on avg R² score):")
for i, (name, score) in enumerate(ranked_models, 1):
print(f"{i}. {name}: R² = {score:.4f}")
# Top 2 Models from the Winning Group
top_models = [name for name, score in ranked_models if ("Tree" in name or "XGBoost" in name)] if "Tree" in winning_group else \
[name for name, score in ranked_models if ("Ridge" in name or "Lasso" in name)]
top_2_models = top_models[:2]
print("\nTop 2 Models from Winning Group:", top_2_models)
# Residual Analysis for Grouped Histogram/Scatter Plots
# Predict on test set
y_pred_ridge = ridge.predict(X_test)
y_pred_lasso = lasso.predict(X_test)
y_pred_dtree = dtree.predict(X_test)
y_pred_xgb = xgb.predict(X_test)
# Residuals
regression_errors = {
"Ridge Regression": y_test - y_pred_ridge,
"Lasso Regression": y_test - y_pred_lasso
}
tree_errors = {
"Decision Tree": y_test - y_pred_dtree,
"XGBoost": y_test - y_pred_xgb
}
# Side-by-side Residual Error Histograms with KDE for Model Groups
# Create subplots for side-by-side visualization
fig, axes = plt.subplots(1, 2, figsize=(16, 6), sharey=True)
# Plot for Regression Models
for model, residuals in regression_errors.items():
sns.histplot(residuals, bins=50, kde=True, label=model, ax=axes[0], alpha=0.5)
axes[0].axvline(0, color='black', linestyle='--', linewidth=1)
axes[0].set_title("\n\nHistogram of Residual Errors: Regression Models", fontsize=13)
axes[0].set_xlabel("Residual Error")
axes[0].set_ylabel("Frequency")
axes[0].legend(title="Model")
axes[0].grid(True, linestyle="--", alpha=0.6)
# Plot for Tree-Based Models
for model, residuals in tree_errors.items():
sns.histplot(residuals, bins=50, kde=True, label=model, ax=axes[1], alpha=0.5)
axes[1].axvline(0, color='black', linestyle='--', linewidth=1)
axes[1].set_title("\n\nHistogram of Residual Errors: Tree-Based Models", fontsize=13)
axes[1].set_xlabel("Residual Error")
axes[1].legend(title="Model")
axes[1].grid(True, linestyle="--", alpha=0.6)
plt.tight_layout()
plt.show()
# Side-by-side Scatter Plots: Residuals vs Actual for Model Groups
# Create subplots for side-by-side visualization
fig, axes = plt.subplots(1, 2, figsize=(16, 6), sharey=True)
# Plot for Regression Models
for model, residuals in regression_errors.items():
axes[0].scatter(y_test, residuals, alpha=0.5, label=model)
axes[0].axhline(0, color='black', linestyle='--', linewidth=1)
axes[0].set_title("\n\nResiduals vs Actual: Regression Models", fontsize=13)
axes[0].set_xlabel("Actual ScaleScoreAvg")
axes[0].set_ylabel("Residuals (Actual - Predicted)")
axes[0].legend(title="Model")
axes[0].grid(True, linestyle="--", alpha=0.6)
# Plot for Tree-Based Models
for model, residuals in tree_errors.items():
axes[1].scatter(y_test, residuals, alpha=0.5, label=model)
axes[1].axhline(0, color='black', linestyle='--', linewidth=1)
axes[1].set_title("\n\nResiduals vs Actual: Tree-Based Models", fontsize=13)
axes[1].set_xlabel("Actual ScaleScoreAvg")
axes[1].legend(title="Model")
axes[1].grid(True, linestyle="--", alpha=0.6)
plt.tight_layout()
plt.show()
Shape of X_train: (519869, 84) Shape of X_test: (129968, 84) Process successfully! Final Model Comparison Report ANOVA across all models: F = 238910.4692, p = 0.0000 Reject H₀: At least one model's predictive accuracy differs significantly from the others. T-test (Tree-Based vs Regression): T = 181.1951, p = 0.0000 Reject H₀: Tree-based models have significantly different predictive accuracy than regression models. Tree-based models outperform regression models (H₁ supported). Winning Model Group: Tree-Based Models Model Ranking (based on avg R² score): 1. XGBoost: R² = 0.8618 2. Decision Tree: R² = 0.8487 3. Ridge Regression: R² = 0.2946 4. Lasso Regression: R² = 0.2910 Top 2 Models from Winning Group: ['XGBoost', 'Decision Tree']
In [ ]:
Feature Importance (Tree-Based Models: Decision Tree & XGBoost)¶
In [10]:
# Call function: split_data
X_train, X_test, y_train, y_test = split_data()
# Load tree-based models
dtree = joblib.load("saved_models/decision_tree_model_tuned.joblib")
xgb = joblib.load("saved_models/XGBoost_Regressor_Tuned.joblib")
# Extract feature importances
feature_importance_tree = {
"dtree": dtree.feature_importances_,
"xgb": xgb.feature_importances_,
}
# Create DataFrame
feature_importance_tree_df = pd.DataFrame(feature_importance_tree, index=X_train.columns)
# Compute mean importance and sort
feature_importance_tree_df["Mean_Importance"] = feature_importance_tree_df.mean(axis=1)
top_tree_features = feature_importance_tree_df.sort_values(by="Mean_Importance", ascending=False).head(15)
# Display results
print("\nFeature Importance (Tree-Based Models)")
print(top_tree_features.to_string())
# Plot
plt.figure(figsize=(10, 6))
top_tree_features["Mean_Importance"].plot(kind='barh', color='forestgreen', edgecolor='black')
plt.xlabel("Mean Feature Importance")
plt.ylabel("Feature")
plt.title("\n\nTop 15 Influential Features (Tree-Based Models)")
plt.gca().invert_yaxis()
plt.grid(axis="x", linestyle="--", alpha=0.7)
plt.tight_layout()
plt.show()
Shape of X_train: (519869, 84)
Shape of X_test: (129968, 84)
Process successfully!
Feature Importance (Tree-Based Models)
dtree xgb Mean_Importance
Incidents 0.385484 0.191493 0.288489
AvgDuration 0.312013 0.243002 0.277507
PctEnrollment 0.254819 0.042667 0.148743
AssessmentName_Smarter Balanced Summative Assessment 0.032142 0.107642 0.069892
AssessmentName_Delaware System of Student Assessment (DeSSA) 0.001464 0.039766 0.020615
AssessmentName_SAT School-Day (Spring) 0.001355 0.033250 0.017302
ContentArea_Mathematics 0.000102 0.031046 0.015574
category_Suspended 0.000910 0.021201 0.011056
ContentArea_Essay 0.003676 0.015484 0.009580
AssessmentName_DeSSA Alternate Assessment 0.000272 0.018355 0.009313
ContentArea_English Language Arts 0.000021 0.017583 0.008802
ContentArea_Computer Science 0.000286 0.014528 0.007407
category_Expulsion (Permanent Removal from School) 0.000035 0.012712 0.006374
SchoolYear_2023 0.000208 0.010592 0.005400
SchoolYear_2022 0.001944 0.008735 0.005340
In [ ]:
Feature Coefficients (Regression Models: Ridge & Lasso)¶
In [11]:
# Call function: split_data
X_train, X_test, y_train, y_test = split_data()
# Load regression models
lasso = joblib.load("saved_models/lasso_regression_model_tuned.joblib")
ridge = joblib.load("saved_models/ridge_model_tuned.joblib")
# Extract absolute coefficients
feature_coefficients = {
"ridge": np.abs(ridge.coef_),
"lasso": np.abs(lasso.coef_),
}
# Create DataFrame
feature_coefficients_df = pd.DataFrame(feature_coefficients, index=X_train.columns)
# Compute mean coefficient magnitude and sort
feature_coefficients_df["Mean_Coefficient"] = feature_coefficients_df.mean(axis=1)
top_coeff_features = feature_coefficients_df.sort_values(by="Mean_Coefficient", ascending=False).head(15)
# Display results
print("\nFeature Coefficients (Regression Models)")
print(top_coeff_features.to_string())
# Plot
plt.figure(figsize=(10, 6))
top_coeff_features["Mean_Coefficient"].plot(kind='barh', color='steelblue', edgecolor='black')
plt.xlabel("Mean Absolute Coefficient")
plt.ylabel("Feature")
plt.title("\n\nTop 15 Influential Features (Regression Models)")
plt.gca().invert_yaxis()
plt.grid(axis="x", linestyle="--", alpha=0.7)
plt.tight_layout()
plt.show()
Shape of X_train: (519869, 84)
Shape of X_test: (129968, 84)
Process successfully!
Feature Coefficients (Regression Models)
ridge lasso Mean_Coefficient
Incidents 1.682993 1.511145 1.597069
AvgDuration 0.996360 0.000000 0.498180
PctEnrollment 0.283685 0.281220 0.282453
AssessmentName_Smarter Balanced Summative Assessment 0.247759 0.236320 0.242040
AssessmentName_SAT School-Day (Spring) 0.155035 0.163387 0.159211
DistrictCode_68 0.318391 0.000000 0.159195
AssessmentName_Delaware System of Student Assessment (DeSSA) 0.144578 0.123425 0.134002
DistrictCode_9614 0.242363 0.000000 0.121182
ContentArea_Essay 0.077516 0.076836 0.077176
DistrictCode_9606 0.132770 0.000000 0.066385
DistrictCode_9605 0.119697 0.000000 0.059848
DistrictCode_15 0.082731 0.027224 0.054978
category_Expulsion (Permanent Removal from School) 0.056744 0.035242 0.045993
DistrictCode_88 0.087946 0.000000 0.043973
DistrictCode_82 0.087017 0.000000 0.043509
In [ ]:
Winning group: Tree-based - decision tree and XGBoost - t-test¶
In [12]:
# Load models
dtree = joblib.load("saved_models/decision_tree_model_tuned.joblib")
xgb = joblib.load("saved_models/XGBoost_Regressor_Tuned.joblib")
# Compute Cross-Validated R² Scores
cv = 3 # 3-Fold Cross Validation
scorer = make_scorer(r2_score)
dtree_scores = cross_val_score(dtree, X_train, y_train, cv=cv, scoring=scorer)
xgb_scores = cross_val_score(xgb, X_train, y_train, cv=cv, scoring=scorer)
# Print T-Test Results
stat, p_value = ttest_ind(dtree_scores, xgb_scores, equal_var=False)
mean_dtree = np.mean(dtree_scores)
mean_xgb = np.mean(xgb_scores)
mean_diff = mean_xgb - mean_dtree
# Print T-Test Results
print("\n=== T-Test: Decision Tree vs. XGBoost ===")
print(f"Decision Tree Mean R²: {mean_dtree:.4f}")
print(f"XGBoost Mean R²: {mean_xgb:.4f}")
print(f"Mean Difference: {mean_diff:.4f}")
print(f"T-Statistic: {stat:.4f}")
print(f"P-Value: {p_value:.4f}")
if p_value < 0.05:
print("XGBoost significantly outperforms Decision Tree (p < 0.05).")
else:
print("No significant difference between XGBoost and Decision Tree.")
=== T-Test: Decision Tree vs. XGBoost === Decision Tree Mean R²: 0.8487 XGBoost Mean R²: 0.8618 Mean Difference: 0.0131 T-Statistic: -23.9538 P-Value: 0.0010 XGBoost significantly outperforms Decision Tree (p < 0.05).
In [ ]:
Winning group: Tree-based - decision tree and XGBoost - ANOVA¶
In [13]:
# Load Models
dtree = joblib.load("saved_models/decision_tree_model_tuned.joblib")
xgb = joblib.load("saved_models/XGBoost_Regressor_Tuned.joblib")
# Compute R² Scores Using Cross-Validation
cv = 3
scorer = make_scorer(r2_score)
dtree_scores = cross_val_score(dtree, X_train, y_train, cv=cv, scoring=scorer)
xgb_scores = cross_val_score(xgb, X_train, y_train, cv=cv, scoring=scorer)
# Create DataFrame for ANOVA
df_anova = pd.DataFrame({
'R2_Score': np.concatenate([dtree_scores, xgb_scores]),
'Model': ['Decision Tree'] * len(dtree_scores) + ['XGBoost'] * len(xgb_scores)
})
# Perform One-Way ANOVA Test
model = ols('R2_Score ~ C(Model)', data=df_anova).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
# Output Results
print("\n=== One-Way ANOVA: Decision Tree vs. XGBoost ===")
print(anova_table)
p_value = anova_table['PR(>F)'][0]
if p_value < 0.05:
print("XGBoost significantly outperforms Decision Tree (p < 0.05).")
else:
print("No significant difference between XGBoost and Decision Tree.")
=== One-Way ANOVA: Decision Tree vs. XGBoost ===
sum_sq df F PR(>F)
C(Model) 0.000259 1.0 573.78397 0.000018
Residual 0.000002 4.0 NaN NaN
XGBoost significantly outperforms Decision Tree (p < 0.05).
In [ ]:
In [5]:
# Call function: split_data
X_train, X_test, y_train, y_test = split_data()
# Load Tuned Top 2 Tree-Based Models
dtree = joblib.load("saved_models/decision_tree_model_tuned.joblib")
xgb = joblib.load("saved_models/XGBoost_Regressor_Tuned.joblib")
# Generate Predictions
y_pred_dtree = dtree.predict(X_test)
y_pred_xgb = xgb.predict(X_test)
#: Build Residuals DataFrame
residuals_df = pd.DataFrame({
'Actual': y_test,
'Decision Tree': y_pred_dtree,
'XGBoost': y_pred_xgb,
'Residual_DT': y_test - y_pred_dtree,
'Residual_XGB': y_test - y_pred_xgb
})
# Add School Year column back (if encoded)
if 'School Year' not in residuals_df.columns and any('SchoolYear_' in col for col in X_test.columns):
school_year_cols = [col for col in X_test.columns if 'SchoolYear_' in col]
residuals_df['School Year'] = (
X_test[school_year_cols].idxmax(axis=1).str.replace("SchoolYear_", "")
)
else:
residuals_df['School Year'] = "Unknown"
# Extract residuals for plotting
tree_errors = {
'Decision Tree': residuals_df['Residual_DT'],
'XGBoost': residuals_df['Residual_XGB']
}
# Residuals vs Actual + Actual vs Predicted (Side-by-Side)
fig, axes = plt.subplots(1, 2, figsize=(18, 6))
# Residuals vs Actual
axes[0].scatter(residuals_df['Actual'], residuals_df['Residual_DT'], alpha=0.5, label='Decision Tree')
axes[0].scatter(residuals_df['Actual'], residuals_df['Residual_XGB'], alpha=0.5, label='XGBoost')
axes[0].axhline(0, color='black', linestyle='--')
axes[0].set_title("Residuals vs Actual", fontsize=13)
axes[0].set_xlabel("Actual ScaleScoreAvg")
axes[0].set_ylabel("Residual (Actual - Predicted)")
axes[0].legend()
axes[0].grid(True, linestyle='--', alpha=0.6)
# Actual vs Predicted
sns.scatterplot(data=residuals_df, x='Actual', y='Decision Tree', label="Decision Tree", alpha=0.6, ax=axes[1])
sns.scatterplot(data=residuals_df, x='Actual', y='XGBoost', label="XGBoost", alpha=0.6, ax=axes[1])
axes[1].plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--')
axes[1].set_title("Actual vs Predicted", fontsize=13)
axes[1].set_xlabel("Actual")
axes[1].set_ylabel("Predicted")
axes[1].legend()
axes[1].grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()
# Histograms of Residual Errors (Side-by-Side)
fig, axes = plt.subplots(1, len(tree_errors), figsize=(6 * len(tree_errors), 5), sharey=True)
if len(tree_errors) == 1:
axes = [axes]
for ax, (model, residuals) in zip(axes, tree_errors.items()):
sns.histplot(residuals, bins=50, kde=True, ax=ax, color='skyblue')
ax.axvline(0, color='black', linestyle='--')
ax.set_title(f"{model} Residuals", fontsize=12)
ax.set_xlabel("Residual Error")
ax.set_ylabel("Frequency")
ax.grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()
# Residuals vs Actual (Model-Specific)
fig, axs = plt.subplots(1, 2, figsize=(14, 6), sharey=True)
# Decision Tree
sns.scatterplot(x=y_test, y=tree_errors["Decision Tree"], ax=axs[0],
alpha=0.4, color='blue', edgecolor='k', linewidth=0.2)
axs[0].axhline(0, color='black', linestyle='--')
axs[0].set_title("Residuals vs Actual: Decision Tree")
axs[0].set_xlabel("Actual ScaleScoreAvg")
axs[0].set_ylabel("Residual (Actual - Predicted)")
axs[0].grid(True, linestyle='--', alpha=0.6)
# XGBoost
sns.scatterplot(x=y_test, y=tree_errors["XGBoost"], ax=axs[1],
alpha=0.4, color='orange', edgecolor='k', linewidth=0.2)
axs[1].axhline(0, color='black', linestyle='--')
axs[1].set_title("Residuals vs Actual: XGBoost")
axs[1].set_xlabel("Actual ScaleScoreAvg")
axs[1].grid(True, linestyle='--', alpha=0.6)
plt.suptitle("Residual Analysis: Tree-Based Models", fontsize=15)
plt.tight_layout()
plt.show()
# School Year Trend Analysis (if available)
if 'School Year' in residuals_df.columns and residuals_df['School Year'].nunique() > 1:
trend_summary = residuals_df.groupby('School Year').mean(numeric_only=True).reset_index()
print("\n=== Data Used for Trend Graphs ===")
display(trend_summary[['School Year', 'Actual', 'Decision Tree', 'XGBoost', 'Residual_DT', 'Residual_XGB']])
fig, axes = plt.subplots(1, 2, figsize=(18, 6), sharex=True)
# Performance Over Time
axes[0].plot(trend_summary['School Year'], trend_summary['Actual'], marker='o', label='Actual', color='black')
axes[0].plot(trend_summary['School Year'], trend_summary['Decision Tree'], marker='o', linestyle='--', label='Decision Tree')
axes[0].plot(trend_summary['School Year'], trend_summary['XGBoost'], marker='o', linestyle='--', label='XGBoost')
axes[0].set_title("Predicted Performance by School Year", fontsize=13)
axes[0].set_xlabel("School Year")
axes[0].set_ylabel("ScaleScoreAvg")
axes[0].legend()
axes[0].tick_params(axis='x', rotation=45)
axes[0].grid(True, linestyle='--', alpha=0.6)
# Residuals Over Time
axes[1].plot(trend_summary['School Year'], trend_summary['Residual_DT'], marker='o', linestyle='--', label='Residuals: Decision Tree')
axes[1].plot(trend_summary['School Year'], trend_summary['Residual_XGB'], marker='o', linestyle='--', label='Residuals: XGBoost')
axes[1].axhline(0, color='black', linestyle='--')
axes[1].set_title("Residual Trends by School Year", fontsize=13)
axes[1].set_xlabel("School Year")
axes[1].set_ylabel("Residual (Actual - Predicted)")
axes[1].legend()
axes[1].tick_params(axis='x', rotation=45)
axes[1].grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()
Shape of X_train: (519869, 84) Shape of X_test: (129968, 84) Process successfully!
=== Data Used for Trend Graphs ===
| School Year | Actual | Decision Tree | XGBoost | Residual_DT | Residual_XGB | |
|---|---|---|---|---|---|---|
| 0 | 2019 | 0.305992 | 0.299860 | 0.305479 | 0.006132 | 0.000512 |
| 1 | 2020 | 0.161507 | 0.162667 | 0.161055 | -0.001160 | 0.000451 |
| 2 | 2021 | 0.213464 | 0.225232 | 0.218139 | -0.011768 | -0.004674 |
| 3 | 2022 | 0.104016 | 0.102489 | 0.101364 | 0.001527 | 0.002653 |
| 4 | 2023 | 0.165420 | 0.167931 | 0.167270 | -0.002511 | -0.001849 |
In [ ]:
In [14]:
# Load the saved XGBoost Tuned model
model_filename = "saved_models/XGBoost_Regressor_Tuned.joblib"
loaded_model = joblib.load(model_filename)
print(f"Model loaded: {model_filename}")
# Load the preprocessed and encoded base dataset
base_data_file = '4_ready_training-Copy1.csv' # Fully processed and encoded dataset
base_df = pd.read_csv(base_data_file)
print(f"Base data loaded: {base_data_file}")
# Proceed with the full dataset
print(f"Using entire dataset: {base_df.shape[0]} rows")
original_columns = base_df.columns.tolist()
identifier_columns = ['School Code', 'District Code', 'Assessment Name', 'ContentArea', 'Grade']
# Ensure all identifier columns exist; if missing, add with default 'Unknown'
for col in identifier_columns:
if col not in base_df.columns:
print(f"Warning: {col} not found in base data. Adding default value.")
base_df[col] = 'Unknown'
# Fit Target Scaler on Real Training Data
target_scaler = StandardScaler()
real_target = base_df['ScaleScoreAvg'].values.reshape(-1, 1)
target_scaler.fit(real_target)
print("Target scaler fitted on real training target.")
# Bootstrap Sampling for Simulation
num_simulated_rows = 20000
# Bootstrap sample from full dataset
simulated_data = base_df.sample(n=num_simulated_rows, replace=True, random_state=42).reset_index(drop=True)
# Ensure all original columns are present; add missing ones with defaults
for col in original_columns:
if col not in simulated_data.columns:
simulated_data[col] = 'Unknown' if col in identifier_columns else 0
# Remove duplicate columns
simulated_data = simulated_data.loc[:, ~simulated_data.columns.duplicated()]
# Save simulated data for reuse
simulated_data_file = 'test_xgboost_model_tuned.csv'
simulated_data.to_csv(simulated_data_file, index=False)
print(f"Simulated data saved to '{simulated_data_file}'")
# Hyperparameter Tuning
expected_features = pd.Index(loaded_model.feature_names_in_).unique()
simulated_data_encoded = simulated_data.reindex(columns=expected_features, fill_value=0)
X = simulated_data_encoded
y = simulated_data['ScaleScoreAvg']
# define the parameter grid for Randomized Search
param_dist = {
'n_estimators': [100, 150, 200, 250, 300, 400, 500],
'max_depth': [3, 4, 5, 6, 7, 8, 10],
'learning_rate': [0.01, 0.03, 0.05, 0.07, 0.1],
'subsample': [0.6, 0.7, 0.8, 0.9, 1.0],
'colsample_bytree': [0.6, 0.7, 0.8, 0.9, 1.0],
'reg_alpha': [1e-8, 1e-4, 1e-2, 0.1, 0.5, 1.0],
'reg_lambda': [1e-8, 1e-4, 1e-2, 0.1, 0.5, 1.0],
'random_state': [42]
}
# Set up cross-validation and model
xgb = XGBRegressor()
kf = KFold(n_splits=3, shuffle=True, random_state=42)
# Perform Randomized Search
random_search = RandomizedSearchCV(
estimator=xgb,
param_distributions=param_dist,
n_iter=10,
scoring='r2',
cv=kf,
verbose=1,
random_state=42,
n_jobs=-1
)
# Fit RandomizedSearchCV
random_search.fit(X, y)
# Output best parameters
print("\nBest Hyperparameters (RandomizedSearchCV):")
print(random_search.best_params_)
# Train model with best parameters
model = random_search.best_estimator_
model.fit(X, y)
# Make Predictions and Inverse Transform
simulated_data_encoded = simulated_data.reindex(columns=expected_features, fill_value=0)
predicted_scaled = model.predict(simulated_data_encoded)
predicted_scalescoreavg = target_scaler.inverse_transform(predicted_scaled.reshape(-1, 1)).flatten()
# Add noise to simulate more realistic target values
noise_multiplier = 0.05
noise = np.random.normal(0, noise_multiplier * np.std(predicted_scalescoreavg), size=predicted_scalescoreavg.shape)
simulated_scalescoreavg = np.clip(predicted_scalescoreavg + noise, 0, 100)
# Combine Predictions with Original Columns
predictions_df = simulated_data.copy()
predictions_df['Predicted_ScaleScoreAvg'] = predicted_scalescoreavg
predictions_df['ScaleScoreAvg'] = simulated_scalescoreavg
# Evaluate Model Performance
y_actual = predictions_df['ScaleScoreAvg']
y_pred = predictions_df['Predicted_ScaleScoreAvg']
r2 = r2_score(y_actual, y_pred)
mse = mean_squared_error(y_actual, y_pred)
mae = mean_absolute_error(y_actual, y_pred)
if np.var(y_pred) > 0 and np.var(y_actual) > 0:
correlation, p_value = pearsonr(y_actual, y_pred)
else:
correlation, p_value = float('nan'), float('nan')
print("\n=== Prediction Performance on Simulated Full Dataset ===")
print(f"R² Score: {r2:.3f}")
print(f"MSE: {mse:.6f}")
print(f"MAE: {mae:.6f}")
print(f"Pearson Correlation: {correlation:.3f}, P-value: {p_value:.3f}")
# Visualize Actual vs Predicted Distributions
plt.figure(figsize=(8, 6))
sns.histplot(y_actual, color='blue', label='Actual', kde=True)
sns.histplot(y_pred, color='red', label='Predicted', kde=True)
plt.title("Distribution of Actual vs Predicted Values (Simulated Data)")
plt.xlabel("ScaleScoreAvg")
plt.ylabel("Frequency")
plt.legend()
plt.tight_layout()
plt.show()
Model loaded: saved_models/XGBoost_Regressor_Tuned.joblib
Base data loaded: 4_ready_training-Copy1.csv
Using entire dataset: 649837 rows
Target scaler fitted on real training target.
Simulated data saved to 'test_xgboost_model_tuned.csv'
Fitting 3 folds for each of 10 candidates, totalling 30 fits
Best Hyperparameters (RandomizedSearchCV):
{'subsample': 0.9, 'reg_lambda': 0.1, 'reg_alpha': 0.1, 'random_state': 42, 'n_estimators': 250, 'max_depth': 6, 'learning_rate': 0.05, 'colsample_bytree': 0.8}
=== Prediction Performance on Simulated Full Dataset ===
R² Score: 0.854
MSE: 0.013996
MAE: 0.048477
Pearson Correlation: 0.955, P-value: 0.000
In [ ]:
In [ ]: