## ----include = FALSE--------------------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
options(width = 90)
)
## ----eval = F---------------------------------------------------------------------------
# install.packages("stgam", dependencies = TRUE)
# remotes::install_github("lexcomber/stgam")
## ----warning=F, message=F---------------------------------------------------------------
# load the packages
library(stgam)
library(cols4all) # for nice shading in graphs and maps
library(cowplot) # for managing plots
library(dplyr) # for data manipulation
library(ggplot2) # for plotting and mapping
library(glue) # for model construction
library(mgcv) # for GAMs
library(sf) # for spatial data
library(doParallel) # for parallelising operations
library(purrr) # for model construction
library(tidyr) # for model construction
# load the data
data(productivity)
data(us_data)
## ----locationplot, message = F, warning=F, fig.height = 4, fig.width = 7, fig.cap = "The US States and the geoemtric centroids used as locations."----
ggplot() + geom_sf(data = us_data, fill = NA) +
geom_point(data = productivity |> filter(year == "1970"), aes(x = X, y = Y)) +
theme_bw() + ylab("") + xlab("")
## ---------------------------------------------------------------------------------------
head(productivity)
## ---------------------------------------------------------------------------------------
# define intercept term
productivity <- productivity |> mutate(Intercept = 1)
# create the SVC
svc.gam = gam(privC ~ 0 +
Intercept + s(X, Y, bs = 'gp', by = Intercept) +
unemp + s(X, Y, bs = "gp", by = unemp) +
pubC + s(X, Y, bs = "gp", by = pubC),
data = productivity |> filter(year == "1970"))
## ----ch2gamcheck, fig.height = 7, fig.width = 7, fig.cap = "The GAM GP SVC diagnostic plots."----
# check
gam.check(svc.gam)
# model summary
summary(svc.gam)
## ---------------------------------------------------------------------------------------
get_b2<- productivity |> filter(year == "1970") |> mutate(Intercept = 0, unemp = 0, pubC = 1)
res <- productivity |> filter(year == "1970") |>
mutate(b2 = predict(svc.gam, newdata = get_b2))
## ---------------------------------------------------------------------------------------
get_b0 <- productivity |> filter(year == "1970") |> mutate(Intercept = 1, unemp = 0, pubC = 0)
res <- res |> mutate(b0 = predict(svc.gam, newdata = get_b0))
get_b1 <- productivity |> filter(year == "1970") |> mutate(Intercept = 0, unemp = 1, pubC = 0)
res <- res |> mutate(b1 = predict(svc.gam, newdata = get_b1))
## ----eval = T---------------------------------------------------------------------------
res |> select(b0, b1, b2) |>
apply(2, summary) |> round(2)
## ---------------------------------------------------------------------------------------
terms = c("Intercept", "unemp", "pubC")
res <- calculate_vcs(model = svc.gam,
terms = terms,
input_data = productivity |> filter(year == "1970"))
summary(res[, paste0("b_",terms)])
## ----ch2svccoefs, fig.height = 7, fig.width = 7, fig.cap = "The spatially varying coefficient (SVC) estimates."----
# join the data
map_results <-
us_data |> left_join(res |> select(GEOID, b_Intercept, b_unemp, b_pubC), by = "GEOID")
# plot the insignificant coefficient estimates
tit =expression(paste(""*beta[0]*""))
p1 <-
ggplot(data = map_results, aes(fill=b_Intercept)) +
geom_sf() +
scale_fill_continuous_c4a_div(palette="brewer.spectral",name=tit) +
coord_sf() +
ggtitle("Intercept: not significant")
tit =expression(paste(""*beta[1]*""))
p2 <-
ggplot(data = map_results, aes(fill=b_unemp)) +
geom_sf() +
scale_fill_continuous_c4a_div(palette="brewer.spectral",name=tit) +
coord_sf() +
ggtitle("Unemployment: not significant")
# plot the significant pubC coefficient estimates
tit =expression(paste(""*beta[2]*" "))
p3 <-
ggplot(data = map_results, aes(fill=b_pubC)) +
geom_sf() +
scale_fill_continuous_c4a_div(palette="brewer.prgn",name=tit) +
coord_sf() +
ggtitle("Public captial: significant")
plot_grid(p1, p2, p3, ncol = 1)
## ---------------------------------------------------------------------------------------
# create the TVC
tvc.gam = gam(privC ~ 0 +
Intercept + s(year, bs = 'gp', by = Intercept) +
unemp + s(year, bs = "gp", by = unemp) +
pubC + s(year, bs = "gp", by = pubC),
data = productivity)
## ----eval = F---------------------------------------------------------------------------
# gam.check(tvc.gam)
# summary(tvc.gam)
## ---------------------------------------------------------------------------------------
terms = c("Intercept", "unemp", "pubC")
res <- calculate_vcs(model = tvc.gam,
terms = terms,
input_data = productivity)
summary(res[, paste0("b_",terms)])
## ----ch2tvccoefs, eval = T, echo = T, fig.height = 3, fig.width = 7, message=F, warning=F, fig.cap = "Trends in the temporally varying coefficient estimates."----
res |>
filter(state == "ARIZONA") |>
select(year, b_Intercept, b_unemp, b_pubC) |>
pivot_longer(-year) |>
ggplot(aes(x = year, y = value)) +
geom_point() + geom_line() +
facet_wrap(~name, scale = "free") +
theme_bw() + xlab("Year") + ylab("")
## ---------------------------------------------------------------------------------------
stvc.gam = gam(privC ~ 0 +
Intercept + s(X, Y, year, bs = 'gp', by = Intercept) +
unemp + s(X, Y, year, bs = "gp", by = unemp) +
pubC + s(X, Y, year, bs = "gp", by = pubC),
data = productivity)
## ----eval = F---------------------------------------------------------------------------
# gam.check(stvc.gam)
# summary(stvc.gam)
## ---------------------------------------------------------------------------------------
terms = c("Intercept", "unemp", "pubC")
res <- calculate_vcs(model = stvc.gam,
terms = terms,
input_data = productivity)
summary(res[, paste0("b_",terms)])
## ----eval = T---------------------------------------------------------------------------
res |>
select(year, b_Intercept, b_unemp, b_pubC) |>
group_by(year) |>
summarise(med_b0 = median(b_Intercept),
med_b1 = median(b_unemp),
med_b2 = median(b_pubC))
## ----ch2stvccoefsbox, echo = T, fig.height = 4, fig.width = 7, fig.cap = "The temporal variation of the Public capital coefficient estimates over 17 years.", fig.pos = 'h'----
# inputs to plot
res |> select(starts_with("b"), year) |>
mutate(year = "All Years") -> tmp
cols = c(c4a("tableau.red_gold", n = 17, reverse = T), "grey")
tit =expression(paste(""*beta[`Private Capital`]*""))
# plot
res |> select(starts_with("b"), year) |>
rbind(tmp) |>
mutate(year = factor(year)) |>
ggplot(aes(y = year, x = b_pubC, fill = year)) +
geom_boxplot(outlier.alpha = 0.1) +
scale_fill_manual(values=cols, guide = "none") +
theme_bw() + xlab(tit) + ylab("Time")
## ----ch2stvccoefsmap, message = F, warning = F, fig.height = 4, fig.width = 7, fig.cap = "The spatial variation of the Unemployment coefficient estimates over time."----
tit =expression(paste(""*beta[`Public Capital`]*""))
# join the data
map_results <-
us_data |> left_join(res |> select(GEOID, year, b_Intercept, b_unemp, b_pubC), by = "GEOID")
# create the plot
map_results |>
ggplot() + geom_sf(aes(fill = b_pubC), col = NA) +
scale_fill_binned_c4a_seq(palette="scico.lajolla", name = tit) +
facet_wrap(~year) +
theme_bw() + xlab("") + ylab("") +
theme(
strip.background =element_rect(fill="white"),
strip.text = element_text(size = 8, margin = margin()),
legend.position = "inside", legend.position.inside = c(0.7, 0.1),
legend.direction = "horizontal",
legend.key.width = unit(1, "cm"),
axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())