---
title: "Getting Started with saeHB.Spatial.Beta"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting Started with saeHB.Spatial.Beta}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  message = FALSE
)
```

## Introduction

The `saeHB.Spatial.Beta` package provides several functions to estimate small area proportions using Hierarchical Bayesian (HB) methods under spatial and non-spatial models for area-level random effects. This package is specifically designed to accommodate survey design effects (DEFF) for sampling variances.

In this vignette, we will demonstrate a complete analytical workflow:

1. **Preparation and Data Loading**, including calculating the variance and RSE of the direct estimates to serve as a baseline.
2. **Fitting the HB Beta Deff Non-Spatial Model**.
3. **Spatial Autocorrelation Diagnostic** on the random effects using Moran's I.
4. **Fitting the HB Beta Deff Spatial Models** (SAR and Leroux CAR).
5. **Evaluation and Comparison** of estimation precision using Relative Standard Error (RSE).

## Model Specifications

Before diving into the workflow, let's briefly review the underlying models. Let $\hat{\theta}_i$ denote the direct estimator of the proportion for area $i$ ($i = 1, \dots, m$), and $\theta_i$ denote the true proportion parameter. Under the Hierarchical Bayesian framework, the Beta spatial model accommodating survey design effects is specified as follows:

**1. Sampling Model:**
$$\hat{\theta}_i \mid \theta_i \sim \text{Beta}(a_i, b_i)$$
To accommodate the effective sample size from complex surveys, the shape parameters are strictly defined using the area's sample size ($n_i$) and survey design effect ($d_i$). Thus, the parameters are formulated as:
$$a_i = \theta_i \left(\frac{n_i}{d_i} - 1\right) \quad \text{and} \quad b_i = (1 - \theta_i) \left(\frac{n_i}{d_i} - 1\right)$$
With this specification, the expected value holds $E(\hat{\theta}_i \mid \theta_i) = \theta_i$. The sampling variance of the direct estimator can be expressed as $\psi_i = \left[ \frac{\theta_i (1-\theta_i)}{n_i} \right] d_i$, indicating that the sampling variance directly depends on the proportion parameter $\theta_i$ and the design effect $d_i$.

**2. Linking Model:**
The true proportion $\theta_i$ is modeled using a logit link function:
$$\text{logit}(\theta_i) = \mathbf{x}_i^T \boldsymbol{\beta} + v_i$$
Here, $\mathbf{x}_i^T$ represents the vector of auxiliary variables, $\boldsymbol{\beta}$ the regression coefficients, and $v_i$ the area-specific random effect.

**3. Spatial Random Effects:**
To account for spatial dependency, the random effects vector $\mathbf{v} = (v_1, \dots, v_m)^T$ can be modeled using two different spatial structures:

a.  **SAR Model:** The random effects follow a Simultaneous Autoregressive process, mathematically defined as:

    $$\mathbf{v} = \rho \mathbf{W} \mathbf{v} + \mathbf{u}, \quad \mathbf{u} \sim N_m(\mathbf{0}, \sigma_u^2 \mathbf{I})$$

    Assuming the matrix $(\mathbf{I} - \rho \mathbf{W})$ is non-singular, it can be rewritten explicitly as $\mathbf{v} = (\mathbf{I} - \rho \mathbf{W})^{-1} \mathbf{u}$, which implies $\mathbf{v} \sim N_m(\mathbf{0}, \mathbf{G})$ with the covariance dispersion matrix defined as:

    $$\mathbf{G} = \sigma_u^2 [(\mathbf{I} - \rho \mathbf{W})^T (\mathbf{I} - \rho \mathbf{W})]^{-1}$$

    Where $\rho$ is the spatial autocorrelation parameter, $\mathbf{W}$ is a row-standardized spatial weights matrix, and $\mathbf{u}$ is the independent random error vector.

b.  **Leroux CAR Model:** The random effects follow a Conditional Autoregressive structure, $\mathbf{v} \sim N_m(\mathbf{0}, \mathbf{Q}(\rho)^{-1})$, with the precision matrix defined as:

    $$\mathbf{Q}(\rho) = \frac{1}{\sigma_v^2} [(1 - \rho)\mathbf{I} + \rho(\mathbf{D} - \mathbf{W})]$$

    Where $\mathbf{W}$ is the binary spatial adjacency matrix, $\mathbf{D}$ is a diagonal matrix containing the number of neighbors for each area, and $\rho \in [0,1)$ controls the spatial dependence.

## Step 1: Preparation and Data Loading

First, load the package along with the provided synthetic dataset (`databeta`). We will also load `ggplot2` for visualization. After loading the data, we calculate the variance and the Relative Standard Error (RSE) of the direct estimates to serve as our baseline for comparison.

```{r setup}
library(saeHB.Spatial.Beta)
library(ggplot2)

# Load data
data("databeta")

# Calculate Variance of Direct Estimator for proportion data considering DEFF
# var(y) = [y * (1 - y) / n_i] * deff
var_direct <- (databeta$y * (1 - databeta$y) / databeta$n_i) * databeta$deff

# Calculate Relative Standard Error (RSE) of Direct Estimation
databeta$rse_direct <- (sqrt(var_direct) / databeta$y) * 100
```

## Step 2: Fitting the HB Beta Deff Non-Spatial Model

We begin by fitting a baseline non-spatial model that accommodates the survey design effect. Here, we use the default Markov Chain Monte Carlo (MCMC) iterations.

*Note: By default, the function runs MCMC with a predefined number of iterations and burn-in periods. For practical applications, you may need to adjust parameters such as iter.mcmc, burn.in, thin, and chains to ensure proper mixing and convergence.*

```{r ns-model, results='hide'}
mod_ns_deff <- betadeff_nonspatial(
  formula = y ~ x1 + x2,
  deff = "deff",
  n_i = "n_i",
  data = databeta
)
```

Extract the RSE for the non-spatial estimates:

```{r ns-rse}
rse_ns_deff <- (mod_ns_deff$est$Est.Error / mod_ns_deff$est$Estimate) * 100
```

## Step 3: Spatial Autocorrelation Diagnostic (Moran's I)

To determine whether a spatial model is warranted, we evaluate the spatial autocorrelation of the random effects ($v$) obtained from the non-spatial model. We use a row-standardized spatial weight matrix for Moran's I testing.

```{r moran-weights}
data("weight_mat")
W_listw <- spdep::mat2listw(weight_mat, style = "W")

# Extract the mean of the random effects (v)
v_ns_deff <- as.numeric(mod_ns_deff$randeff$Estimate)
```

Since our dataset has a relatively small number of areas ($m = 36$), we use the Monte Carlo permutation approach with 999 permutations, which computes the p-value empirically by randomly permuting the observed values. The analytical randomisation approach can be used as an alternative (`mc = FALSE`).

```{r moran-test}
set.seed(123) 
moran_result <- moran_test(x = v_ns_deff, listw = W_listw, mc = TRUE, nsim = 999)
print(moran_result)
```

A significant p-value confirms that the non-spatial model left unexplained spatial structure, heavily justifying the use of spatial models.

## Step 4: Fitting the HB Beta Deff Spatial Models

We will now fit two spatial models: the Simultaneous Autoregressive (SAR) model and the Leroux Conditional Autoregressive (CAR) model. The SAR model requires a row-standardized weight matrix, while the Leroux CAR model requires a binary adjacency matrix.

```{r spatial-models, results='hide'}
# 1. Fit Spatial SAR Model
# Load the spatial weight matrix for the SAR model
data("weight_mat")
mod_sar_deff <- betadeff_sar(
  formula = y ~ x1 + x2,
  deff = "deff",
  n_i = "n_i",
  proxmat = weight_mat,
  data = databeta
)
# 2. Fit Spatial Leroux CAR Model
# Load the binary adjacency matrix for the Leroux CAR model
data("adjacency_mat")
mod_leroux_deff <- betadeff_lerouxcar(
  formula = y ~ x1 + x2,
  deff = "deff",
  n_i = "n_i",
  proxmat = adjacency_mat,
  data = databeta
)
```

Extract the RSE for both spatial models:

```{r spatial-rse}
rse_sar_deff <- (mod_sar_deff$est$Est.Error / mod_sar_deff$est$Estimate) * 100
rse_leroux_deff <- (mod_leroux_deff$est$Est.Error / mod_leroux_deff$est$Estimate) * 100
```

## Step 5: Evaluation and Comparison

We compare the performance of the Direct Estimator, the HB Beta Deff Non-Spatial Model, and the HB Beta Deff Spatial Models (SAR and Leroux CAR). A lower RSE indicates a more reliable and precise estimate.

First, let's look at the summary statistics of the RSEs provided by the SAE models compared to the direct estimates.

```{r rse-summary, echo=FALSE}
calc_stats <- function(x) {
  c(
    Minimum = min(x, na.rm = TRUE),
    `First Quartile` = unname(quantile(x, 0.25, na.rm = TRUE)),
    Median = median(x, na.rm = TRUE),
    Mean = mean(x, na.rm = TRUE),
    `Third Quartile` = unname(quantile(x, 0.75, na.rm = TRUE)),
    Maximum = max(x, na.rm = TRUE)
  )
}

rse_summary <- data.frame(
  Direct = calc_stats(databeta$rse_direct),
  `HB Beta Deff Non-Spatial` = calc_stats(rse_ns_deff),
  `HB Beta Deff Spatial Leroux CAR` = calc_stats(rse_leroux_deff),
  `HB Beta Deff Spatial SAR` = calc_stats(rse_sar_deff),
  check.names = FALSE
)

knitr::kable(
  rse_summary, 
  digits = 2, 
  align = "c",
  caption = "Descriptive Statistics of Relative Standard Error (RSE) (%)"
)
```

We can also visualize this comparison across all areas using a line and point plot to observe the fluctuation of RSE values.

```{r comparison, fig.width=8, fig.height=5}
# Combine RSEs into a single data frame for plotting
df_rse <- data.frame(
  Area = seq_along(databeta$y),
  Direct = databeta$rse_direct,
  Non_Spatial = rse_ns_deff,
  Spatial_SAR = rse_sar_deff,
  Spatial_Leroux = rse_leroux_deff
)

# Order by Direct RSE for better visualization
df_rse <- df_rse[order(df_rse$Direct), ]
df_rse$Area_Index <- seq_len(nrow(df_rse))

# Plotting the RSE Comparison
ggplot(df_rse, aes(x = Area_Index)) +
  # Direct Estimation
  geom_line(aes(y = Direct, color = "Direct"), linewidth = 0.8, alpha = 0.6) +
  geom_point(aes(y = Direct, color = "Direct"), size = 2, alpha = 0.6) +
  
  # Non-Spatial
  geom_line(aes(y = Non_Spatial, color = "HB Beta Deff Non-Spatial"), linewidth = 0.8, alpha = 0.8) +
  geom_point(aes(y = Non_Spatial, color = "HB Beta Deff Non-Spatial"), size = 2, alpha = 0.8) +
  
  # Spatial SAR
  geom_line(aes(y = Spatial_SAR, color = "HB Beta Deff Spatial SAR"), linewidth = 1) +
  geom_point(aes(y = Spatial_SAR, color = "HB Beta Deff Spatial SAR"), size = 2) +
  
  # Spatial Leroux CAR
  geom_line(aes(y = Spatial_Leroux, color = "HB Beta Deff Spatial Leroux CAR"), linewidth = 1) +
  geom_point(aes(y = Spatial_Leroux, color = "HB Beta Deff Spatial Leroux CAR"), size = 2) +
  
  scale_color_manual(
    name = "Estimator",
    values = c("Direct" = "#E69F00", 
               "HB Beta Deff Non-Spatial" = "#56B4E9", 
               "HB Beta Deff Spatial SAR" = "#009E73",
               "HB Beta Deff Spatial Leroux CAR" = "#D55E00")
  ) +
  labs(
    title = "Comparison of Relative Standard Error (RSE)",
    subtitle = "Lower RSE indicates higher precision",
    x = "Area (Ordered by Direct RSE)",
    y = "RSE (%)"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")
```