Looking at spatial lags and anomolies to see if they work out like temporal anomolies/fixed effect models

spatial_causality.R

library(spdep)
library(spatialreg)
library(SpatialTools)
library(dplyr)
library(nlme)
library(mgcv)
library(purrr)
library(stringr)
library(ggplot2)
library(broom)
library(broom.mixed)

#setup space
#set.seed(31415)
nc <- 20 #num cells
my_coords <- expand.grid(1:nc, 1:nc) %>%
  rename(x = Var1, y = Var2) %>%
  st_as_sf(coords = c("x", "y"), remove = FALSE)

dist_mat <- dist1(as.matrix(st_coordinates(my_coords)))
wmat <- 1/dist_mat
diag(wmat) <- 0
weights_list <- mat2listw(wmat)

#Covariance function for GP
cov_fun <- function(d, etasq = 1, l = 3){
  etasq*exp(-(d/(2*l))^2)
}


#a helpful function for visualization
make_persp <- function(nc = 20, y = rmvnorm(1, rep(0, nrow(k2)), k2),
                       pal_cols = colorRampPalette( c("brown", "brown", "green", "lightgreen") )){
  z <- matrix(y, ncol=nc)
  
  jet.colors <- pal_cols
  color <- jet.colors(length(y))
  zfacet <- z[-1, -1] + z[-1, -nc] + z[-nc, -1] + z[-nc, -nc]
  facetcol <- cut(zfacet, length(y))
  
  persp(x=1:nc, y=1:nc, z=matrix(y, ncol=nc), col=color[facetcol],
        theta = 10, phi = 25, xlab="", ylab="", zlab="", box=FALSE, axes=FALSE,
        mar=c(0,0,0,0) )
}

#our values for the underlying latent spatial processes 

k2 <- cov_fun(dist_mat, 1, 2)


#put it into data
dat <- my_coords %>%
  mutate(
    latent = rmvnorm(1, rep(0, nrow(k2)), k2),
    cause_1 = rnorm(length(latent), latent),
    cause_2 = rnorm(length(latent), -2*latent),
    latent_uncor = rmvnorm(1, rep(0, nrow(k2)), k2),
    cause_3 = rnorm(length(latent_uncor), latent_uncor),
    response = rnorm(length(latent), cause_1 + cause_2 + cause_3, 2)
    )

#Spatial Neighborhoods for points
#https://mgimond.github.io/Spatial/spatial-autocorrelation-in-r.html#morans-i-as-a-function-of-a-distance-band
S.dist  <-  dnearneigh(dat, 0, 3)  
lw <- nb2listw(S.dist, style="W",zero.policy=T) 


dat <- dat %>%
  mutate(sp_lag_cause_1 = lag.listw(lw, cause_1),
         sp_lag_response = lag.listw(lw, response),
         cause_1_anomoly = cause_1 - sp_lag_cause_1
  )

# #show us!
# make_persp(y = dat$latent)
# make_persp(y = dat$cause_1)
# make_persp(y = dat$cause_2)
# make_persp(y = dat$cause_3)
# make_persp(y = dat$response)


mod_true <- lm(response ~ cause_1 + cause_2 + cause_3, data = dat)
mod_bad <- lm(response ~ cause_1, data = dat)
mod_splag <- lm(response ~ cause_1 + sp_lag_cause_1, data = dat)
mod_anom_splag <- lm(response ~ cause_1_anomoly + sp_lag_cause_1, data = dat)

mod_gls <- gls(response ~ cause_1, 
               correlation = corGaus(form = ~ x + y),
               data = dat)

mod_gam <- gam(response ~ cause_1 +
                 s(x , y, bs = "gp"),
               method = "REML",
               data = dat)

mod_anom_splag_gam <- gam(response ~ cause_1_anomoly + sp_lag_cause_1 +
                            s(x , y, bs = "gp"),
                          method = "REML", 
                          data = dat)


mod_lagrslm <-lagsarlm(response ~ cause_1,
                        data = dat,
                        listw = lw)

mod_list <- list(mod_true = mod_true, mod_bad = mod_bad, mod_splag = mod_splag,
                 mod_anom_splag = mod_anom_splag, mod_anom_splag_gam = mod_anom_splag_gam,
                 mod_gls = mod_gls, mod_gam = mod_gam)#,
#                 mod_lagrslm = mod_lagrslm)
  

map_df(mod_list, tidy, .id = "model", parametric = TRUE) %>%
  filter(str_detect(term, "cause_1")) %>%
  ggplot(aes(x = term, y = estimate, ymin = estimate-2*std.error, ymax = estimate + 2*std.error,
             color = model)) +
  geom_pointrange(position = position_dodge(width = 0.3)) +
  geom_hline(yintercept = 0, lty = 2) +
  coord_flip() +
  theme_bw()


rsq <- function(mod){
    o <- residuals(mod) + fitted(mod)
    r <- residuals(mod)
    data.frame(rsq = summary(lm(r ~ o))$r.squared)
}

map_df(mod_list, rsq, .id = "model") %>%
  ggplot(aes(x = model, y = rsq, fill = model)) +
  geom_col(position = position_dodge(width = 0.3)) +
  coord_flip() +
  theme_bw()