library("gridExtra")
library("ggplot2")
library("dplyr")
library("Rcpp")
library("Gmisc")
library("grid")
library("spdep")
library("sf")
library("spsur")

set.seed(2610)

## ----spcformula------------------------------
spcformula <- WAGE83 | WAGE81 ~ UN83 + NMR83 + SMSA | UN80 + NMR80 + SMSA


## ----data_NCOVR------------------------------
data("NCOVR", package = "spsur")


## ----ncovrlw, warning = FALSE, message = FALSE----
co <- sf::st_coordinates(sf::st_centroid(NCOVR.sf))
ncovrlw <- nb2listw(knn2nb(knearneigh(co, k = 10, longlat = TRUE)), style = "W")


## ----ncovrformula----------------------------
ncovrformula <- HR80 | DV80 | FP79 ~ 
  PS80 + UE80 | PS80 + UE80 + SOUTH | PS80


## ----plotncovr1, warning = FALSE, echo = FALSE, fig.cap = "\\label{Fig-quantile} Spatial distribution of homicide rates in US counties in 1980. (HR80)"----
q <- quantile(NCOVR.sf$HR80)
NCOVR.sf$Quantile<- as.factor((NCOVR.sf$HR80 > q[2]) + (NCOVR.sf$HR80 > q[3]) +
  (NCOVR.sf$HR80 >=  q[4]) + 1)
ggplot(data = NCOVR.sf) +
  geom_sf(aes(fill = Quantile), color = "black", size = .2) +
  theme_bw(base_size = 6) +
  scale_fill_manual(values = c("#FFFEDE", "#FFDFA2", "#FFA93F", "#D5610D"))+
  theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))


## ----LMNCOVR---------------------------------
LMs.ncovr <- lmtestspsur(formula = ncovrformula, data = NCOVR.sf, 
  listw = ncovrlw)
pr.LMs.ncovr <- sapply(LMs.ncovr, print)


## ----spcsurslm-------------------------------
spcsur.slm <- spsurml(formula = spcformula, data = spc, type = "slm", 
  method = "eigen", listw = Wspc)


## ----summary_slm-----------------------------
summary(spcsur.slm)


## ----sdm-------------------------------------
ncovrformulaD <-  ~ PS80 + UE80 | PS80 | PS80 
mlcontrol <- list(fdHess = TRUE, trace = FALSE)
ncovrsur.sdm <- spsurml(formula = ncovrformula, data = NCOVR.sf, 
  type = "sdm", Durbin = ncovrformulaD, listw = ncovrlw,
  method = "LU", control = mlcontrol)
summary(ncovrsur.sdm)


## ----LRtest----------------------------------
spcsur.slm <- spsurml(formula = spcformula, data = spc, type = "slm", 
  listw = Wspc, control = list(trace = FALSE))
spcsur.sdm <- spsurml(formula = spcformula, data = spc, type = "sdm", 
  listw = Wspc, control = list(trace = FALSE))
anova(spcsur.slm, spcsur.sdm)


## ----slm3sls---------------------------------
ncovrsur.slm.3sls <- spsur3sls(formula = ncovrformula, data = NCOVR.sf, 
type = "slm", listw = ncovrlw, trace = FALSE)


## ----print-slm3sls---------------------------
print(ncovrsur.slm.3sls)


## ----waldbeta2-------------------------------
R1 <- matrix(c(1, 0, 0, 0, -1, 0, 0, 0), nrow = 1)
b1 <- matrix(0, ncol = 1)
wald_betas(spcsur.slm, R = R1, b = b1)


## ----estrest---------------------------------
R1 <- matrix(c(1, 0, 0, 0, -1, 0, 0, 0), nrow = 1)
b1 <- matrix(0, ncol = 1)
spcsur.slm.restricted <- spsurml(formula = spcformula, data = spc, 
  type = "slm", listw = Wspc, R = R1, b = b1, control = list(trace = FALSE))
print(spcsur.slm.restricted)


## ----walddeltas------------------------------
R2 <- matrix(c(1, -1), nrow = 1)
b2 <- matrix(0, ncol = 1)
wald_deltas(spcsur.slm.restricted, R = R2, b = b2)


## ----impacts.slm-----------------------------
ncovrw <- as(spdep::listw2mat(ncovrlw), "CsparseMatrix")
ncovrtrw <- spatialreg::trW(ncovrw, type = "MC")
ncovrsur.sdm.impacts <- impactspsur(ncovrsur.sdm, tr = ncovrtrw, R = 1000)


## ----summary_impacts_sdm---------------------
summary(ncovrsur.sdm.impacts[[1]], zstats = TRUE, short = TRUE)


## ----summary_impacts_sdm_2, results = FALSE----
summary(ncovrsur.sdm.impacts[[2]], zstats = TRUE, short = TRUE)


## ----summary_impacts_slx---------------------
ncovrsur.slx <- spsurml(formula = ncovrformula, data = NCOVR.sf, 
  type = "slx", listw = ncovrlw, control = list(trace = FALSE))
ncovrsur.slx.impacts <- impactspsur(ncovrsur.slx)


## ----geometry-hexagon------------------------
sfc <- st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1), c(0, 0)))))
hexs <- st_make_grid(sfc, cellsize = .049, square = FALSE)
hexs.sf <- st_sf(hexs)


## ----DGP-------------------------------------
W <- nb2listw(poly2nb(as(hexs.sf, "Spatial"), queen = FALSE))


## ----dgp-------------------------------------
Sigma <- matrix(c(1, 0.5, 0.5, 1), nrow = 2)
Betas <- c(1, 2, 3, 2, -1, 0.5)
rho <- c(0.2, 0.8)
set.seed(123)
dgp.spatial.sur <- dgp_spsur(Sigma = Sigma, Betas = Betas, Thetas = NULL, 
  lambda = NULL, rho = rho, Tm = 1, G = 2, N = 537, p = 3, 
  listw = W)


## ----sur-sim---------------------------------
sur.dgp.sim <- spsurml(Y = dgp.spatial.sur$Y, X = dgp.spatial.sur$X, 
  type = "sim", G = 2, Tm = 1, N = 537, p = 3, 
  listw = W, control = list(trace = FALSE))
res <- residuals(sur.dgp.sim)


## ----plot-dgp-residuals, warning = FALSE, echo = FALSE, fig.cap = "\\label{Fig-dgp-err} Spatial pattern of the SUR-SIM residuals. $\\epsilon_1$ (left) $\\epsilon_2$ (right)"----
hexs.sf$res.eq1 <- res[[1]]
hexs.sf$res.eq2 <- res[[2]]
p1 <- ggplot(data = hexs.sf) + geom_sf(aes(fill = res.eq1), color = "gray", size = .1) +
  theme_bw(base_size = 6) + theme(legend.position = "none")
p2 <- ggplot(data = hexs.sf) + geom_sf(aes(fill = res.eq2), color = "gray", size = .1) +
  theme_bw(base_size = 6) + theme(legend.position = "none")
gridExtra::grid.arrange(p1, p2, nrow = 1)


## ----sdm_sim---------------------------------
sur.dgp.slm <- spsurml(Y = dgp.spatial.sur$Y, X = dgp.spatial.sur$X, 
  type = "slm", G = 2, Tm = 1, N = 537, p = 3, listw = W,
  control = list(trace = FALSE))
summary(sur.dgp.slm)


## ----boxplot-Mobility-1, warning = FALSE, echo = FALSE, fig.height = 7, fig.width = 10, fig.cap = "\\label{Fig:boxmaps} (A) Weekly intra-province changes in mobility compared to pre-COVID-19 week (February 14-20); (B) Weekly inter-province changes in mobility compared to pre-COVID-19 week (February 14-20). In orange, total lockdown, weeks 6 and 7 (from March 28 to April 10). The red line marks the mobility level = 1 (the reference week)."----
m <- c(0.1, 0.1, 0.1, 0.1)
plot1 <- ggplot(spain.covid, aes(x = as.factor(time), y = Within), fill = color) +
  geom_hline(yintercept = 1, color = "red", size = 1) +
  ylim(0.20, 1.21) +
  ggtitle("(A)") +
  ylab("Intra-Province-Mobility") +
  xlab("Weeks") +
  geom_rect(data = spain.covid, aes(NULL, NULL, xmin = 3.5, xmax = 5.5,
    ymin = -Inf, ymax = +Inf), alpha = 0.01, fill = "gray") +
  geom_rect(data = spain.covid, aes(NULL, NULL, xmin = 5.5, xmax = 7.5,
    ymin = -Inf, ymax = +Inf), alpha = 0.01, fill = "orange") +
  geom_rect(data = spain.covid, aes(NULL, NULL, xmin = 7.5, xmax = 13.5,
    ymin = -Inf, ymax = +Inf), fill = alpha("gray", 0.01)) +
  geom_boxplot() + 
  geom_jitter(color = "black", size = 0.2, alpha = 1) +
  theme_bw() + 
  theme(plot.margin = unit(m, "cm"))

plot2 <- ggplot(spain.covid, aes(x = as.factor(time), y = Exits), fill = color) +
  geom_hline(yintercept = 1, color = "red", size = 1) +
  ylim(0.20, 1.21) +
  ggtitle("(B)") +
  ylab("Inter-Province-Mobility") +
  xlab("Weeks") +
  geom_rect(data = spain.covid, aes(NULL, NULL, xmin = 3.5, xmax = 5.5,
    ymin = -Inf, ymax = +Inf), alpha = 0.01, fill = "gray") +
  geom_rect(data = spain.covid, aes(NULL, NULL, xmin = 5.5, xmax = 7.5,
    ymin = -Inf, ymax = +Inf), alpha = 0.01, fill = "orange") +
  geom_rect(data = spain.covid, aes(NULL, NULL, xmin = 7.5, xmax = 13.5,
    ymin = -Inf, ymax = +Inf), fill = alpha("gray", 0.01)) +
  geom_boxplot() + 
  geom_jitter(color = "black", size = 0.2, alpha = 1) +
  theme_bw() + 
  theme(plot.margin = unit(m, "cm"))

grid.arrange(plot1, plot2, layout_matrix = rbind(c(1, 1, 2, 2), c(NA, 3, 3, NA)),
  ncol = 4, nrow = 2)


## ----boxplot-Mobility-2, warning = FALSE, echo = FALSE, strip.echo = FALSE, fig.height = 6, fig.width = 8, fig.cap = "\\label{Fig:boxmaps-2} (A) Weekly incidence, total positive PCR cases per 100, 000 inhabitants (in logs). In gray, partial lockdown weeks 4-5 and 8-13 (from March 14-27 and April 11th-23). In orange, total lockdown, weeks 6 and 7 (from March 28 to April 10). The blue line marks zero incidence. (B) Quartiles of provincial COVID-19 incidence for week 6 (in log). Similar spatial pattern in oher weeks."----
m <- c(0.1, 0.1, 0.1, 0.1)
plot3 <- ggplot(spain.covid, aes(x = as.factor(time), y = Incidence), fill = color) +
  geom_rect(data = spain.covid, aes(NULL, NULL, xmin = 3.5, xmax = 5.5), 
            ymin = -Inf, ymax = +Inf, alpha = 0.010, fill = alpha("grey", 0.1)) +
  geom_rect(data = spain.covid, aes(NULL, NULL, xmin = 5.5, xmax = 7.5), 
            ymin = -Inf, ymax = +Inf, alpha = 0.010, fill = alpha("orange", 0.1)) +
  geom_rect(data = spain.covid, aes(NULL, NULL, xmin = 7.5, xmax = 13.5), 
            ymin = -Inf, ymax = +Inf, alpha = 0.010, fill = alpha("gray", 0.1)) +
  geom_hline(yintercept = 0, color  = "blue") +
  geom_boxplot() + 
  geom_jitter(color = "black", size = 0.2, alpha = 1) +
  ggtitle("(A)") +
  ylab("Weekly-province-incidence") +
  xlab("Weeks") +
  theme_bw() + 
  theme(plot.margin = unit(m, "cm"))
spain.covid.sf$Incidence <- spain.covid$Incidence[spain.covid$time == 6]
q <- quantile(spain.covid.sf$Incidence)
Incidence.plot <- spain.covid.sf %>%
  mutate(Incidence = Incidence, 
    Incidence = case_when(Incidence < q[2] ~ "1st", 
    Incidence > q[2] & Incidence <=  q[3]  ~ "2nd", 
    Incidence > q[3] & Incidence <=  q[4] ~ "3rd", 
    Incidence > q[4] ~ "4th")) %>%
  ggplot() +
  geom_sf(aes(fill = Incidence), color = "black", size = .2) +
  scale_fill_manual(name = "Quartile", values = c("#FFFEDE", "#FFDFA2", "#FFA93F", "#D5610D")) +
  ggtitle("(B)") +
  theme_bw(base_size = 7) +
  theme(axis.text = element_blank(), 
    legend.position = "bottom",  plot.margin = unit(m, "cm"), 
    plot.title = element_text(size = 13))
grid.arrange(plot3, Incidence.plot, layout_matrix = rbind(c(1, 1, 2, 2), c(NA, 3, 3, NA)),
  ncol = 4, nrow = 2)


## ----data, warning = FALSE-------------------
data("spain.covid", package = "spsur")


## ----formula---------------------------------
formula <- Within | Exits ~ Emergence + EmergenceTotal + Density + 
  Old65 + Essential + Incidence


## ----SUR-SIM---------------------------------
Tm <- 17
covid.sim <- spsurml(formula = formula, data = spain.covid, type = "sim", 
                     Tm = Tm, control = list(trace = FALSE))
summary(covid.sim)


## ----W, echo = 2:3, strip.white = FALSE--------
listw <- spdep::poly2nb(spain.covid.sf, queen = TRUE)
listw <- spdep::nb2listw(listw, style = "W", zero.policy = TRUE)


## ----lm-sim----------------------------------
covid.lmtest <- lmtestspsur(formula = formula, data = spain.covid, 
Tm = Tm, listw = listw, zero.policy = TRUE)
pr.covid.lmtest <- sapply(covid.lmtest, broom::tidy)
print(as.data.frame(t(pr.covid.lmtest)))

## ----all-models, results = FALSE-------------
formulaD <- ~ Essential + Incidence
covid.slx <- spsurml(formula = formula, data = spain.covid, type = "slx", 
  Tm = Tm, Durbin = formulaD, listw = listw, zero.policy = TRUE)
covid.slm <- spsurml(formula = formula, data = spain.covid, type = "slm", 
  Tm = Tm, listw = listw, zero.policy = TRUE)
covid.sem <- spsurml(formula = formula, data = spain.covid, type = "sem", 
  Tm = Tm, listw = listw, zero.policy = TRUE)
covid.sdm <- spsurml(formula = formula, data = spain.covid, type = "sdm", 
  Tm = Tm, Durbin = formulaD, listw = listw, zero.policy = TRUE)
covid.sdem <- spsurml(formula = formula, data = spain.covid, type = "sdem", 
  Tm = Tm, Durbin = formulaD, listw = listw, zero.policy = TRUE)
covid.sarar <- spsurml(formula = formula, data = spain.covid, 
  type = "sarar", Tm = Tm, listw = listw, zero.policy = TRUE)
covid.gnm <- spsurml(formula = formula, data = spain.covid, type = "gnm", 
  Tm = Tm, Durbin = formulaD, listw = listw, zero.policy = TRUE)


## ----logLik----------------------------------
anova(covid.sim, covid.slx, covid.slm, covid.sem, covid.sdm, covid.sdem, 
  covid.sarar, covid.gnm, lrtest = FALSE)


## ----lrtest-3--------------------------------
anova(covid.sdem, covid.gnm, lrtest = TRUE)


## ----SUR-GNM---------------------------------
print(covid.gnm)


## ----Walds1----------------------------------
R <- matrix(0, ncol = Tm + 1, nrow = 1)
R[2] <- 1
R[11] <- -1
wald_betas(covid.gnm, R = R, b = 0)


## ----Walds2----------------------------------
R <- matrix(0, ncol = Tm + 1, nrow = 1)
R[3] <- 1
R[12] <- -1
wald_betas(covid.gnm, R = R, b = 0)


## ----SUR-GNM2--------------------------------
R <- matrix(0, ncol = Tm + 1, nrow = 1)
R[3]  <- 1 
R[12] <- -1
covid.gnm.r <- spsur::spsurml(formula = formula, data = spain.covid, 
  type = "gnm", Tm = Tm, Durbin = formulaD, R = R, b = 0, 
  listw = listw, zero.policy = TRUE, control = list(trace = FALSE))
print(covid.gnm.r)

## ----SUR-GNM3--------------------------------
formula <- Within | Exits ~ 
  Emergence + EmergenceTotal + Density + Essential + Incidence |
  Emergence + EmergenceTotal + Density + Old65 + Essential + Incidence
formulaD <- ~ Essential + Incidence | Essential
R <- matrix(0, ncol = 16, nrow = 1)
R[3]  <- 1 
R[11] <- -1
covid.gnm.rr <- spsurml(formula = formula, data = spain.covid, 
  type = "gnm", Tm = Tm, Durbin = formulaD, R = R, b = 0, listw = listw,
  zero.policy = TRUE, control = list(trace = FALSE))
print(covid.gnm.rr)


## ----wald_delta_1----------------------------
R2 <- matrix(c(1, 0, -1, 0), nrow = 1)
b2 <- matrix(0, ncol = 1)
wald_deltas(covid.gnm.rr, R = R2, b = b2)


## ----wald_delta_2----------------------------
R2 <- matrix(c(0, 1, 0, -1), nrow = 1)
b2 <- matrix(0, ncol = 1)
wald_deltas(covid.gnm.rr, R = R2, b = b2)


## ----impacts---------------------------------
impacts.covid.gnm.rr <- impactspsur(covid.gnm.rr, listw = listw, R = 1000)


## --------------------------------------------
summary(impacts.covid.gnm.rr[[1]], zstats = FALSE, short = TRUE)


## --------------------------------------------
summary(impacts.covid.gnm.rr[[2]], zstats = FALSE, short = TRUE)


## --------------------------------------------
sessionInfo()