---
title: "Model Diagnostics and Residual Analysis"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Model Diagnostics and Residual Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  eval = TRUE
)
```

## Introduction

In Small Area Estimation (SAE), model-based estimators borrow strength across domains and time to improve estimation precision over direct survey estimators. However, introducing model assumptions carries the risk of model misspecification, bias, and inappropriate shrinkage.

To address these challenges, **fastsae** provides a unified diagnostic framework centered on the `diagnose()` function and its dedicated S3 `autoplot()` visualization method. It evaluates both **frequentist C++ EBLUP** models (`eblup_fh`, `eblup_sfh`, `eblup_stfh`, `eblup_bhf`) and **Bayesian INLA HB** models (`hb_area`).

---

## Core Diagnostic Criteria

The diagnostic workflow evaluates four fundamental dimensions of SAE model quality:

1. **Precision Gain & RSE Reduction**:
   Quantifies the improvement in Relative Standard Error ($\text{RSE} = \frac{\sqrt{\text{MSE}}}{\hat{\theta}} \times 100\%$) compared to direct survey estimates. Evaluates whether areas satisfy the official statistics reliability threshold (typically $\text{RSE} < 25\%$ or $20\%$).
2. **Goodness-of-Fit Statistic (Brown et al., 2001)**:
   Computes the Chi-square diagnostic statistic $W$ over the sampled domains:
   $$W = \sum_{d=1}^D \frac{(y_d - \hat{\theta}_d)^2}{V_d + \text{MSE}(\hat{\theta}_d)}$$
   and compares it with $\chi^2(D)$. Only large values indicate that the model estimates differ from the direct estimates by more than their uncertainty, so the p-value is upper-tailed (as in `emdi`). The statistic ignores the positive correlation between $y_d$ and $\hat{\theta}_d$, so it is conservative: under a correct model $W$ is usually well below $D$.
3. **Bias Diagnostic (Brown et al., 2001)**:
   The regression of the direct estimates on the model estimates,
   $$y_d = \alpha + \beta \hat{\theta}_d + e_d,$$
   is reported descriptively (with `autoplot(type = "calibration")`) to look for systematic deviations. It is not used as a formal test of $\alpha = 0, \beta = 1$: for shrinkage estimators such as the EBLUP the slope exceeds 1 even when the model is correct (about $1/\gamma_d$ for the variation not explained by the covariates), so such a test would reject correct models.
4. **Residual Spatial Autocorrelation (Moran's I)**:
   Verifies whether spatial autocorrelation has been fully captured by the model random effects. Significant residual spatial clustering indicates that unmodeled spatial patterns remain.

---

## Practical Example: Diagnosing a Fay-Herriot Model

We demonstrate the diagnostic workflow using the built-in `mys` dataset (Mean Years of Schooling in 42 regencies):

### 1. Fit Model and Run `diagnose()`

```{r diag_fit, message = FALSE, warning = FALSE}
library(fastsae)
library(ggplot2)

data("mys")
data("mys_proxmat")

# Fit Area-Level Fay-Herriot model
fit_fh <- eblup_fh(
  formula = y ~ x1 + x2 + x3,
  vardir = ~vardir,
  data = na.omit(mys),
  method = "REML",
  print_result = FALSE
)

# Run comprehensive diagnostic checks
diag_fh <- diagnose(fit_fh, W = mys_proxmat[!is.na(mys$y), !is.na(mys$y)])
print(diag_fh)
```

### 2. Inspect Diagnostic Components

The resulting `fastsae_diagnose` object provides structured access to each diagnostic domain:

```{r diag_precision}
# Efficiency gains
diag_fh$precision

# Calibration test results
diag_fh$brown_test

# Moran's I spatial test on residuals
diag_fh$spatial_test
```

---

## Visual Diagnostics with `autoplot()`

`fastsae` implements dedicated ggplot2 visualization types for `fastsae_diagnose` objects:

### 1. Calibration Plot (`type = "calibration"`)

Plots direct survey estimates against model predictions along with the 1:1 identity line and the fitted calibration regression:

```{r plot_calib}
autoplot(diag_fh, type = "calibration")
```

### 2. Precision Comparison (`type = "rse"`)

Visualizes the domain-by-domain reduction in Relative Standard Error:

```{r plot_rse}
autoplot(diag_fh, type = "rse")
```

### 3. Residual Diagnostics (`type = "residuals"` & `type = "qq"`)

Inspect standardized residuals against fitted values to detect heteroscedasticity, non-linearities, or outliers:

```{r plot_res}
autoplot(diag_fh, type = "residuals")
autoplot(diag_fh, type = "qq")
```

### 4. Combined Diagnostic Dashboard (`type = "all"`)

Generates a unified diagnostic view summarizing calibration and precision gains:

```{r plot_all}
autoplot(diag_fh, type = "all")
```

---

## Diagnosing Bayesian Spatio-Temporal Models

The exact same `diagnose()` function applies directly to Bayesian models fitted with `hb_area()`:

```{r diag_hb, message = FALSE, warning = FALSE, eval = requireNamespace("INLA", quietly = TRUE)}
data("sim_area", package = "fastsae")

# Fit Spatial EBP model with BYM2 prior
fit_hb <- hb_area(
  formula = y_gaussian ~ x1 + x2,
  data = sim_area,
  vardir = "vardir",
  spatial = "bym2",
  W = create_weights(cbind(sim_area$x_coord, sim_area$y_coord), method = "knn", k = 4),
  print_result = FALSE
)

# Diagnose Bayesian model
diag_hb <- diagnose(fit_hb)
diag_hb$precision
diag_hb$status
```

---

## References

- Brown, G., Chambers, R., Heady, P., & Heasman, D. (2001). Evaluation criteria for small area estimation methods. *Statistics in Transition*, 5(2), 185–200.
- Fay, R. E., & Herriot, R. A. (1979). Estimates of income for small places: An application of James-Stein procedures to Census data. *Journal of the American Statistical Association*, 74(366), 269–277.
- Rao, J. N. K., & Molina, I. (2015). *Small Area Estimation* (2nd ed.). John Wiley & Sons.