spatial_diffusion.R
In dsem: Dynamic Structural Equation Models

## ----include = FALSE, warning=FALSE, message=FALSE----------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
# Install locally
#  devtools::install_local( R'(C:\Users\James.Thorson\Desktop\Git\dsem)', force=TRUE )
# Build
#  setwd(R'(C:\Users\James.Thorson\Desktop\Git\dsem)'); devtools::build_rmd("vignettes/spatial_diffusion.Rmd")

## ----simulate, echo=TRUE, message=FALSE---------------------------------------
library(igraph)
library(Matrix)

# Simulation
adjacency_graph = make_graph( ~ A - B - C - D - E )
A = as.matrix( adjacency_graph )

# Diffusion rate
Dprime = 1 * A
diag(Dprime) = -1 * rowSums(Dprime)

# Movement transition matrix
M = expm( Dprime )

# set seed for reproducibility
set.seed(101)

# Simulate densities
n_times = 100
n_burnin = 100
x_ti = matrix( NA, nrow=n_times+n_burnin, ncol = nrow(M) )
x_ti[1,] = rnorm(n=nrow(M), mean = 0, sd = 1 )
for( t in 2:nrow(x_ti) ){
  x_ti[t,] = (x_ti[t-1,] %*% M)[1,] + rnorm(n=nrow(M), mean = 0, sd = 0.1)
}

# Subset to times after burn-in
x_ti = x_ti[ n_burnin+seq_len(n_times), ]

## ----fit, echo=TRUE, message=FALSE--------------------------------------------
library(dsem)

# Specify SEM
sem = "
  # Spatial correlation
  A -> B, 0, d0
  B -> C, 0, d0
  C -> D, 0, d0
  D -> E, 0, d0
  E -> D, 0, d0
  D -> C, 0, d0
  C -> B, 0, d0
  B -> A, 0, d0

  # Spatio-temporal diffusion
  A -> B, 1, d
  B -> C, 1, d
  C -> D, 1, d
  D -> E, 1, d
  E -> D, 1, d
  D -> C, 1, d
  C -> B, 1, d
  B -> A, 1, d

  # Self-limitation
  A -> A, 1, rho
  B -> B, 1, rho
  C -> C, 1, rho
  D -> D, 1, rho
  E -> E, 1, rho
"

# Fit
colnames(x_ti) = c("A","B","C","D","E")
fit = dsem(
  tsdata = ts(x_ti),
  sem = sem
)

## ----predict, echo=TRUE, message=FALSE----------------------------------------
# Calculate total effect
effect = total_effect(fit, n_lags = 3)

# Calculate predicted movement
Mhat = array( subset(effect,lag==1)$total_effect, dim(M) )
dimnames(Mhat) = dimnames(M)

# Display predicted movement
knitr::kable( Mhat, digits=2)

## ----display, echo=TRUE, message=FALSE----------------------------------------
knitr::kable( as.matrix(M), digits=2)