Building Actuarial Analytics with Snowflake and Anaconda

Now that we’ve cleaned our data and engineered our features, we’ve visualized patterns. We can actually predict health insurance coverage from demographics.

from sklearn.preprocessing import LabelEncoder, StandardScaler
from sklearn.model_selection import train_test_split

categorical_features = ['education', 'marital_status', 'race', 'income']
numerical_features = ['age', 'income_log', 'family_size', 'employed',
                     'income_per_capita', 'age_squared', 'is_senior', 'is_young_adult']

# Encode categorical variables
df_encoded = df_model.copy()
label_encoders = {}

for col in categorical_features:
    le = LabelEncoder()
    df_encoded[col + '_encoded'] = le.fit_transform(df_encoded[col])
    label_encoders[col] = le

# Create feature matrix
feature_cols = numerical_features + [col + '_encoded' for col in categorical_features]
X = df_encoded[feature_cols]
y = df_encoded['hicoverage']

# Encode target variable
le_target = LabelEncoder()
y_encoded = le_target.fit_transform(y)

A note for how we’re handling income_log and income, two representations of income as both a categorical feature and a numerical feature, this may feel like feature duplication, but it actually highlights how different representations can capture different patterns. The continuous relationship, captured in numerical income, is a smooth gradual relationship. As income increases the probability of coverage increases gradually. This is what income_log (the continuous numerical feature) captures. Alternatively, income in categorical features will show patterns in how discrete thresholds may trigger changes in coverage eligibility. Different models benefit from different representations.

categorical_features = ['education', 'marital_status', 'race', 'income']
numerical_features = ['age', 'income_log', 'family_size', 'employed',
                     'income_per_capita', 'age_squared', 'is_senior', 'is_young_adult']

# Encode categorical variables
df_encoded = df_model.copy()
label_encoders = {}

for col in categorical_features:
    le = LabelEncoder()
    df_encoded[col + '_encoded'] = le.fit_transform(df_encoded[col])
    label_encoders[col] = le

# Create feature matrix
feature_cols = numerical_features + [col + '_encoded' for col in categorical_features]
X = df_encoded[feature_cols]
y = df_encoded['hicoverage']

# Encode target variable
le_target = LabelEncoder()
y_encoded = le_target.fit_transform(y)

# Split data
X_train, X_test, y_train, y_test = train_test_split(
    X, y_encoded, test_size=0.2, random_state=42, stratify=y_encoded
)

print(f"Training set: {len(X_train)} samples")
print(f"Test set: {len(X_test)} samples")
print(f"Number of features: {X_train.shape[1]}")
print(f"Coverage types: {le_target.classes_}")

# Scale features
scaler = StandardScaler()
#X_train_scaled = scaler.fit_transform(X_train)
#X_test_scaled = scaler.transform(X_test)

# Initialize models
models = {
    'Logistic Regression': LogisticRegression(max_iter=1000, random_state=42, multi_class='multinomial'),
    'Random Forest': RandomForestClassifier(n_estimators=100, max_depth=10, random_state=42),
    'Gradient Boosting': GradientBoostingClassifier(n_estimators=100, max_depth=5, random_state=42)
}

# Train and evaluate models
results = {}

for model_name, model in models.items():
    print(f"\n--- Training {model_name} ---")
    
    # Use scaled data for Logistic Regression, original for tree-based
    if 'Logistic' in model_name:
        model.fit(X_train_scaled, y_train)
        y_pred = model.predict(X_test_scaled)
        y_pred_proba = model.predict_proba(X_test_scaled)
    else:
        model.fit(X_train, y_train)
        y_pred = model.predict(X_test)
        y_pred_proba = model.predict_proba(X_test)
    
    # Calculate metrics
    accuracy = accuracy_score(y_test, y_pred)
    
    print(f"Accuracy: {accuracy:.4f}")
    print("\nClassification Report:")
    print(classification_report(y_test, y_pred, target_names=le_target.classes_))
    
    results[model_name] = {
        'model': model,
        'accuracy': accuracy,
        'y_pred': y_pred,
        'y_pred_proba': y_pred_proba,
        'y_test': y_test
    }

# Plot model comparison
plot_model_comparison(results, le_target)

# Return best model and encoders
best_model_name = max(results, key=lambda x: results[x]['accuracy'])
print(f"\n*** Best Model: {best_model_name} (Accuracy: {results[best_model_name]['accuracy']:.4f}) ***")

return results[best_model_name]['model'], le_target, label_encoders, scaler, feature_cols

We’re labeling our features, encoding them as numbers, and ensuring that we’re applying the same encoding as new data arrives. Next, we create the feature mix. We create a matrix for every variable we use to make predictions.

X = df_encoded[feature_cols]
y = df_encoded['hicoverage']

X is the design matrix for predictors and y is the response vector for outcome. Each row is a person, each column is a feature, and y tells us their actual coverage type. Then, we encode the target so we can translate the result.

Next, we split our data, following the golden rule: Never test on data you trained on. We set random_state=42 to help with debugging. If accuracy changes, it’s due to code changes and not random split variation—with a special nod to the answer to the universe: 42—and we stratify y_encoded to make sure that each parameter will match distributions between our train and test.

We scale for the distance-based algorithms: logistic, neural networks, support vector machine (SVM). In tree-based models, like random forest and gradient boosting, we don’t care about feature scaling. The splits work the same regardless of the distance between calculations. We create an if / else block to only scale for the logistic models.