Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Looking at spatial lags and anomolies to see if they work out like temporal anomolies/fixed effect models
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.