A modular R framework for fitting and validating heteroscedastic regression models of the form:
ΔDES_t = X'β + σ(Z) · ε_t, where σ(Z) = √exp(Z'γ)
This models both the mean (via covariates X) and variance (via covariates Z) of changes in Dietary Energy Supply.
des_model/
├── R/
│ ├── 01_model_core.R # Core model fitting and prediction
│ ├── 02_diagnostics.R # Residual diagnostics and GOF metrics
│ ├── 03_specifications.R # Pre-defined and custom model specs
│ ├── 04_data_preparation.R # Data loading and transformation
│ ├── 05_visualization.R # Coefficient plots and marginal effects
│ └── 06_cross_validation.R # Out-of-sample validation
├── run_analysis.R # Main analysis script
├── data/ # Data directory (add your data here)
└── output/ # Results output
├── figures/
└── tables/
# Load all modules
source("R/01_model_core.R")
source("R/02_diagnostics.R")
source("R/03_specifications.R")
source("R/04_data_preparation.R")
source("R/05_visualization.R")
source("R/06_cross_validation.R")
# Load and prepare data
main_df <- load_main_data("data/your_data.csv")
train <- prepare_training_data(main_df)
# Fit baseline model
model_data <- prepare_model_data(train, spec_full)
fit <- fit_hetero(model_data)
# View results
print(fit)
summary(fit)
# Run diagnostics
diagnose_residuals(fit)
goodness_of_fit(fit)Pre-defined specifications make it easy to test alternative models:
# Available specifications
spec_full # Full model (all covariates)
spec_conflict # Conflict variables only
spec_climate # Climate variables only
spec_economy # Economic variables only
spec_democracy # Democracy variables only
spec_no_quad # No quadratic terms
spec_homoscedastic # Constant variance (no Z covariates)Create custom specifications:
my_spec <- create_spec(
name = "My Custom Model",
x_vars = c("conflict_l1", "warm_d_l1", "econ_d_prop_l1"),
z_vars = c("ytbest_l1", "lgdppc_l1"),
x_labels = c("Intercept", "Conflict", "Warming", "GDP Growth"),
z_labels = c("Intercept", "BRD", "Log GDP")
)
# Or modify an existing spec
modified_spec <- modify_spec(
base_spec = spec_full,
remove_x = c("pop_d_prop_l1"),
add_z = c("secprop_l1"),
name = "Modified Full"
)# Prepare data
model_data <- prepare_model_data(train, spec_full, normalize = TRUE)
# Fit model (with outlier trimming)
fit <- fit_hetero(model_data, trim_quantile = 0.90)
# Access results
fit$coefficients # Coefficient table
fit$beta # Mean equation coefficients
fit$gamma # Variance equation coefficients
fit$loglik # Log-likelihood
fit$vcov # Variance-covariance matrix
# Predictions
pred <- predict(fit) # Returns mu, sigma, lower, upper# Residual diagnostics (creates diagnostic plots)
diag <- diagnose_residuals(fit, plot = TRUE)
print(diag)
# Goodness of fit metrics
gof <- goodness_of_fit(fit)
print(gof)
# Includes: RMSE, MAE, R², AIC, BIC, coverage, CRPS
# Compare multiple models
comparison <- compare_models(fit1, fit2, fit3)
print(comparison)
plot(comparison)
# Likelihood ratio test
lr_test(fit_restricted, fit_full)Heteroscedasticity diagnostics:
# Comprehensive heteroscedasticity tests
het_diag <- diagnose_heteroscedasticity(fit, plot = TRUE)
print(het_diag)
# Output includes:
# - Breusch-Pagan test (H0: homoscedasticity)
# - White's test (H0: homoscedasticity)
# - LR test comparing heteroscedastic vs homoscedastic model
# - Variance ratios across predicted σ quantiles
# - Check if Pearson residuals have constant variance (they should!)
# Visual check: squared residuals by covariate
diagnose_het_by_covariate(fit)Run all standard robustness specifications:
# Get all robustness specs
robust_specs <- standard_robustness_specs()
# Fit all and compare
results <- run_all_specs(robust_specs, train)
print(results$comparison)
plot(results$comparison)Leave-one-out variable analysis:
loo_specs <- create_loo_specs(spec_full)
loo_results <- run_all_specs(loo_specs, train)Subgroup analysis:
# Split by DES level
des_splits <- split_by_des(train, threshold = 2500)
fit_high <- fit_hetero(prepare_model_data(des_splits$high_des, spec_full))
fit_low <- fit_hetero(prepare_model_data(des_splits$low_des, spec_full))
# Split by time period
time_splits <- split_by_year(train, split_year = 2000)DES measurement error robustness (smoothing):
# DES is known to have measurement error. Test robustness to 3-year and 5-year
# moving average smoothing of DES before computing changes.
# Quick: prepare training data with 3-year MA smoothed DES
train_smoothed <- prepare_training_data(main_df, des_smooth = 3)
# Comprehensive: compare raw vs multiple smoothing windows
des_robust <- robustness_des_smoothing(
main_df,
spec_full,
smooth_windows = c(1, 3, 5) # 1 = raw, 3 = 3-year MA, 5 = 5-year MA
)
# View comparison
print(des_robust$comparison)
# Check how a specific coefficient changes with smoothing
compare_smoothing_effect(des_robust, "BRD_MA5")Time-series cross-validation:
cv_results <- ts_cross_validate(train, spec_full, n_folds = 5, test_years = 3)
print(cv_results)Leave-one-country-out:
loco_results <- loco_cross_validate(train, spec_full)
print(loco_results)
plot(loco_results)Compare CV across models:
cv_comparison <- compare_cv(train, robust_specs, cv_type = "ts")
print(cv_comparison)Calibration assessment:
calib <- calibration_assessment(fit)
print(calib)
plot(calib)# Coefficient plots
plot_coefficients(fit, equation = "mean")
plot_dual_coefficients(fit)
# Compare across models
plot_coefficient_comparison(list(fit1, fit2, fit3), equation = "mean")
# Predictions
plot_predictions(fit)
plot_by_covariate(fit, "conflict_l1")
# Marginal effects
me <- marginal_effect(fit, "conflict_l1")
plot(me)
# Average marginal effects
ame <- average_marginal_effect(fit, "best", marginal_values = seq(-2000, 2000, 100))
plot(ame)
# Save all plots for a model
save_model_plots(fit, output_dir = "output/figures")-
Fit baseline model
fit_baseline <- fit_hetero(prepare_model_data(train, spec_full))
-
Run diagnostics
diagnose_residuals(fit_baseline) goodness_of_fit(fit_baseline)
-
Test alternative specifications
robust_results <- run_all_specs(standard_robustness_specs(), train) print(robust_results$comparison)
-
Likelihood ratio tests
lr_test(robust_results$fits$homoscedastic, robust_results$fits$full) lr_test(robust_results$fits$no_quad, robust_results$fits$full)
-
Cross-validation
cv_ts <- ts_cross_validate(train, spec_full) cv_loco <- loco_cross_validate(train, spec_full)
-
Calibration
calib <- calibration_assessment(fit_baseline) plot(calib)
-
Subgroup analysis
# By DES level, time period, region, etc. -
Export results
coef_tables <- coefficient_table(robust_results$fits) fwrite(robust_results$comparison$comparison, "output/model_comparison.csv")
If you're primarily interested in prediction (not causal inference), several alternatives may outperform the baseline heteroscedastic model:
| Priority | Model | Why |
|---|---|---|
| Interpretability + uncertainty | Current (hetero) or GAMLSS | Clear coefficients |
| Best point predictions | Gradient Boosting | Captures nonlinearities |
| Robust to outliers | Quantile Regression | Median more robust |
| Full predictive distribution | Quantile RF | Nonparametric intervals |
| Uncertainty quantification | Bayesian (brms) | Full posterior |
Dedicated module for comparing the two main approaches:
| Aspect | Heteroscedastic Model | Quantile Regression |
|---|---|---|
| Assumption | Normal errors with varying variance | Distribution-free |
| Point prediction | Conditional mean | Conditional median (or any quantile) |
| Intervals | Based on estimated σ(Z) | Direct quantile estimates |
| Robustness | Sensitive to outliers | Robust to outliers |
| Coefficients | One set for mean, one for variance | Different coefficients at each quantile |
| Interpretation | "Effect on average change" | "Effect at different parts of distribution" |
# Full comparison with cross-validation
comparison <- full_hetero_qreg_comparison(train, spec = spec_full, run_cv = TRUE)
# View summary
print(comparison)
# Access individual components
comparison$predictions$metrics # In-sample metrics
comparison$cv_results$summary # CV metrics
comparison$coefficients$qreg # QR coefficients at each quantile
# Plots
comparison$plots$coefficients # Coefficients across quantiles
comparison$plots$intervals # Interval comparison
comparison$plots$width # Interval width comparison
comparison$plots$cv # CV results boxplotThe coefficient comparison plot shows QR coefficients (blue) across quantiles τ = 0.05, 0.10, ..., 0.95, with the heteroscedastic mean coefficient (dashed red line) for comparison.
Interpretation:
- Flat blue line ≈ red line: Effects are homogeneous across the distribution; heteroscedastic model is appropriate
- Sloped blue line: Effects differ at different quantiles (e.g., conflict affects food-insecure countries more than food-secure ones)
- Blue line crosses zero: Effect changes sign across the distribution
Prefer Heteroscedastic when:
- Residuals are approximately normal
- You want an interpretable variance model
- Prediction intervals should reflect modeled uncertainty drivers
Prefer Quantile Regression when:
- Heavy-tailed or skewed residuals
- Outliers are common
- You want to understand heterogeneous effects across the distribution
- You don't want to assume a parametric distribution
# 1. Quantile Regression - robust, no distributional assumptions
qr_fit <- fit_quantile_reg(train)
qr_pred <- predict(qr_fit, test)
plot_qreg_coefficients(qr_fit) # How effects vary by quantile
# 2. GAMLSS - flexible parametric, add smooth terms
gamlss_fit <- fit_gamlss(train)
gamlss_pred <- predict_gamlss(gamlss_fit, test)
# 3. GAMLSS with smooth terms (nonlinear relationships)
gamlss_smooth <- fit_gamlss_smooth(train, smooth_vars = c("conflict_l1", "lgdppc_l1"))
# 4. Quantile Random Forest - nonparametric prediction intervals
qrf_fit <- fit_qrf(train)
qrf_pred <- predict(qrf_fit, test)
# 5. Gradient Boosting with quantile loss
gbm_fit <- fit_gbm_quantile(train)
gbm_pred <- predict(gbm_fit, test)
# 6. Bayesian heteroscedastic regression (requires brms/Stan)
bayes_fit <- fit_bayesian_hetero(train)
bayes_pred <- predict_bayesian(bayes_fit, test)# In-sample comparison
all_preds <- list(
hetero = predict(fit_full),
gamlss = gamlss_pred,
qreg = qr_pred,
qrf = qrf_pred
)
compare_prediction_models(all_preds, train$des_d)
# Cross-validation comparison
cv_results <- cv_compare_models(train, c("hetero", "gamlss", "qreg", "qrf"))
print(cv_results$summary)# For alternative models
install.packages(c("gamlss", "quantreg", "quantregForest", "lightgbm", "brms"))| Variable | Description |
|---|---|
conflict_l1 |
Rolling sum of transformed BRD (5-year) |
warm_d_l1 |
Change in TX90 (5-year) |
econ_d_prop_l1 |
Proportional GDP/capita change (3-year) |
dem_d_l1 |
Change in democracy index (10-year) |
dem_d_2_l1 |
Squared democracy change |
pop_d_prop_l1 |
Proportional population change (5-year) |
| Variable | Description |
|---|---|
ytbest_l1 |
Yeo-Johnson transformed BRD |
tx90pgs_l1 |
TX90 (days >90th percentile temp) |
lgdppc_l1 |
Log GDP per capita |
v2x_polyarchy_l1 |
Electoral Democracy Index |
v2x_polyarchy_2_l1 |
Squared democracy index |
lpopulation_l1 |
Log population |
# Core dependencies
install.packages(c(
"data.table", "tidyverse", "tsibble", "slider", "zoo",
"ggplot2", "patchwork", "ggdist",
"moments", "MASS"
))
# For quantile regression comparison
install.packages("quantreg")
# For alternative models (optional)
install.packages(c("gamlss", "quantregForest", "lightgbm", "brms"))- All predictors are standardized by default for comparability
- Outlier trimming (default: 90th percentile) improves robustness
- The model uses
nlm()for optimization with analytical Hessian - Cross-validation uses expanding window for time-series