CorrelVariogram_pt2.rmd

---
title: "Spatial Dependence, Part 2: Correlograms & Semivariograms with CSV Data (in R)"
author: "Daniel M. Parker (w/ help from ChatGPT)"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: false
    number_sections: false
    code_folding: "hide"
    theme: cosmo
    highlight: textmate
---

```{r setup, message=FALSE, warning=FALSE}
knitr::opts_chunk$set(message = FALSE, warning = FALSE)
library(sf)
library(sp)
library(spdep)
library(ggplot2)
library(gstat)
```

# 1) Overview

This Part 2 tutorial shows how to compute **spatial correlograms** (Moran’s I; Geary’s C) and a **semivariogram** from a simple **CSV** of points:

- CSV must contain at least: **X**, **Y**, and **one numeric value column**.
- This repo includes two dummy datasets:
  - `dummy_strong_spatial.csv` with **COVIDincid** (strictly positive; strong spatial trend)
  - `dummy_no_spatial.csv` with **KuruPrev** (weak/no spatial trend)

We’ll:
1. Load and project the data
2. Choose **distance bins** (max lag = **1/3 of the maximum pairwise distance**)
3. Compute **Moran’s I** and **Geary’s C** by distance ring
4. Compute a **semivariogram** with `gstat`
5. Plot and compare

# 2) Load a CSV and project to meters

> Set `file_name` to either `"dummy_strong_spatial.csv"` (COVIDincid) or `"dummy_no_spatial.csv"` (KuruPrev).  
> If you use your own CSV, ensure it has columns `X`, `Y`, and a numeric **value** column.

```{r}
# Choose which CSV to use
file_name <- "dummy_strong_spatial.csv"   # or "dummy_no_spatial.csv" (or your own path)

# Read CSV
df <- read.csv(file_name)

# Inspect columns
str(df)
head(df)

# Convert lon/lat (EPSG:4326) to sf POINTS
pts <- st_as_sf(df, coords = c("X","Y"), crs = 4326)

# Project to a suitable projected CRS in meters (adjust EPSG as needed)
# Example: UTM zone 33N
pts <- st_transform(pts, 32633)

# Choose which value column to analyze:
# If using the included datasets, value column is named COVIDincid or KuruPrev.
candidate_cols <- intersect(c("COVIDincid","KuruPrev","Value"), names(pts))
stopifnot(length(candidate_cols) >= 1)
value_col <- candidate_cols[1]
value_col

coords <- st_coordinates(pts)     # matrix of projected x,y (meters)
x      <- pts[[value_col]]        # numeric variable
summary(x)
```

# 3) Define distance bins

We’ll pick bins by:
- computing the **maximum pairwise distance** among points,
- setting **max lag = 1/3** of that maximum,
- choosing a **fixed number of bins** (e.g., 10).

Then we’ll visualize the **distribution of pairwise distances** with a histogram and use the bin edges as breaks.

```{r}
# Pairwise Euclidean distances (meters)
dists <- as.numeric(dist(coords))

# Max lag and bins
max_distance <- max(dists) / 3
n_bins <- 10
edges  <- seq(0, max_distance, length.out = n_bins + 1)

bins <- data.frame(
  bin_id = 1:n_bins,
  d_min  = head(edges, -1),
  d_max  = tail(edges, -1)
)
bins
```

### Histogram of pairwise distances with bin edges

```{r, fig.height=3.8}
### Histogram of pairwise distances with chosen bin edges overlaid
hist(dists,
     main = "Pairwise distances",
     xlab = "Distance (meters)",
     ylab = "Pair count")
abline(v = edges, lty = 3)  # draw the chosen bin edges as vertical lines
```

# 4) Build binary weights for a distance ring

Here we define how to treat pairs of points within a given distance bin: if two points fall inside the same ring, we code them as **neighbors (1)**; if they fall outside the ring (across bins), we code them as **not neighbors (0)**. This creates a binary weights list for each ring.

```{r}
# Build binary weights for a distance ring [d_min, d_max)
ring_listw <- function(coords, d_min, d_max, zero_ok = TRUE) {
  # neighbors if separated by a distance within [d_min, d_max)
  nb <- spdep::dnearneigh(coords, d1 = d_min, d2 = d_max, longlat = FALSE)
  lw <- spdep::nb2listw(nb, style = "B", zero.policy = zero_ok)
  list(nb = nb, lw = lw)
}
```

# 5) Moran’s I correlogram (with permutation p-values)

```{r}
moran_ring <- function(x, coords, d_min, d_max, nsim = 499, zero_ok = TRUE) {
  w <- ring_listw(coords, d_min, d_max, zero_ok)
  # If no neighbor links exist in this ring, return NAs
  if (sum(unlist(w$lw$weights), na.rm = TRUE) == 0) {
    return(list(I = NA_real_, p_mc = NA_real_, n_links = 0))
  }
  mc <- spdep::moran.mc(x, w$lw, nsim = nsim,
                        zero.policy = zero_ok, alternative = "two.sided")
  list(I = as.numeric(mc$statistic),
       p_mc = mc$p.value,
       n_links = sum(unlist(w$lw$weights)))
}

set.seed(123)
res_I <- do.call(rbind, lapply(seq_len(nrow(bins)), function(i) {
  b <- bins[i,]
  r <- moran_ring(x, coords, b$d_min, b$d_max, nsim = 499, zero_ok = TRUE)
  data.frame(bin_id = b$bin_id, d_min = b$d_min, d_max = b$d_max,
             I = r$I, p_mc = r$p_mc, n_links = r$n_links)
}))
res_I$mid <- (res_I$d_min + res_I$d_max)/2
res_I
```

### Plot Moran’s I correlogram

```{r, fig.height=4}
ggplot(res_I, aes(mid, I)) +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_point() + geom_line() +
  labs(x = "Distance (bin midpoint, m)", y = "Moran's I",
       title = paste0("Moran's I by distance (", value_col, ")")) +
  theme_minimal()
```

# 6) Geary’s C correlogram (optional)

```{r}
geary_ring <- function(x, coords, d_min, d_max, nsim = 499, zero_ok = TRUE) {
  w <- ring_listw(coords, d_min, d_max, zero_ok)
  if (sum(unlist(w$lw$weights), na.rm = TRUE) == 0) {
    return(list(C = NA_real_, p_mc = NA_real_, n_links = 0))
  }
  mc <- spdep::geary.mc(x, w$lw, nsim = nsim,
                        zero.policy = zero_ok, alternative = "two.sided")
  list(C = as.numeric(mc$statistic),
       p_mc = mc$p.value,
       n_links = sum(unlist(w$lw$weights)))
}

res_C <- do.call(rbind, lapply(seq_len(nrow(bins)), function(i) {
  b <- bins[i,]
  r <- geary_ring(x, coords, b$d_min, b$d_max, nsim = 499, zero_ok = TRUE)
  data.frame(bin_id = b$bin_id, d_min = b$d_min, d_max = b$d_max,
             C = r$C, p_mc = r$p_mc, n_links = r$n_links)
}))
res_C$mid <- (res_C$d_min + res_C$d_max)/2
res_C
```

### Plot Geary’s C correlogram

```{r, fig.height=4}
ggplot(res_C, aes(mid, C)) +
  geom_hline(yintercept = 1, linetype = 2) +
  geom_point() + geom_line() +
  labs(x = "Distance (bin midpoint, m)", y = "Geary's C",
       title = paste0("Geary's C by distance (", value_col, ")")) +
  theme_minimal()
```

# 7) Semivariogram (gstat)

The **semivariogram** uses **half the average squared difference** between values at a given lag distance. It is **variance-like** (non-negative, in squared units) and usually increases with distance toward a **sill**.

```{r, fig.height=4}
# gstat works with sp classes; convert sf to Spatial
pts_sp <- as(pts, "Spatial")
vg <- gstat::variogram(as.formula(paste0(value_col, " ~ 1")), data = pts_sp)
plot(vg, main = paste0("Semivariogram (", value_col, ")"))
```

# 8) Notes & tips

- **Max lag choice:** using 1/3 of the maximum pairwise distance is a rule of thumb; increase/decrease based on domain knowledge and sampling density.
- **Bin count:** with more points, you can use more bins; with few points, use fewer bins to avoid empty rings.
- **Moran vs Geary vs Variogram:** Moran’s I is correlation-like (can go negative → dispersion). Geary’s C and semivariograms are dissimilarity/variance-like (non-negative; Geary’s C centered at 1).
- **CRS:** always project to meters before distance-based analyses.

# 9) Try the other dataset

Change `file_name` at the top to:
- `"dummy_no_spatial.csv"` to use **KuruPrev** (little/no spatial dependence), or
- back to `"dummy_strong_spatial.csv"` for **COVIDincid** (clear trend).

Then re-knit and compare the correlograms and semivariograms.