R Statistical Modeling Best Practices¶

Core Principles¶

Model Selection: Choose appropriate models for the data and question
Validation: Use proper validation techniques to assess model performance
Interpretability: Ensure models are interpretable and explainable
Robustness: Build models that are robust to outliers and violations
Documentation: Document model assumptions, limitations, and results

Linear Models¶

Linear Regression¶

# R/01-linear-models.R

#' Comprehensive linear regression analysis
#'
#' @param data Data frame
#' @param formula Regression formula
#' @param validation_method Validation method
#' @return Linear regression results
linear_regression_analysis <- function(data, formula, validation_method = "holdout") {
  # Fit the model
  model <- lm(formula, data = data)

  # Model summary
  model_summary <- summary(model)

  # Model diagnostics
  diagnostics <- perform_model_diagnostics(model, data)

  # Model validation
  validation_results <- validate_model(model, data, validation_method)

  # Model selection
  selection_results <- perform_model_selection(data, formula)

  results <- list(
    model = model,
    summary = model_summary,
    diagnostics = diagnostics,
    validation = validation_results,
    selection = selection_results
  )

  return(results)
}

#' Perform model diagnostics
#'
#' @param model Fitted model
#' @param data Original data
#' @return Diagnostic results
perform_model_diagnostics <- function(model, data) {
  diagnostics <- list(
    residuals = residuals(model),
    fitted_values = fitted(model),
    leverage = hatvalues(model),
    cooks_distance = cooks.distance(model),
    dffits = dffits(model),
    dfbetas = dfbetas(model)
  )

  # Normality tests
  diagnostics$normality_tests <- list(
    shapiro_test = shapiro.test(diagnostics$residuals),
    anderson_darling = nortest::ad.test(diagnostics$residuals)
  )

  # Heteroscedasticity tests
  diagnostics$heteroscedasticity_tests <- list(
    breusch_pagan = lmtest::bptest(model),
    white_test = lmtest::bptest(model, ~ fitted(model) + I(fitted(model)^2))
  )

  # Autocorrelation tests
  diagnostics$autocorrelation_tests <- list(
    durbin_watson = lmtest::dwtest(model)
  )

  # Multicollinearity
  diagnostics$multicollinearity <- calculate_multicollinearity(model)

  return(diagnostics)
}

#' Calculate multicollinearity measures
#'
#' @param model Fitted model
#' @return Multicollinearity measures
calculate_multicollinearity <- function(model) {
  # Variance Inflation Factors
  vif_values <- car::vif(model)

  # Condition Index
  X <- model.matrix(model)
  eigen_values <- eigen(t(X) %*% X)$values
  condition_index <- sqrt(max(eigen_values) / eigen_values)

  # Tolerance
  tolerance <- 1 / vif_values

  return(list(
    vif = vif_values,
    condition_index = condition_index,
    tolerance = tolerance
  ))
}

#' Validate model using various methods
#'
#' @param model Fitted model
#' @param data Original data
#' @param method Validation method
#' @return Validation results
validate_model <- function(model, data, method = "holdout") {
  if (method == "holdout") {
    return(holdout_validation(model, data))
  } else if (method == "cross_validation") {
    return(cross_validation(model, data))
  } else if (method == "bootstrap") {
    return(bootstrap_validation(model, data))
  }
}

#' Holdout validation
#'
#' @param model Fitted model
#' @param data Original data
#' @return Holdout validation results
holdout_validation <- function(model, data) {
  # Split data
  set.seed(123)
  train_indices <- sample(1:nrow(data), 0.7 * nrow(data))
  train_data <- data[train_indices, ]
  test_data <- data[-train_indices, ]

  # Fit model on training data
  train_model <- lm(formula(model), data = train_data)

  # Predict on test data
  predictions <- predict(train_model, newdata = test_data)
  actual <- test_data[[as.character(formula(model)[[2]])]]

  # Calculate metrics
  mse <- mean((actual - predictions)^2)
  rmse <- sqrt(mse)
  mae <- mean(abs(actual - predictions))
  r_squared <- 1 - sum((actual - predictions)^2) / sum((actual - mean(actual))^2)

  return(list(
    mse = mse,
    rmse = rmse,
    mae = mae,
    r_squared = r_squared,
    predictions = predictions,
    actual = actual
  ))
}

#' Cross-validation
#'
#' @param model Fitted model
#' @param data Original data
#' @param k Number of folds
#' @return Cross-validation results
cross_validation <- function(model, data, k = 10) {
  set.seed(123)
  folds <- create_folds(1:nrow(data), k = k)

  cv_results <- list()

  for (i in 1:k) {
    train_indices <- unlist(folds[-i])
    test_indices <- folds[[i]]

    train_data <- data[train_indices, ]
    test_data <- data[test_indices, ]

    # Fit model
    cv_model <- lm(formula(model), data = train_data)

    # Predict
    predictions <- predict(cv_model, newdata = test_data)
    actual <- test_data[[as.character(formula(model)[[2]])]]

    # Calculate metrics
    mse <- mean((actual - predictions)^2)
    mae <- mean(abs(actual - predictions))

    cv_results[[i]] <- list(mse = mse, mae = mae)
  }

  # Aggregate results
  mean_mse <- mean(sapply(cv_results, function(x) x$mse))
  mean_mae <- mean(sapply(cv_results, function(x) x$mae))
  sd_mse <- sd(sapply(cv_results, function(x) x$mse))
  sd_mae <- sd(sapply(cv_results, function(x) x$mae))

  return(list(
    mean_mse = mean_mse,
    mean_mae = mean_mae,
    sd_mse = sd_mse,
    sd_mae = sd_mae,
    cv_results = cv_results
  ))
}

#' Create k-fold cross-validation folds
#'
#' @param indices Data indices
#' @param k Number of folds
#' @return List of fold indices
create_folds <- function(indices, k = 10) {
  n <- length(indices)
  fold_size <- floor(n / k)
  remainder <- n %% k

  folds <- list()
  start <- 1

  for (i in 1:k) {
    end <- start + fold_size - 1
    if (i <= remainder) {
      end <- end + 1
    }

    folds[[i]] <- indices[start:end]
    start <- end + 1
  }

  return(folds)
}

Generalized Linear Models¶

# R/02-generalized-linear-models.R

#' Comprehensive GLM analysis
#'
#' @param data Data frame
#' @param formula GLM formula
#' @param family GLM family
#' @return GLM analysis results
glm_analysis <- function(data, formula, family = gaussian()) {
  # Fit the model
  model <- glm(formula, data = data, family = family)

  # Model summary
  model_summary <- summary(model)

  # Model diagnostics
  diagnostics <- perform_glm_diagnostics(model, data)

  # Model validation
  validation_results <- validate_glm_model(model, data)

  # Model comparison
  comparison_results <- compare_glm_models(data, formula, family)

  results <- list(
    model = model,
    summary = model_summary,
    diagnostics = diagnostics,
    validation = validation_results,
    comparison = comparison_results
  )

  return(results)
}

#' Perform GLM diagnostics
#'
#' @param model Fitted GLM model
#' @param data Original data
#' @return GLM diagnostic results
perform_glm_diagnostics <- function(model, data) {
  diagnostics <- list(
    residuals = residuals(model),
    fitted_values = fitted(model),
    deviance_residuals = residuals(model, type = "deviance"),
    pearson_residuals = residuals(model, type = "pearson"),
    leverage = hatvalues(model),
    cooks_distance = cooks.distance(model)
  )

  # Deviance analysis
  diagnostics$deviance_analysis <- perform_deviance_analysis(model)

  # Goodness of fit tests
  diagnostics$goodness_of_fit <- perform_goodness_of_fit_tests(model)

  return(diagnostics)
}

#' Perform deviance analysis
#'
#' @param model Fitted GLM model
#' @return Deviance analysis results
perform_deviance_analysis <- function(model) {
  # Null deviance
  null_deviance <- model$null.deviance

  # Residual deviance
  residual_deviance <- model$deviance

  # Deviance explained
  deviance_explained <- (null_deviance - residual_deviance) / null_deviance

  # Deviance residuals
  deviance_residuals <- residuals(model, type = "deviance")

  return(list(
    null_deviance = null_deviance,
    residual_deviance = residual_deviance,
    deviance_explained = deviance_explained,
    deviance_residuals = deviance_residuals
  ))
}

#' Perform goodness of fit tests
#'
#' @param model Fitted GLM model
#' @return Goodness of fit test results
perform_goodness_of_fit_tests <- function(model) {
  # Hosmer-Lemeshow test (for logistic regression)
  if (model$family$family == "binomial") {
    hosmer_lemeshow <- ResourceSelection::hoslem.test(model$y, fitted(model))
  } else {
    hosmer_lemeshow <- NULL
  }

  # Pearson chi-square test
  pearson_chi_square <- sum(residuals(model, type = "pearson")^2)
  pearson_p_value <- 1 - pchisq(pearson_chi_square, model$df.residual)

  return(list(
    hosmer_lemeshow = hosmer_lemeshow,
    pearson_chi_square = pearson_chi_square,
    pearson_p_value = pearson_p_value
  ))
}

#' Compare different GLM models
#'
#' @param data Data frame
#' @param formula Base formula
#' @param family GLM family
#' @return Model comparison results
compare_glm_models <- function(data, formula, family) {
  # Fit different models
  models <- list()

  # Full model
  models$full <- glm(formula, data = data, family = family)

  # Stepwise selection
  models$stepwise <- step(models$full, direction = "both", trace = FALSE)

  # AIC-based selection
  models$aic <- step(models$full, direction = "both", k = 2, trace = FALSE)

  # BIC-based selection
  models$bic <- step(models$full, direction = "both", k = log(nrow(data)), trace = FALSE)

  # Compare models
  comparison <- data.frame(
    model = names(models),
    aic = sapply(models, AIC),
    bic = sapply(models, BIC),
    deviance = sapply(models, function(x) x$deviance),
    df_residual = sapply(models, function(x) x$df.residual)
  )

  return(list(
    models = models,
    comparison = comparison
  ))
}

Mixed Effects Models¶

Linear Mixed Effects¶

# R/03-mixed-effects-models.R

#' Comprehensive linear mixed effects analysis
#'
#' @param data Data frame
#' @param formula Mixed effects formula
#' @param random_effects Random effects specification
#' @return Mixed effects analysis results
lme_analysis <- function(data, formula, random_effects) {
  # Fit the model
  model <- lme4::lmer(formula, data = data, REML = TRUE)

  # Model summary
  model_summary <- summary(model)

  # Model diagnostics
  diagnostics <- perform_lme_diagnostics(model, data)

  # Model validation
  validation_results <- validate_lme_model(model, data)

  # Model comparison
  comparison_results <- compare_lme_models(data, formula, random_effects)

  results <- list(
    model = model,
    summary = model_summary,
    diagnostics = diagnostics,
    validation = validation_results,
    comparison = comparison_results
  )

  return(results)
}

#' Perform LME diagnostics
#'
#' @param model Fitted LME model
#' @param data Original data
#' @return LME diagnostic results
perform_lme_diagnostics <- function(model, data) {
  diagnostics <- list(
    residuals = residuals(model),
    fitted_values = fitted(model),
    random_effects = ranef(model),
    conditional_residuals = residuals(model, type = "response"),
    marginal_residuals = residuals(model, type = "pearson")
  )

  # Normality tests for residuals
  diagnostics$residual_normality <- shapiro.test(diagnostics$residuals)

  # Normality tests for random effects
  diagnostics$random_effects_normality <- lapply(diagnostics$random_effects, function(x) {
    if (ncol(x) == 1) {
      shapiro.test(x[, 1])
    } else {
      NULL
    }
  })

  # Heteroscedasticity tests
  diagnostics$heteroscedasticity <- lmtest::bptest(lm(residuals(model) ~ fitted(model)))

  return(diagnostics)
}

#' Validate LME model
#'
#' @param model Fitted LME model
#' @param data Original data
#' @return LME validation results
validate_lme_model <- function(model, data) {
  # Cross-validation
  cv_results <- cross_validate_lme(model, data)

  # Bootstrap validation
  bootstrap_results <- bootstrap_validate_lme(model, data)

  # Model fit statistics
  fit_statistics <- calculate_lme_fit_statistics(model)

  return(list(
    cross_validation = cv_results,
    bootstrap = bootstrap_results,
    fit_statistics = fit_statistics
  ))
}

#' Cross-validate LME model
#'
#' @param model Fitted LME model
#' @param data Original data
#' @param k Number of folds
#' @return Cross-validation results
cross_validate_lme <- function(model, data, k = 10) {
  set.seed(123)
  folds <- create_folds(1:nrow(data), k = k)

  cv_results <- list()

  for (i in 1:k) {
    train_indices <- unlist(folds[-i])
    test_indices <- folds[[i]]

    train_data <- data[train_indices, ]
    test_data <- data[test_indices, ]

    # Fit model
    cv_model <- lme4::lmer(formula(model), data = train_data, REML = TRUE)

    # Predict
    predictions <- predict(cv_model, newdata = test_data, allow.new.levels = TRUE)
    actual <- test_data[[as.character(formula(model)[[2]])]]

    # Calculate metrics
    mse <- mean((actual - predictions)^2)
    mae <- mean(abs(actual - predictions))

    cv_results[[i]] <- list(mse = mse, mae = mae)
  }

  # Aggregate results
  mean_mse <- mean(sapply(cv_results, function(x) x$mse))
  mean_mae <- mean(sapply(cv_results, function(x) x$mae))
  sd_mse <- sd(sapply(cv_results, function(x) x$mse))
  sd_mae <- sd(sapply(cv_results, function(x) x$mae))

  return(list(
    mean_mse = mean_mse,
    mean_mae = mean_mae,
    sd_mse = sd_mse,
    sd_mae = sd_mae,
    cv_results = cv_results
  ))
}

#' Calculate LME fit statistics
#'
#' @param model Fitted LME model
#' @return Fit statistics
calculate_lme_fit_statistics <- function(model) {
  # AIC and BIC
  aic <- AIC(model)
  bic <- BIC(model)

  # Log-likelihood
  log_lik <- logLik(model)

  # R-squared
  r_squared <- calculate_lme_r_squared(model)

  return(list(
    aic = aic,
    bic = bic,
    log_likelihood = log_lik,
    r_squared = r_squared
  ))
}

#' Calculate LME R-squared
#'
#' @param model Fitted LME model
#' @return R-squared values
calculate_lme_r_squared <- function(model) {
  # Marginal R-squared
  marginal_r_squared <- r.squaredGLMM(model)[1]

  # Conditional R-squared
  conditional_r_squared <- r.squaredGLMM(model)[2]

  return(list(
    marginal = marginal_r_squared,
    conditional = conditional_r_squared
  ))
}

Time Series Models¶

ARIMA Models¶

# R/04-time-series-models.R

#' Comprehensive ARIMA analysis
#'
#' @param data Time series data
#' @param order ARIMA order
#' @param seasonal_order Seasonal ARIMA order
#' @return ARIMA analysis results
arima_analysis <- function(data, order = NULL, seasonal_order = NULL) {
  # Auto-select order if not provided
  if (is.null(order)) {
    order <- auto_arima_order(data)
  }

  # Fit the model
  model <- arima(data, order = order, seasonal = seasonal_order)

  # Model summary
  model_summary <- summary(model)

  # Model diagnostics
  diagnostics <- perform_arima_diagnostics(model, data)

  # Model validation
  validation_results <- validate_arima_model(model, data)

  # Forecasting
  forecast_results <- forecast_arima_model(model, data)

  results <- list(
    model = model,
    summary = model_summary,
    diagnostics = diagnostics,
    validation = validation_results,
    forecast = forecast_results
  )

  return(results)
}

#' Auto-select ARIMA order
#'
#' @param data Time series data
#' @return ARIMA order
auto_arima_order <- function(data) {
  # Use auto.arima from forecast package
  auto_model <- forecast::auto.arima(data, seasonal = FALSE, stepwise = FALSE, approximation = FALSE)

  return(auto_model$arma[1:3])
}

#' Perform ARIMA diagnostics
#'
#' @param model Fitted ARIMA model
#' @param data Original data
#' @return ARIMA diagnostic results
perform_arima_diagnostics <- function(model, data) {
  diagnostics <- list(
    residuals = residuals(model),
    fitted_values = fitted(model),
    acf_residuals = acf(residuals(model), plot = FALSE),
    pacf_residuals = pacf(residuals(model), plot = FALSE)
  )

  # Ljung-Box test
  diagnostics$ljung_box <- Box.test(residuals(model), type = "Ljung-Box")

  # Normality test
  diagnostics$normality <- shapiro.test(residuals(model))

  # ARCH test
  diagnostics$arch_test <- FinTS::ArchTest(residuals(model))

  return(diagnostics)
}

#' Validate ARIMA model
#'
#' @param model Fitted ARIMA model
#' @param data Original data
#' @return ARIMA validation results
validate_arima_model <- function(model, data) {
  # Out-of-sample validation
  n <- length(data)
  train_size <- floor(0.8 * n)

  train_data <- data[1:train_size]
  test_data <- data[(train_size + 1):n]

  # Fit model on training data
  train_model <- arima(train_data, order = model$arma[1:3])

  # Forecast
  forecast_values <- forecast::forecast(train_model, h = length(test_data))

  # Calculate metrics
  mse <- mean((test_data - forecast_values$mean)^2)
  mae <- mean(abs(test_data - forecast_values$mean))
  mape <- mean(abs((test_data - forecast_values$mean) / test_data)) * 100

  return(list(
    mse = mse,
    mae = mae,
    mape = mape,
    forecast_values = forecast_values$mean,
    actual_values = test_data
  ))
}

#' Forecast ARIMA model
#'
#' @param model Fitted ARIMA model
#' @param data Original data
#' @param h Forecast horizon
#' @return Forecast results
forecast_arima_model <- function(model, data, h = 12) {
  # Generate forecasts
  forecast_values <- forecast::forecast(model, h = h)

  # Calculate prediction intervals
  prediction_intervals <- forecast_values$lower[, 1:2]
  prediction_intervals <- cbind(prediction_intervals, forecast_values$upper[, 1:2])
  colnames(prediction_intervals) <- c("lower_80", "lower_95", "upper_80", "upper_95")

  return(list(
    forecast = forecast_values$mean,
    prediction_intervals = prediction_intervals,
    forecast_object = forecast_values
  ))
}

Model Selection and Comparison¶

Model Selection Framework¶

# R/05-model-selection.R

#' Comprehensive model selection
#'
#' @param data Data frame
#' @param response_variable Response variable name
#' @param candidate_models List of candidate models
#' @return Model selection results
model_selection <- function(data, response_variable, candidate_models) {
  # Fit all candidate models
  fitted_models <- fit_candidate_models(data, response_variable, candidate_models)

  # Calculate selection criteria
  selection_criteria <- calculate_selection_criteria(fitted_models)

  # Perform model comparison
  model_comparison <- compare_models(fitted_models, selection_criteria)

  # Select best model
  best_model <- select_best_model(fitted_models, selection_criteria)

  results <- list(
    fitted_models = fitted_models,
    selection_criteria = selection_criteria,
    model_comparison = model_comparison,
    best_model = best_model
  )

  return(results)
}

#' Fit candidate models
#'
#' @param data Data frame
#' @param response_variable Response variable name
#' @param candidate_models List of candidate models
#' @return Fitted models
fit_candidate_models <- function(data, response_variable, candidate_models) {
  fitted_models <- list()

  for (model_name in names(candidate_models)) {
    model_spec <- candidate_models[[model_name]]

    tryCatch({
      if (model_spec$type == "linear") {
        fitted_models[[model_name]] <- lm(model_spec$formula, data = data)
      } else if (model_spec$type == "glm") {
        fitted_models[[model_name]] <- glm(model_spec$formula, data = data, family = model_spec$family)
      } else if (model_spec$type == "lme") {
        fitted_models[[model_name]] <- lme4::lmer(model_spec$formula, data = data)
      }
    }, error = function(e) {
      warning(paste("Failed to fit model", model_name, ":", e$message))
    })
  }

  return(fitted_models)
}

#' Calculate selection criteria
#'
#' @param fitted_models Fitted models
#' @return Selection criteria
calculate_selection_criteria <- function(fitted_models) {
  criteria <- data.frame(
    model = names(fitted_models),
    aic = sapply(fitted_models, AIC),
    bic = sapply(fitted_models, BIC),
    log_likelihood = sapply(fitted_models, function(x) as.numeric(logLik(x))),
    r_squared = sapply(fitted_models, function(x) {
      if (inherits(x, "lm")) {
        summary(x)$r.squared
      } else {
        NA
      }
    }),
    adj_r_squared = sapply(fitted_models, function(x) {
      if (inherits(x, "lm")) {
        summary(x)$adj.r.squared
      } else {
        NA
      }
    }),
    stringsAsFactors = FALSE
  )

  return(criteria)
}

#' Compare models
#'
#' @param fitted_models Fitted models
#' @param selection_criteria Selection criteria
#' @return Model comparison results
compare_models <- function(fitted_models, selection_criteria) {
  comparison <- list()

  # AIC comparison
  comparison$aic <- selection_criteria[order(selection_criteria$aic), ]

  # BIC comparison
  comparison$bic <- selection_criteria[order(selection_criteria$bic), ]

  # Likelihood ratio tests
  comparison$likelihood_ratio_tests <- perform_likelihood_ratio_tests(fitted_models)

  # Information criteria differences
  comparison$aic_differences <- calculate_information_criteria_differences(selection_criteria, "aic")
  comparison$bic_differences <- calculate_information_criteria_differences(selection_criteria, "bic")

  return(comparison)
}

#' Perform likelihood ratio tests
#'
#' @param fitted_models Fitted models
#' @return Likelihood ratio test results
perform_likelihood_ratio_tests <- function(fitted_models) {
  model_names <- names(fitted_models)
  n_models <- length(model_names)

  lrt_results <- data.frame(
    model1 = character(0),
    model2 = character(0),
    lr_statistic = numeric(0),
    p_value = numeric(0),
    stringsAsFactors = FALSE
  )

  for (i in 1:(n_models - 1)) {
    for (j in (i + 1):n_models) {
      model1 <- fitted_models[[model_names[i]]]
      model2 <- fitted_models[[model_names[j]]]

      if (inherits(model1, "lm") && inherits(model2, "lm")) {
        lrt <- anova(model1, model2)
        lr_statistic <- lrt$F[2]
        p_value <- lrt$`Pr(>F)`[2]

        lrt_results <- rbind(lrt_results, data.frame(
          model1 = model_names[i],
          model2 = model_names[j],
          lr_statistic = lr_statistic,
          p_value = p_value,
          stringsAsFactors = FALSE
        ))
      }
    }
  }

  return(lrt_results)
}

#' Calculate information criteria differences
#'
#' @param selection_criteria Selection criteria
#' @param criterion Information criterion
#' @return Information criteria differences
calculate_information_criteria_differences <- function(selection_criteria, criterion) {
  min_value <- min(selection_criteria[[criterion]], na.rm = TRUE)
  differences <- selection_criteria[[criterion]] - min_value

  return(differences)
}

#' Select best model
#'
#' @param fitted_models Fitted models
#' @param selection_criteria Selection criteria
#' @return Best model
select_best_model <- function(fitted_models, selection_criteria) {
  # Select model with minimum AIC
  best_aic_idx <- which.min(selection_criteria$aic)
  best_aic_model <- names(fitted_models)[best_aic_idx]

  # Select model with minimum BIC
  best_bic_idx <- which.min(selection_criteria$bic)
  best_bic_model <- names(fitted_models)[best_bic_idx]

  return(list(
    best_aic = best_aic_model,
    best_bic = best_bic_model,
    best_aic_model = fitted_models[[best_aic_model]],
    best_bic_model = fitted_models[[best_bic_model]]
  ))
}

TL;DR Runbook¶

Quick Start¶

# 1. Linear regression
linear_results <- linear_regression_analysis(data, y ~ x1 + x2 + x3)

# 2. GLM analysis
glm_results <- glm_analysis(data, y ~ x1 + x2, family = binomial())

# 3. Mixed effects
lme_results <- lme_analysis(data, y ~ x1 + x2 + (1|group))

# 4. ARIMA
arima_results <- arima_analysis(ts_data)

# 5. Model selection
candidate_models <- list(
  model1 = list(type = "linear", formula = y ~ x1),
  model2 = list(type = "linear", formula = y ~ x1 + x2)
)
selection_results <- model_selection(data, "y", candidate_models)

Essential Patterns¶

# Model fitting and validation
fit_and_validate <- function(data, formula, model_type = "linear") {
  # Fit model
  if (model_type == "linear") {
    model <- lm(formula, data = data)
  } else if (model_type == "glm") {
    model <- glm(formula, data = data, family = binomial())
  }

  # Validate model
  validation <- validate_model(model, data)

  # Return results
  return(list(model = model, validation = validation))
}