---
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. This package is specifically designed to accommodate survey design effects (DEFF) and handle non-sampled areas.

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

1. Calculating Direct Estimates.
2. Fitting the HB Beta Deff Non-Spatial Model.
3. Checking spatial autocorrelation on the random effects using Moran's I.
4. Fitting the HB Beta Deff Spatial Models (SAR and Leroux CAR).
5. Comparing the 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$$
Or equivalently in a conditional distribution form: $\text{logit}(\theta_i) \mid \boldsymbol{\beta}, \sigma_v^2 \sim N(\mathbf{x}_i^T \boldsymbol{\beta}, \sigma_v^2)$. Here, $\mathbf{x}_i^T$ represents the vector of auxiliary variables, $\boldsymbol{\beta}$ the regression coefficients, and $v_i$ the area-specific spatial 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 real-world applications with complex datasets, you may need to adjust parameters such as iter.mcmc, burn.in, and thin to ensure 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 if a spatial model is necessary, we evaluate the spatial autocorrelation of the random effects ($v$) from the non-spatial model. 

We use the **Monte Carlo Permutation** approach. This method is highly recommended because it computes the p-value empirically by shuffling the spatial units, making it robust. The Analytical approach (randomization) can be used as an alternative (`mc = FALSE`) for very large datasets where computational speed is a priority.

```{r moran-test}
# Load the spatial weight matrix for the diagnostic test
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)

# Monte Carlo Permutation Approach
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 Non-Spatial Model, and the Spatial Models. 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),
  `Non-Spatial` = calc_stats(rse_ns_deff),
  `Spatial SAR` = calc_stats(rse_sar_deff),
  `Spatial Leroux` = calc_stats(rse_leroux_deff),
  check.names = FALSE
)

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

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"), linewidth = 1) +
  geom_point(aes(y = Spatial_Leroux, color = "HB Beta Deff Spatial Leroux"), 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" = "#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")
```