Types of Basketball Injuries and Risk Factors

Common NBA Injuries

Primary Risk Factors

Load Management and Tracking Data

Wearable Technology and Tracking Systems

Key Load Monitoring Metrics

Load Management Strategies

Python: Machine Learning Injury Prediction Model

Feature Engineering and Predictive Modeling

This example demonstrates building a gradient boosting classifier to predict injury risk using player tracking data and workload metrics.

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier
from sklearn.metrics import classification_report, roc_auc_score, roc_curve
from sklearn.metrics import confusion_matrix, precision_recall_curve
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime, timedelta

# Load player tracking and injury data
def load_and_prepare_data():
    """
    Load player data with tracking metrics, workload, and injury outcomes
    """
    # Example data structure
    data = pd.read_csv('player_tracking_data.csv')

    # Features expected in dataset:
    # - player_id, date, age, position
    # - minutes_played, distance_total, high_speed_distance
    # - player_load, jump_count, decel_events, accel_events
    # - avg_speed, max_speed
    # - days_since_injury, previous_injury_count
    # - back_to_back (binary), travel_hours
    # - sleep_hours, hrv_score, wellness_score
    # - injury_next_7days (target: 0=no injury, 1=injury)

    return data

def engineer_features(df):
    """
    Create advanced features for injury prediction
    """
    df = df.sort_values(['player_id', 'date'])

    # Calculate rolling workload metrics
    for days in [7, 14, 28]:
        df[f'load_{days}d'] = df.groupby('player_id')['player_load'].transform(
            lambda x: x.rolling(days, min_periods=1).mean()
        )
        df[f'minutes_{days}d'] = df.groupby('player_id')['minutes_played'].transform(
            lambda x: x.rolling(days, min_periods=1).sum()
        )

    # Acute:Chronic Workload Ratio (ACWR)
    df['acwr'] = df['load_7d'] / df['load_28d']
    df['acwr'] = df['acwr'].fillna(1.0)

    # Workload changes (week-to-week)
    df['load_change_pct'] = df.groupby('player_id')['player_load'].pct_change(periods=7)

    # Cumulative load monotony (variation coefficient)
    df['load_monotony'] = df.groupby('player_id')['player_load'].transform(
        lambda x: x.rolling(7, min_periods=1).mean() / (x.rolling(7, min_periods=1).std() + 0.1)
    )

    # High-intensity work ratio
    df['high_intensity_ratio'] = df['high_speed_distance'] / (df['distance_total'] + 0.1)

    # Exposure time features
    df['minutes_cumulative_14d'] = df['minutes_14d']
    df['games_played_7d'] = df.groupby('player_id')['minutes_played'].transform(
        lambda x: (x.rolling(7, min_periods=1).count())
    )

    # Recovery markers
    df['recovery_score'] = (df['sleep_hours'] / 8.0) * (df['hrv_score'] / 100.0) * (df['wellness_score'] / 10.0)

    # Days since last high-load game
    high_load_threshold = df['player_load'].quantile(0.75)
    df['high_load_game'] = (df['player_load'] > high_load_threshold).astype(int)
    df['days_since_high_load'] = df.groupby('player_id').apply(
        lambda x: (x['date'] - x[x['high_load_game'] == 1]['date'].shift()).dt.days
    ).reset_index(level=0, drop=True)
    df['days_since_high_load'] = df['days_since_high_load'].fillna(99)

    # Age-related risk
    df['age_risk_score'] = np.where(df['age'] > 30, (df['age'] - 30) * 0.5, 0)

    # Injury history interaction
    df['history_load_interaction'] = df['previous_injury_count'] * df['acwr']

    return df

def create_risk_zones(acwr):
    """
    Categorize ACWR into risk zones
    """
    if acwr < 0.8:
        return 'detraining'
    elif 0.8 <= acwr <= 1.3:
        return 'optimal'
    elif 1.3 < acwr <= 1.5:
        return 'moderate_risk'
    else:
        return 'high_risk'

def build_injury_prediction_model(df):
    """
    Train gradient boosting model for injury prediction
    """
    # Select features
    feature_cols = [
        'age', 'minutes_played', 'distance_total', 'high_speed_distance',
        'player_load', 'jump_count', 'decel_events', 'accel_events',
        'load_7d', 'load_14d', 'load_28d', 'minutes_7d', 'minutes_28d',
        'acwr', 'load_change_pct', 'load_monotony', 'high_intensity_ratio',
        'games_played_7d', 'recovery_score', 'days_since_high_load',
        'days_since_injury', 'previous_injury_count', 'age_risk_score',
        'history_load_interaction', 'back_to_back', 'travel_hours',
        'sleep_hours', 'hrv_score', 'wellness_score'
    ]

    # Remove rows with missing target or excessive missing features
    df_model = df.dropna(subset=['injury_next_7days'])
    df_model = df_model.dropna(subset=feature_cols, thresh=len(feature_cols)-3)
    df_model[feature_cols] = df_model[feature_cols].fillna(df_model[feature_cols].median())

    X = df_model[feature_cols]
    y = df_model['injury_next_7days']

    # Split data temporally (train on earlier dates, test on later)
    split_date = df_model['date'].quantile(0.75)
    train_mask = df_model['date'] < split_date

    X_train, X_test = X[train_mask], X[~train_mask]
    y_train, y_test = y[train_mask], y[~train_mask]

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

    # Train Gradient Boosting Classifier
    # Note: Injury data is typically highly imbalanced (few injuries)
    injury_rate = y_train.mean()
    scale_pos_weight = (1 - injury_rate) / injury_rate

    gb_model = GradientBoostingClassifier(
        n_estimators=200,
        learning_rate=0.05,
        max_depth=5,
        min_samples_split=20,
        min_samples_leaf=10,
        subsample=0.8,
        random_state=42
    )

    gb_model.fit(X_train_scaled, y_train)

    # Predictions
    y_pred = gb_model.predict(X_test_scaled)
    y_pred_proba = gb_model.predict_proba(X_test_scaled)[:, 1]

    # Evaluation
    print("Gradient Boosting Model Performance")
    print("=" * 50)
    print(classification_report(y_test, y_pred, target_names=['No Injury', 'Injury']))
    print(f"\nROC-AUC Score: {roc_auc_score(y_test, y_pred_proba):.3f}")

    # Feature importance
    feature_importance = pd.DataFrame({
        'feature': feature_cols,
        'importance': gb_model.feature_importances_
    }).sort_values('importance', ascending=False)

    print("\nTop 10 Most Important Features:")
    print(feature_importance.head(10))

    return gb_model, scaler, feature_cols, y_test, y_pred_proba

def plot_model_performance(y_test, y_pred_proba):
    """
    Visualize model performance metrics
    """
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # ROC Curve
    fpr, tpr, thresholds = roc_curve(y_test, y_pred_proba)
    auc = roc_auc_score(y_test, y_pred_proba)

    axes[0].plot(fpr, tpr, label=f'ROC Curve (AUC = {auc:.3f})', linewidth=2)
    axes[0].plot([0, 1], [0, 1], 'k--', label='Random Classifier')
    axes[0].set_xlabel('False Positive Rate')
    axes[0].set_ylabel('True Positive Rate')
    axes[0].set_title('ROC Curve - Injury Prediction')
    axes[0].legend()
    axes[0].grid(alpha=0.3)

    # Precision-Recall Curve
    precision, recall, pr_thresholds = precision_recall_curve(y_test, y_pred_proba)

    axes[1].plot(recall, precision, linewidth=2)
    axes[1].set_xlabel('Recall')
    axes[1].set_ylabel('Precision')
    axes[1].set_title('Precision-Recall Curve')
    axes[1].grid(alpha=0.3)

    plt.tight_layout()
    plt.savefig('injury_model_performance.png', dpi=300, bbox_inches='tight')
    plt.show()

def calculate_risk_score(player_data, model, scaler, feature_cols):
    """
    Calculate injury risk score for a player
    """
    # Prepare features
    X = player_data[feature_cols].values.reshape(1, -1)
    X_scaled = scaler.transform(X)

    # Predict probability
    risk_prob = model.predict_proba(X_scaled)[0, 1]

    # Convert to risk categories
    if risk_prob < 0.1:
        risk_level = 'Low'
        color = 'green'
    elif risk_prob < 0.25:
        risk_level = 'Moderate'
        color = 'yellow'
    elif risk_prob < 0.4:
        risk_level = 'High'
        color = 'orange'
    else:
        risk_level = 'Very High'
        color = 'red'

    return {
        'risk_probability': risk_prob,
        'risk_level': risk_level,
        'color': color,
        'recommendations': generate_recommendations(player_data, risk_level)
    }

def generate_recommendations(player_data, risk_level):
    """
    Generate actionable recommendations based on risk assessment
    """
    recommendations = []

    if player_data['acwr'] > 1.5:
        recommendations.append("ACWR elevated - consider load reduction or rest")

    if player_data['back_to_back'] == 1 and risk_level in ['High', 'Very High']:
        recommendations.append("High risk on back-to-back - recommend rest")

    if player_data['sleep_hours'] < 7:
        recommendations.append("Insufficient sleep - prioritize recovery")

    if player_data['days_since_injury'] < 30:
        recommendations.append("Recent return from injury - monitor closely")

    if player_data['minutes_played'] > 35:
        recommendations.append("High minutes - consider rotation adjustment")

    if not recommendations:
        recommendations.append("Maintain current training load and recovery protocols")

    return recommendations

# Example usage
if __name__ == "__main__":
    # Load and prepare data
    df = load_and_prepare_data()
    df = engineer_features(df)

    # Build model
    model, scaler, features, y_test, y_pred_proba = build_injury_prediction_model(df)

    # Visualize performance
    plot_model_performance(y_test, y_pred_proba)

    # Example: Assess risk for specific player
    player_today = df[df['player_id'] == 'player_001'].iloc[-1]
    risk_assessment = calculate_risk_score(player_today, model, scaler, features)

    print(f"\nPlayer Risk Assessment:")
    print(f"Risk Probability: {risk_assessment['risk_probability']:.1%}")
    print(f"Risk Level: {risk_assessment['risk_level']}")
    print(f"Recommendations:")
    for rec in risk_assessment['recommendations']:
        print(f"  - {rec}")

R: Survival Analysis for Injury Risk

Time-to-Injury Modeling with Cox Proportional Hazards

Survival analysis models the time until an injury event occurs, accounting for players who remain injury-free (censored observations). This approach is particularly valuable for understanding how risk factors influence injury timing.

# Load required libraries
library(survival)
library(survminer)
library(dplyr)
library(ggplot2)
library(tidyr)
library(splines)
library(car)

# Load player tracking and injury data
load_player_data <- function() {
  # Data structure:
  # - player_id: unique identifier
  # - start_date: observation start
  # - end_date: observation end or injury date
  # - injury_event: 1 if injury occurred, 0 if censored (season ended)
  # - age, position, height, weight
  # - avg_minutes_per_game, avg_player_load
  # - acwr_mean, acwr_sd (variability in workload ratio)
  # - previous_injuries (count)
  # - sleep_hours_avg, hrv_avg

  data <- read.csv("player_injury_survival_data.csv")
  return(data)
}

# Calculate time-to-event
prepare_survival_data <- function(data) {
  data <- data %>%
    mutate(
      # Calculate follow-up time in days
      follow_up_days = as.numeric(difftime(end_date, start_date, units = "days")),

      # Risk categories
      acwr_category = case_when(
        acwr_mean < 0.8 ~ "Detraining",
        acwr_mean >= 0.8 & acwr_mean <= 1.3 ~ "Optimal",
        acwr_mean > 1.3 & acwr_mean <= 1.5 ~ "Moderate Risk",
        acwr_mean > 1.5 ~ "High Risk"
      ),
      acwr_category = factor(acwr_category,
                             levels = c("Optimal", "Detraining", "Moderate Risk", "High Risk")),

      # Age categories
      age_group = case_when(
        age < 25 ~ "Young (<25)",
        age >= 25 & age < 30 ~ "Prime (25-29)",
        age >= 30 ~ "Veteran (30+)"
      ),
      age_group = factor(age_group, levels = c("Prime (25-29)", "Young (<25)", "Veteran (30+)")),

      # Workload categories
      high_workload = ifelse(avg_minutes_per_game > 32, "High Load", "Normal Load"),

      # Previous injury history
      injury_history = ifelse(previous_injuries > 0, "Prior Injury", "No Prior Injury")
    )

  return(data)
}

# Fit Cox Proportional Hazards Model
fit_cox_model <- function(data) {
  # Create survival object
  surv_obj <- Surv(time = data$follow_up_days, event = data$injury_event)

  # Fit multivariable Cox model
  cox_model <- coxph(
    surv_obj ~ age + position +
      avg_minutes_per_game + avg_player_load +
      acwr_mean + acwr_sd +
      previous_injuries +
      sleep_hours_avg + hrv_avg,
    data = data
  )

  # Print model summary
  print(summary(cox_model))

  # Test proportional hazards assumption
  ph_test <- cox.zph(cox_model)
  print(ph_test)

  return(cox_model)
}

# Fit model with categorical predictors
fit_cox_categorical <- function(data) {
  surv_obj <- Surv(time = data$follow_up_days, event = data$injury_event)

  cox_cat <- coxph(
    surv_obj ~ age_group + position + acwr_category +
      high_workload + injury_history + sleep_hours_avg,
    data = data
  )

  print(summary(cox_cat))
  return(cox_cat)
}

# Calculate hazard ratios with confidence intervals
extract_hazard_ratios <- function(cox_model) {
  hr_df <- data.frame(
    variable = names(coef(cox_model)),
    HR = exp(coef(cox_model)),
    lower_CI = exp(confint(cox_model)[, 1]),
    upper_CI = exp(confint(cox_model)[, 2]),
    p_value = summary(cox_model)$coefficients[, "Pr(>|z|)"]
  )

  hr_df <- hr_df %>%
    mutate(
      significant = ifelse(p_value < 0.05, "*", ""),
      HR_text = sprintf("%.2f (%.2f-%.2f)%s", HR, lower_CI, upper_CI, significant)
    )

  print("Hazard Ratios (95% CI):")
  print(hr_df %>% select(variable, HR_text, p_value))

  return(hr_df)
}

# Plot survival curves by risk category
plot_survival_curves <- function(data) {
  surv_obj <- Surv(time = data$follow_up_days, event = data$injury_event)

  # Fit survival curves by ACWR category
  fit_acwr <- survfit(surv_obj ~ acwr_category, data = data)

  # Plot with ggsurvplot
  p1 <- ggsurvplot(
    fit_acwr,
    data = data,
    conf.int = TRUE,
    pval = TRUE,
    risk.table = TRUE,
    risk.table.height = 0.25,
    title = "Injury-Free Survival by ACWR Category",
    xlab = "Days",
    ylab = "Probability of Remaining Injury-Free",
    legend.title = "ACWR Category",
    legend.labs = levels(data$acwr_category),
    palette = c("#00BA38", "#619CFF", "#F8766D", "#C77CFF"),
    ggtheme = theme_minimal()
  )

  print(p1)

  # Plot by age group
  fit_age <- survfit(surv_obj ~ age_group, data = data)

  p2 <- ggsurvplot(
    fit_age,
    data = data,
    conf.int = TRUE,
    pval = TRUE,
    risk.table = TRUE,
    risk.table.height = 0.25,
    title = "Injury-Free Survival by Age Group",
    xlab = "Days",
    ylab = "Probability of Remaining Injury-Free",
    legend.title = "Age Group",
    ggtheme = theme_minimal()
  )

  print(p2)

  # Plot by injury history
  fit_history <- survfit(surv_obj ~ injury_history, data = data)

  p3 <- ggsurvplot(
    fit_history,
    data = data,
    conf.int = TRUE,
    pval = TRUE,
    risk.table = TRUE,
    risk.table.height = 0.25,
    title = "Injury-Free Survival by Injury History",
    xlab = "Days",
    ylab = "Probability of Remaining Injury-Free",
    legend.title = "Injury History",
    ggtheme = theme_minimal()
  )

  print(p3)
}

# Create hazard ratio forest plot
plot_hazard_ratios <- function(hr_df) {
  # Filter to significant or notable predictors
  hr_plot <- hr_df %>%
    filter(!is.na(HR)) %>%
    mutate(variable = factor(variable, levels = rev(variable)))

  ggplot(hr_plot, aes(x = HR, y = variable)) +
    geom_vline(xintercept = 1, linetype = "dashed", color = "gray50") +
    geom_point(size = 3) +
    geom_errorbarh(aes(xmin = lower_CI, xmax = upper_CI), height = 0.2) +
    scale_x_log10(breaks = c(0.5, 1, 1.5, 2, 3)) +
    labs(
      title = "Hazard Ratios for Injury Risk Factors",
      subtitle = "Cox Proportional Hazards Model",
      x = "Hazard Ratio (95% CI, log scale)",
      y = ""
    ) +
    theme_minimal() +
    theme(
      panel.grid.major.y = element_blank(),
      plot.title = element_text(face = "bold")
    )

  ggsave("hazard_ratio_forest_plot.png", width = 10, height = 6, dpi = 300)
}

# Predict individual player risk
predict_player_risk <- function(cox_model, player_data) {
  # Calculate linear predictor (log hazard ratio)
  linear_pred <- predict(cox_model, newdata = player_data, type = "lp")

  # Calculate risk score (hazard ratio relative to average)
  risk_score <- predict(cox_model, newdata = player_data, type = "risk")

  # Estimate survival probability at specific time points
  surv_prob_30d <- summary(survfit(cox_model, newdata = player_data), times = 30)$surv
  surv_prob_60d <- summary(survfit(cox_model, newdata = player_data), times = 60)$surv
  surv_prob_90d <- summary(survfit(cox_model, newdata = player_data), times = 90)$surv

  results <- data.frame(
    player_id = player_data$player_id,
    risk_score = risk_score,
    prob_injury_free_30d = surv_prob_30d,
    prob_injury_free_60d = surv_prob_60d,
    prob_injury_free_90d = surv_prob_90d,
    injury_prob_30d = 1 - surv_prob_30d,
    injury_prob_60d = 1 - surv_prob_60d,
    injury_prob_90d = 1 - surv_prob_90d
  )

  return(results)
}

# Time-varying covariates model (advanced)
fit_time_varying_model <- function(data_long) {
  # data_long should have multiple rows per player with time-varying ACWR
  # Requires: tstart, tstop, event, acwr_current, other covariates

  surv_tv <- Surv(time = data_long$tstart,
                  time2 = data_long$tstop,
                  event = data_long$event)

  cox_tv <- coxph(
    surv_tv ~ age + acwr_current + avg_player_load + previous_injuries,
    data = data_long
  )

  print(summary(cox_tv))
  return(cox_tv)
}

# Main analysis workflow
main_analysis <- function() {
  # Load and prepare data
  data <- load_player_data()
  data <- prepare_survival_data(data)

  # Descriptive statistics
  cat("\n=== Descriptive Statistics ===\n")
  cat(sprintf("Total players: %d\n", n_distinct(data$player_id)))
  cat(sprintf("Total injuries: %d (%.1f%%)\n",
              sum(data$injury_event),
              100 * mean(data$injury_event)))
  cat(sprintf("Median follow-up: %.0f days\n", median(data$follow_up_days)))

  # Fit Cox models
  cat("\n=== Cox Proportional Hazards Model (Continuous) ===\n")
  cox_cont <- fit_cox_model(data)

  cat("\n=== Cox Proportional Hazards Model (Categorical) ===\n")
  cox_cat <- fit_cox_categorical(data)

  # Extract and plot hazard ratios
  hr_df <- extract_hazard_ratios(cox_cat)
  plot_hazard_ratios(hr_df)

  # Plot survival curves
  plot_survival_curves(data)

  # Example prediction for high-risk player
  high_risk_player <- data.frame(
    player_id = "PLAYER_001",
    age = 32,
    position = "Guard",
    avg_minutes_per_game = 35,
    avg_player_load = 450,
    acwr_mean = 1.6,
    acwr_sd = 0.4,
    previous_injuries = 2,
    sleep_hours_avg = 6.5,
    hrv_avg = 55
  )

  cat("\n=== Example Risk Prediction ===\n")
  risk_pred <- predict_player_risk(cox_cont, high_risk_player)
  print(risk_pred)

  cat("\nInterpretation:")
  cat(sprintf("\n- Risk score: %.2f (%.0f%% higher risk than average player)",
              risk_pred$risk_score,
              (risk_pred$risk_score - 1) * 100))
  cat(sprintf("\n- Probability of injury in next 30 days: %.1f%%",
              risk_pred$injury_prob_30d * 100))
  cat(sprintf("\n- Probability of injury in next 90 days: %.1f%%\n",
              risk_pred$injury_prob_90d * 100))
}

# Run analysis
main_analysis()

Machine Learning Approaches

Advanced Modeling Techniques

Model Evaluation Considerations

Practical Applications for Teams

Integrated Risk Management System

Successful Implementation Examples

Challenges and Limitations

Ethical Considerations

Player Welfare vs. Team Performance

Future Directions