## ----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())