---
title: "Capture the Dominant Spatial Pattern with Two-Dimensional Locations"
author: "Wen-Ting Wang"
output:
    rmarkdown::html_vignette:
    fig_width: 6
    fig_height: 4
vignette: >
  %\VignetteIndexEntry{Capture the Dominant Spatial Pattern with Two-Dimensional Locations}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  dpi = 300
)
```

## Objective
Represent how to use **SpatPCA** for two-dimensional data for capturing the most dominant spatial pattern

## Basic settings
#### Used packages
```{r message=FALSE}
library(SpatPCA)
library(ggplot2)

base_theme <- theme_minimal(base_size = 10, base_family = "Times") +
  theme(legend.position = "bottom")
fill_bar <- guides(fill = guide_colourbar(
    barwidth = 10,
    barheight = 0.5,
    label.position = "bottom")
  )
coltab <- with(
  list(),
  colorRampPalette(c("#3b4cc0", "#f7f7f7", "#b40426"))(128)
)
color_scale_limit <- c(-.8, .8)
```
Selected realizations are displayed as static images below.
#### True spatial pattern (eigenfunction)
- The underlying spatial pattern below indicates realizations will vary dramatically at the center and be almost unchanged at the both ends of the curve.
```{r, out.width = '100%'}
set.seed(1024)
p <- 25
n <- 8
location <-
  matrix(rep(seq(-5, 5, length = p), 2), nrow = p, ncol = 2)
expanded_location <- expand.grid(location[, 1], location[, 2])
unnormalized_eigen_fn <-
  as.vector(exp(-location[, 1] ^ 2) %*% t(exp(-location[, 2] ^ 2)))
true_eigen_fn <-
  unnormalized_eigen_fn / norm(t(unnormalized_eigen_fn), "F")

plot_df <- data.frame(
  location_dim1 = expanded_location[, 1],
  location_dim2 = expanded_location[, 2],
  eigenfunction = true_eigen_fn
)

ggplot(plot_df, aes(location_dim1, location_dim2)) +
  geom_tile(aes(fill = eigenfunction)) +
  scale_fill_gradientn(colours = coltab, limits = color_scale_limit) +
  base_theme +
  labs(title = "True Eigenfunction", fill = "") +
  fill_bar

```

## Experiment
#### Generate 2-D realizations 
- We want to generate 100 random sample based on 
  - The spatial signal for the true spatial pattern is distributed normally with $\sigma=20$
  - The noise follows the standard normal distribution.

```{r}
realizations <- rnorm(n = n, sd = 3) %*% t(true_eigen_fn) + matrix(rnorm(n = n * p^2), n, p^2)
```

#### Animate realizations
- We can see simulated central realizations change in a wide range more frequently than the others.
```{r, out.width = '100%'}
realization_df <- data.frame(
  location_dim1 = expanded_location[, 1],
  location_dim2 = expanded_location[, 2],
  value = realizations[1, ]
)

ggplot(realization_df, aes(location_dim1, location_dim2)) +
  geom_tile(aes(fill = value)) +
  scale_fill_gradientn(colours = coltab, limits = c(-10, 10)) +
  base_theme +
  labs(title = "1st realization", fill = "") +
  fill_bar
```

#### Apply `SpatPCA::spatpca`
We add a candidate set of `tau2` to see how **SpatPCA** obtain a localized smooth pattern.
```{r}
tau2 <- c(0, exp(seq(log(10), log(400), length = 10)))
cv <- spatpca(x = expanded_location, Y = realizations, tau2 = tau2)
eigen_est <- cv$eigenfn
```
#### Compare **SpatPCA** with PCA
The following figure shows that **SpatPCA** can find sparser pattern than PCA, which is close to the true pattern.
```{r, out.width = '100%'}
plot_df <- data.frame(
  location_dim1 = expanded_location[, 1],
  location_dim2 = expanded_location[, 2],
  spatpca = eigen_est[, 1],
  pca = svd(realizations)$v[, 1]
)

plot_df_long <- rbind(
  data.frame(location_dim1 = plot_df$location_dim1,
             location_dim2 = plot_df$location_dim2,
             estimate = "spatpca",
             eigenfunction = plot_df$spatpca),
  data.frame(location_dim1 = plot_df$location_dim1,
             location_dim2 = plot_df$location_dim2,
             estimate = "pca",
             eigenfunction = plot_df$pca)
)

ggplot(plot_df_long, aes(location_dim1, location_dim2)) +
  geom_tile(aes(fill = eigenfunction)) +
  scale_fill_gradientn(colours = coltab, limits = color_scale_limit) +
  base_theme +
  facet_wrap(~estimate) +
  labs(fill = "") +
  fill_bar
```