is_by <- vapply(object[["smooth"]], gratia:::is_by_smooth, logical(1L))
if (any(is_by)) {
S <- vapply(strsplit(S, ":"), `[[`, character(1L), 1L)
}
npara <- 0
nsmooth <- length(S)
if (isTRUE(parametric)) {
terms <- parametric_terms(object)
npara <- length(terms)
p <- vector("list", length = npara)
}
g <- l <- vector("list", length = nsmooth)
for (i in unique(S)) {
eS <- evaluate_smooth(object, smooth = i, n = n, unconditional = unconditional,
overall_uncertainty = overall_uncertainty, dist = dist)
l[S == i] <- split(eS, eS[["smooth"]])
}
l <- l[select]
d <- d[select]
g <- g[select]
if (length(g) == 0L) {
message("Unable to draw any of the model terms.")
return(invisible(g))
}
for (i in seq_along(l)) {
g[[i]] <- draw(l[[i]])
}
if (isTRUE(parametric)) {
for (i in seq_along(terms)) {
p[[i]] <- evaluate_parametric_term(object, term = terms[i])
g[[i + length(g)]] <- draw(p[[i]])
}
}
if (isTRUE(identical(scales, "fixed"))) {
wrapper <- function(x) {
range(x[["est"]] + (2 * x[["se"]]), x[["est"]] -
(2 * x[["se"]]))
}
ylims <- range(unlist(lapply(l, wrapper)))
if (isTRUE(parametric)) {
ylims <- range(ylims, unlist(lapply(p, function(x) range(x[["upper"]],
x[["lower"]]))))
}
gg <- seq_along(g)[c(d == 1L, rep(TRUE, npara))]
for (i in gg) {
g[[i]] <- g[[i]] + lims(y = ylims)# + geom_rug(data = mf, aes(x = ))
}
}
g
}
p <- mydraw.gam(region_gam)
# customize plot labels
p[[1]] <- p[[1]] + labs(title = "% Change Mobility", x ="") + geom_rug(italy_25, mapping = aes(x = homer), inherit.aes = F) + theme_bw() + theme(title = element_text(size = 22), axis.text = element_text(size = 16))
p[[2]] <- p[[2]]  + labs(x = "", title = expression(paste("PM"[2.5], " (", mu,"g/m"^3, ")"))) + geom_rug(italy_25[italy_25$pm25 < 170,], mapping = aes(x = pm25), inherit.aes = F) + theme_bw() + theme(title= element_text(size = 22), axis.text = element_text(size = 16))
p[[3]] <- p[[3]] + geom_rug(italy_25, mapping = aes(x = med_temp), inherit.aes = F) + labs(title = expression(paste("Temperature ", ( degree*C))), y = "", x = "")+ theme_bw() + theme(title = element_text(size = 22), axis.text = element_text(size = 16))
p[[4]] <- p[[4]] + geom_rug(italy_25, mapping = aes(x = med_hum), inherit.aes = F) + labs(title = "% Humidity", y = "", x = "") + theme_bw() + theme(title = element_text(size = 22), axis.text = element_text(size = 16))
# plot custom ordering, full GAM
grid.arrange(arrangeGrob(p[[1]],p[[2]], ncol=1, nrow=2), arrangeGrob(p[[3]], p[[4]], ncol = 1, nrow = 2) , layout_matrix = rbind(c(1,2),c(1,2)))
# modify the draw.gam function to allow manual adjustment of each plot label
mydraw.gam <- function (object, parametric = TRUE, select = NULL, scales = c("free",
"fixed"), align = "hv", axis = "lrtb", n = 100, unconditional = FALSE,
overall_uncertainty = TRUE, dist = 0.1, ...)
{
scales <- match.arg(scales)
S <- smooths(object)
select <- gratia:::check_user_select_smooths(smooths = S, select = select)
d <- gratia:::smooth_dim(object)
take <- d <= 2L
select <- select[take]
S <- S[take]
d <- d[take]
is_re <- vapply(object[["smooth"]], gratia:::is_re_smooth, logical(1L))
is_by <- vapply(object[["smooth"]], gratia:::is_by_smooth, logical(1L))
if (any(is_by)) {
S <- vapply(strsplit(S, ":"), `[[`, character(1L), 1L)
}
npara <- 0
nsmooth <- length(S)
if (isTRUE(parametric)) {
terms <- parametric_terms(object)
npara <- length(terms)
p <- vector("list", length = npara)
}
g <- l <- vector("list", length = nsmooth)
for (i in unique(S)) {
eS <- evaluate_smooth(object, smooth = i, n = n, unconditional = unconditional,
overall_uncertainty = overall_uncertainty, dist = dist)
l[S == i] <- split(eS, eS[["smooth"]])
}
l <- l[select]
d <- d[select]
g <- g[select]
if (length(g) == 0L) {
message("Unable to draw any of the model terms.")
return(invisible(g))
}
for (i in seq_along(l)) {
g[[i]] <- draw(l[[i]])
}
if (isTRUE(parametric)) {
for (i in seq_along(terms)) {
p[[i]] <- evaluate_parametric_term(object, term = terms[i])
g[[i + length(g)]] <- draw(p[[i]])
}
}
if (isTRUE(identical(scales, "fixed"))) {
wrapper <- function(x) {
range(x[["est"]] + (2 * x[["se"]]), x[["est"]] -
(2 * x[["se"]]))
}
ylims <- range(unlist(lapply(l, wrapper)))
if (isTRUE(parametric)) {
ylims <- range(ylims, unlist(lapply(p, function(x) range(x[["upper"]],
x[["lower"]]))))
}
gg <- seq_along(g)[c(d == 1L, rep(TRUE, npara))]
for (i in gg) {
g[[i]] <- g[[i]] + lims(y = ylims)# + geom_rug(data = mf, aes(x = ))
}
}
g
}
p <- mydraw.gam(region_gam)
# customize plot labels
p[[1]] <- p[[1]] + labs(title = "% Change Mobility", x ="") + geom_rug(italy_25, mapping = aes(x = homer), inherit.aes = F) + theme_bw() + theme(title = element_text(size = 22), axis.text = element_text(size = 16))
p[[2]] <- p[[2]]  + labs(x = "", title = expression(paste("PM"[2.5], " (", AQI, ")"))) + geom_rug(italy_25[italy_25$pm25 < 170,], mapping = aes(x = pm25), inherit.aes = F) + theme_bw() + theme(title= element_text(size = 22), axis.text = element_text(size = 16))
p[[3]] <- p[[3]] + geom_rug(italy_25, mapping = aes(x = med_temp), inherit.aes = F) + labs(title = expression(paste("Temperature ", ( degree*C))), y = "", x = "")+ theme_bw() + theme(title = element_text(size = 22), axis.text = element_text(size = 16))
p[[4]] <- p[[4]] + geom_rug(italy_25, mapping = aes(x = med_hum), inherit.aes = F) + labs(title = "% Humidity", y = "", x = "") + theme_bw() + theme(title = element_text(size = 22), axis.text = element_text(size = 16))
# plot custom ordering, full GAM
grid.arrange(arrangeGrob(p[[1]],p[[2]], ncol=1, nrow=2), arrangeGrob(p[[3]], p[[4]], ncol = 1, nrow = 2) , layout_matrix = rbind(c(1,2),c(1,2)))
# load packages and data
library(mgcv)
library(gratia)
library(tidyverse)
library(cowplot)
source("ReadAllData1.R")
# make region and province into factors
italy_25$region <- as.factor(italy_25$region)
italy_25$province_name <- as.factor(italy_25$province_name)
#loop through and save gams
province_names <-unique(italy_25$province_name)
for(province in province_names){
assign( paste(province,"_gam", sep = ""), gam(log(rt)~
s(homer, k=5)+s(pm25, k = 5)+s(med_temp, k = 4)+s(med_hum, k = 4), data = italy_25[italy_25$province_name == province,], family = gaussian))}
region_names <- unique(as.factor(italy_25$region))
#loop through and save gams
for(region in region_names){
assign( paste(region,"_gam", sep = ""), gam(log(rt)~
s(homer, k =5)+s(pm25, k = 5)+s(med_temp, k = 4)+s(med_hum, k = 4), data = italy_25[italy_25$region == region,], family = gaussian))}
region_gam <- gam(log(rt) ~s(homer, k = 5)+s(pm25, k = 5)+s(med_temp, k =4) +s(med_hum, k=4), data = italy_25, family = gaussian)
province_plot_list <- list(home = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), pm25 = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), temperature = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), humidity = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)))
province_se_list <- list(home = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), pm25 = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), temperature = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), humidity = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)))
range_list <- list(home = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), pm25 = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), temperature = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), humidity = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)))
#province_plot_list <- list(home = matrix(NA, 160, length(province_names)), pm25 = matrix(NA, 160, length(province_names)), temperature = matrix(NA, 160, length(province_names)), humidity = matrix(NA, 160, length(province_names)))
j = 0
#get only the columns used for gam
italy_small <- italy_25 %>% select(homer, pm25, med_temp, med_hum, province_name)
for(province in province_names){
j = j+1
for(i in 1:4){
temp_gam <- get(paste(province, "_gam", sep = ""))
#province_plot_list[[i]][,j] <- evaluate_smooth(temp_gam, smooth = unique(smooths(region_gam))[i], newdata = italy_small[italy_small$province_name == province,i])$est
#temp_gam <- get(paste(province, "_gam", sep = ""))
province_plot_list[[i]][,j] <- evaluate_smooth(temp_gam, smooth = unique(smooths(region_gam))[i])$est
province_se_list[[i]][,j] <- evaluate_smooth(temp_gam, smooth = unique(smooths(region_gam))[i])$se
range_list[[i]][,j] <- seq(min(italy_small[italy_small$province_name == province,i], na.rm = TRUE), max(italy_small[italy_small$province_name == province,i], na.rm =  TRUE), , 100)
}
}
dfhome <- data.frame(province_plot_list[[1]])
names(dfhome) <- province_names
dfhome$time = 1:100
dfhome <- pivot_longer(dfhome,1:23 , names_to = "Province", values_to = "y_val")
# do the same for the error
dfhomese <- data.frame(province_se_list[[1]])
names(dfhomese) <- province_names
dfhomese$time = 1:100
dfhomese <- pivot_longer(dfhomese,1:23 , names_to = "Province", values_to = "y_se")
homex <- data.frame(range_list[[1]])
names(homex) <-province_names
homex$time = 1:100
homex <- pivot_longer(homex, 1:23, names_to = "Province", values_to = "x_val")
dfhome_all <- full_join(homex, dfhome)
dfhome_all2 <- full_join(dfhome_all, dfhomese)
dfhome_all2 <- left_join(dfhome_all2, italy_25[,c("province_name", "region")], by = c("Province" = "province_name"))
gg_home <- ggplot(dfhome_all2, aes(x_val, y_val, group = Province, color = region))+geom_line()+ labs(title = "% Change Mobility", y = "Effect", x = "")  +theme_bw()+ theme(text = element_text(size = 14), legend.position = "none")
#pm25 plot
dfpm25 <- data.frame(province_plot_list[[2]])
names(dfpm25) <- province_names
dfpm25$time = 1:100
dfpm25<- pivot_longer(dfpm25,1:23 , names_to = "Province", values_to = "y_val")
pm25x <- data.frame(range_list[[2]])
names(pm25x) <-province_names
pm25x$time = 1:100
pm25x <- pivot_longer(pm25x, 1:23, names_to = "Province", values_to = "x_val")
dfpm25_all <- full_join(pm25x, dfpm25)
dfpm25se <- data.frame(province_se_list[[2]])
names(dfpm25se) <- province_names
dfpm25se$time = 1:100
dfpm25se <- pivot_longer(dfpm25se,1:23 , names_to = "Province", values_to = "y_se")
dfpm25_all2 <- full_join(dfpm25_all, dfpm25se)
dfpm25_all2 <- left_join(dfpm25_all2, italy_25[,c("province_name", "region")], by = c("Province" = "province_name"))
gg_pm25 <- ggplot(dfpm25_all2, aes(x_val, y_val, group= Province, color = region))+geom_line() + labs(x = "", title = expression(paste("PM"[2.5], " (", mu,"g/m"^3, ")")), y = "", color = "Region") +theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
# extract a legend that is laid out horizontally
legend_b <- get_legend(
gg_pm25 +
guides(color = guide_legend(nrow = 2)) +
theme(legend.position = "bottom")
)
#plot for temp
dftemp <- data.frame(province_plot_list[[3]])
names(dftemp) <- province_names
dftemp$time = 1:100
dftemp<- pivot_longer(dftemp,1:23 , names_to = "Province", values_to = "y_val")
tempx <- data.frame(range_list[[3]])
names(tempx) <-province_names
tempx$time = 1:100
tempx <- pivot_longer(tempx, 1:23, names_to = "Province", values_to = "x_val")
dftemp_all <- full_join(tempx, dftemp)
dftemp_all <- left_join(dftemp_all, italy_25[,c("province_name", "region")], by = c("Province" = "province_name"))
gg_temp <- ggplot(dftemp_all, aes(x_val, y_val, group = Province, color = region))+geom_line()+labs(x = "", title = expression(paste("Temperature ", ( degree*C))), y = "Effect")  + theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
# plot for humidity
dfhumid <- data.frame(province_plot_list[[4]])
names(dfhumid) <- province_names
dfhumid$time = 1:100
dfhumid<- pivot_longer(dfhumid,1:23 , names_to = "Province", values_to = "y_val")
humidx <- data.frame(range_list[[4]])
names(humidx) <-province_names
humidx$time = 1:100
humidx <- pivot_longer(humidx, 1:23, names_to = "Province", values_to = "x_val")
dfhumid_all <- full_join(humidx, dfhumid)
dfhumid_all <- left_join(dfhumid_all, italy_25[,c("province_name", "region")], by = c("Province" = "province_name"))
gg_humid <- ggplot(dfhumid_all, aes(x_val, y_val, group = Province, color = region))+geom_line()+labs(x = "", title = "% Humidity", y = "") + theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
# plot all 4
p4 <- plot_grid(gg_home, gg_pm25, gg_temp, gg_humid, ncol = 2)
# add legend
plot_grid(p4, legend_b, ncol = 1,rel_heights = c(1,.1))
# set up all plots for grid plotting
dfhome <- data.frame(region_plot_list[[1]])
dfhome <- data.frame(province_plot_list[[1]])
names(dfhome) <- province_names
dfhome$time = 1:100
dfhome <- pivot_longer(dfhome,1:23 , names_to = "Province", values_to = "y_val")
# do the same for the error
dfhomese <- data.frame(province_se_list[[1]])
names(dfhomese) <- province_names
dfhomese$time = 1:100
dfhomese <- pivot_longer(dfhomese,1:23 , names_to = "Province", values_to = "y_se")
homex <- data.frame(range_list[[1]])
names(homex) <-province_names
homex$time = 1:100
homex <- pivot_longer(homex, 1:23, names_to = "Province", values_to = "x_val")
dfhome_all <- full_join(homex, dfhome)
dfhome_all2 <- full_join(dfhome_all, dfhomese)
dfhome_all2 <- left_join(dfhome_all2, italy_25[,c("province_name", "region")], by = c("Province" = "province_name"))
gg_home <- ggplot(dfhome_all2, aes(x_val, y_val, group = Province, color = region))+geom_line()+ labs(title = "% Change Mobility", y = "Effect", x = "")  +theme_bw()+ theme(text = element_text(size = 14), legend.position = "none")
#pm25 plot
dfpm25 <- data.frame(province_plot_list[[2]])
names(dfpm25) <- province_names
dfpm25$time = 1:100
dfpm25<- pivot_longer(dfpm25,1:23 , names_to = "Province", values_to = "y_val")
pm25x <- data.frame(range_list[[2]])
names(pm25x) <-province_names
pm25x$time = 1:100
pm25x <- pivot_longer(pm25x, 1:23, names_to = "Province", values_to = "x_val")
dfpm25_all <- full_join(pm25x, dfpm25)
dfpm25se <- data.frame(province_se_list[[2]])
names(dfpm25se) <- province_names
dfpm25se$time = 1:100
dfpm25se <- pivot_longer(dfpm25se,1:23 , names_to = "Province", values_to = "y_se")
dfpm25_all2 <- full_join(dfpm25_all, dfpm25se)
dfpm25_all2 <- left_join(dfpm25_all2, italy_25[,c("province_name", "region")], by = c("Province" = "province_name"))
gg_pm25 <- ggplot(dfpm25_all2, aes(x_val, y_val, group= Province, color = region))+geom_line() + labs(x = "", title = expression(paste("PM"[2.5], " (", mu,"g/m"^3, ")")), y = "", color = "Region") +theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
# extract a legend that is laid out horizontally
legend_b <- get_legend(
gg_pm25 +
guides(color = guide_legend(nrow = 2)) +
theme(legend.position = "bottom")
)
#plot for temp
dftemp <- data.frame(province_plot_list[[3]])
names(dftemp) <- province_names
dftemp$time = 1:100
dftemp<- pivot_longer(dftemp,1:23 , names_to = "Province", values_to = "y_val")
tempx <- data.frame(range_list[[3]])
names(tempx) <-province_names
tempx$time = 1:100
tempx <- pivot_longer(tempx, 1:23, names_to = "Province", values_to = "x_val")
dftemp_all <- full_join(tempx, dftemp)
dftemp_all <- left_join(dftemp_all, italy_25[,c("province_name", "region")], by = c("Province" = "province_name"))
gg_temp <- ggplot(dftemp_all, aes(x_val, y_val, group = Province, color = region))+geom_line()+labs(x = "", title = expression(paste("Temperature ", ( degree*C))), y = "Effect")  + theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
# plot for humidity
dfhumid <- data.frame(province_plot_list[[4]])
names(dfhumid) <- province_names
dfhumid$time = 1:100
dfhumid<- pivot_longer(dfhumid,1:23 , names_to = "Province", values_to = "y_val")
humidx <- data.frame(range_list[[4]])
names(humidx) <-province_names
humidx$time = 1:100
humidx <- pivot_longer(humidx, 1:23, names_to = "Province", values_to = "x_val")
dfhumid_all <- full_join(humidx, dfhumid)
dfhumid_all <- left_join(dfhumid_all, italy_25[,c("province_name", "region")], by = c("Province" = "province_name"))
gg_humid <- ggplot(dfhumid_all, aes(x_val, y_val, group = Province, color = region))+geom_line()+labs(x = "", title = "% Humidity", y = "") + theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
# plot all 4
p4 <- plot_grid(gg_home, gg_pm25, gg_temp, gg_humid, ncol = 2)
# add legend
plot_grid(p4, legend_b, ncol = 1,rel_heights = c(1,.1))
region_gam <- gam(log(rt) ~s(homer, k = 5)+s(pm25, k = 5)+s(med_temp, k =4) +s(med_hum, k=4), data = italy_25, family = gaussian)
province_plot_list <- list(home = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), pm25 = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), temperature = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), humidity = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)))
province_se_list <- list(home = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), pm25 = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), temperature = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), humidity = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)))
range_list <- list(home = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), pm25 = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), temperature = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)), humidity = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(province_names)))
#province_plot_list <- list(home = matrix(NA, 160, length(province_names)), pm25 = matrix(NA, 160, length(province_names)), temperature = matrix(NA, 160, length(province_names)), humidity = matrix(NA, 160, length(province_names)))
j = 0
#get only the columns used for gam
italy_small <- italy_25 %>% select(homer, pm25, med_temp, med_hum, province_name)
for(province in province_names){
j = j+1
for(i in 1:4){
temp_gam <- get(paste(province, "_gam", sep = ""))
#province_plot_list[[i]][,j] <- evaluate_smooth(temp_gam, smooth = unique(smooths(region_gam))[i], newdata = italy_small[italy_small$province_name == province,i])$est
#temp_gam <- get(paste(province, "_gam", sep = ""))
province_plot_list[[i]][,j] <- evaluate_smooth(temp_gam, smooth = unique(smooths(region_gam))[i])$est
province_se_list[[i]][,j] <- evaluate_smooth(temp_gam, smooth = unique(smooths(region_gam))[i])$se
range_list[[i]][,j] <- seq(min(italy_small[italy_small$province_name == province,i], na.rm = TRUE), max(italy_small[italy_small$province_name == province,i], na.rm =  TRUE), , 100)
}
}
dfhome <- data.frame(province_plot_list[[1]])
names(dfhome) <- province_names
dfhome$time = 1:100
dfhome <- pivot_longer(dfhome,1:23 , names_to = "Province", values_to = "y_val")
# do the same for the error
dfhomese <- data.frame(province_se_list[[1]])
names(dfhomese) <- province_names
dfhomese$time = 1:100
dfhomese <- pivot_longer(dfhomese,1:23 , names_to = "Province", values_to = "y_se")
homex <- data.frame(range_list[[1]])
names(homex) <-province_names
homex$time = 1:100
homex <- pivot_longer(homex, 1:23, names_to = "Province", values_to = "x_val")
dfhome_all <- full_join(homex, dfhome)
dfhome_all2 <- full_join(dfhome_all, dfhomese)
dfhome_all2 <- left_join(dfhome_all2, italy_25[,c("province_name", "region")], by = c("Province" = "province_name"))
gg_home <- ggplot(dfhome_all2, aes(x_val, y_val, group = Province, color = region))+geom_line()+ labs(title = "% Change Mobility", y = "Effect", x = "")  +theme_bw()+ theme(text = element_text(size = 14), legend.position = "none")
#pm25 plot
dfpm25 <- data.frame(province_plot_list[[2]])
names(dfpm25) <- province_names
dfpm25$time = 1:100
dfpm25<- pivot_longer(dfpm25,1:23 , names_to = "Province", values_to = "y_val")
pm25x <- data.frame(range_list[[2]])
names(pm25x) <-province_names
pm25x$time = 1:100
pm25x <- pivot_longer(pm25x, 1:23, names_to = "Province", values_to = "x_val")
dfpm25_all <- full_join(pm25x, dfpm25)
dfpm25se <- data.frame(province_se_list[[2]])
names(dfpm25se) <- province_names
dfpm25se$time = 1:100
dfpm25se <- pivot_longer(dfpm25se,1:23 , names_to = "Province", values_to = "y_se")
dfpm25_all2 <- full_join(dfpm25_all, dfpm25se)
dfpm25_all2 <- left_join(dfpm25_all2, italy_25[,c("province_name", "region")], by = c("Province" = "province_name"))
gg_pm25 <- ggplot(dfpm25_all2, aes(x_val, y_val, group= Province, color = region))+geom_line() + labs(x = "", title = expression(paste("PM"[2.5], " (", AQI, ")")), y = "", color = "Region") +theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
# extract a legend that is laid out horizontally
legend_b <- get_legend(
gg_pm25 +
guides(color = guide_legend(nrow = 2)) +
theme(legend.position = "bottom")
)
#plot for temp
dftemp <- data.frame(province_plot_list[[3]])
names(dftemp) <- province_names
dftemp$time = 1:100
dftemp<- pivot_longer(dftemp,1:23 , names_to = "Province", values_to = "y_val")
tempx <- data.frame(range_list[[3]])
names(tempx) <-province_names
tempx$time = 1:100
tempx <- pivot_longer(tempx, 1:23, names_to = "Province", values_to = "x_val")
dftemp_all <- full_join(tempx, dftemp)
dftemp_all <- left_join(dftemp_all, italy_25[,c("province_name", "region")], by = c("Province" = "province_name"))
gg_temp <- ggplot(dftemp_all, aes(x_val, y_val, group = Province, color = region))+geom_line()+labs(x = "", title = expression(paste("Temperature ", ( degree*C))), y = "Effect")  + theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
# plot for humidity
dfhumid <- data.frame(province_plot_list[[4]])
names(dfhumid) <- province_names
dfhumid$time = 1:100
dfhumid<- pivot_longer(dfhumid,1:23 , names_to = "Province", values_to = "y_val")
humidx <- data.frame(range_list[[4]])
names(humidx) <-province_names
humidx$time = 1:100
humidx <- pivot_longer(humidx, 1:23, names_to = "Province", values_to = "x_val")
dfhumid_all <- full_join(humidx, dfhumid)
dfhumid_all <- left_join(dfhumid_all, italy_25[,c("province_name", "region")], by = c("Province" = "province_name"))
gg_humid <- ggplot(dfhumid_all, aes(x_val, y_val, group = Province, color = region))+geom_line()+labs(x = "", title = "% Humidity", y = "") + theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
# plot all 4
p4 <- plot_grid(gg_home, gg_pm25, gg_temp, gg_humid, ncol = 2)
# add legend
plot_grid(p4, legend_b, ncol = 1,rel_heights = c(1,.1))
# same analysis by region instead of province
region_plot_list <- list(home = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(region_names)), pm25 = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(region_names)), temperature = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(region_names)), humidity = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(region_names)))
region_se_list <- list(home = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(region_names)), pm25 = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(region_names)), temperature = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(region_names)), humidity = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(region_names)))
range_regionlist <- list(home = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(region_names)), pm25 = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(region_names)), temperature = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(region_names)), humidity = matrix(NA, length(evaluate_smooth(region_gam, smooth = unique(smooths(region_gam))[1])$est), length(region_names)))
j = 0
#get only the columns used for gam
italy_small <- italy_25 %>% select(homer, pm25, med_temp, med_hum, region)
for(region in region_names){
j = j+1
for(i in 1:4){
temp_gam <- get(paste(region, "_gam", sep = ""))
#province_plot_list[[i]][,j] <- evaluate_smooth(temp_gam, smooth = unique(smooths(region_gam))[i], newdata = italy_small[italy_small$region == province,i])$est
#temp_gam <- get(paste(province, "_gam", sep = ""))
region_plot_list[[i]][,j] <- evaluate_smooth(temp_gam, smooth = unique(smooths(region_gam))[i])$est
region_se_list[[i]][,j] <- evaluate_smooth(temp_gam, smooth = unique(smooths(region_gam))[i])$se
range_regionlist[[i]][,j] <- seq(min(italy_small[italy_small$region == region,i], na.rm = TRUE), max(italy_small[italy_small$region == region,i], na.rm =  TRUE), , 100)
}
}
# set up all plots for grid plotting
dfhome <- data.frame(region_plot_list[[1]])
names(dfhome) <- region_names
dfhome$time = 1:100
dfhome<- pivot_longer(dfhome,1:8 , names_to = "region", values_to = "y_val")
homex <- data.frame(range_regionlist[[1]])
names(homex) <-region_names
homex$time = 1:100
homex <- pivot_longer(homex, 1:8, names_to = "region", values_to = "x_val")
dfhome_all <- full_join(homex, dfhome)
dfhomese <- data.frame(region_se_list[[1]])
names(dfhomese) <- region_names
dfhomese$time = 1:100
dfhomese <- pivot_longer(dfhomese,1:8 , names_to = "region", values_to = "y_se")
dfhome_all2 <- full_join(dfhome_all, dfhomese)
gg_home <- ggplot(dfhome_all2, aes(x_val, y_val, color = region))+geom_line(alpha = 0.7) + geom_ribbon(aes(ymin = y_val - y_se, ymax = y_val + y_se, fill = region), alpha = 0.25) + labs(x = "", title = "% Change Mobility", y = "Effect")  + theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
dfpm25 <- data.frame(region_plot_list[[2]])
names(dfpm25) <- region_names
dfpm25$time = 1:100
dfpm25<- pivot_longer(dfpm25,1:8 , names_to = "region", values_to = "y_val")
pm25x <- data.frame(range_regionlist[[2]])
names(pm25x) <-region_names
pm25x$time = 1:100
pm25x <- pivot_longer(pm25x, 1:8, names_to = "region", values_to = "x_val")
dfpm25_all <- full_join(pm25x, dfpm25)
dfpm25se <- data.frame(region_se_list[[2]])
names(dfpm25se) <- region_names
dfpm25se$time = 1:100
dfpm25se <- pivot_longer(dfpm25se,1:8 , names_to = "region", values_to = "y_se")
dfpm25_all2 <- full_join(dfpm25_all, dfpm25se)
gg_pm25 <- ggplot(dfpm25_all2, aes(x_val, y_val, color = region))+geom_line(alpha = 0.7) + geom_ribbon(aes(ymin = y_val - y_se, ymax = y_val + y_se, fill = region), alpha = 0.25) + labs(x = "", title = expression(paste("PM"[2.5], " (", AQI, ")")), y = "", color = "Region", fill = "Region")  + theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
dftemp <- data.frame(region_plot_list[[3]])
names(dftemp) <- region_names
dftemp$time = 1:100
dftemp<- pivot_longer(dftemp,1:8 , names_to = "region", values_to = "y_val")
tempx <- data.frame(range_regionlist[[3]])
names(tempx) <-region_names
tempx$time = 1:100
tempx <- pivot_longer(tempx, 1:8, names_to = "region", values_to = "x_val")
dftemp_all <- full_join(tempx, dftemp)
dftempse <- data.frame(region_se_list[[3]])
names(dftempse) <- region_names
dftempse$time = 1:100
dftempse <- pivot_longer(dftempse,1:8 , names_to = "region", values_to = "y_se")
dftemp_all2 <- full_join(dftemp_all, dftempse)
gg_temp <- ggplot(dftemp_all2, aes(x_val, y_val, color = region))+geom_line(alpha = 0.7) + geom_ribbon(aes(ymin = y_val - y_se, ymax = y_val + y_se, fill = region), alpha = 0.25) + labs(x = "", title = expression(paste("Temperature ", ( degree*C))), y = "Effect") + theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
dfhum <- data.frame(region_plot_list[[4]])
names(dfhum) <- region_names
dfhum$time = 1:100
dfhum<- pivot_longer(dfhum,1:8 , names_to = "region", values_to = "y_val")
humx <- data.frame(range_regionlist[[4]])
names(humx) <-region_names
humx$time = 1:100
humx <- pivot_longer(humx, 1:8, names_to = "region", values_to = "x_val")
dfhum_all <- full_join(humx, dfhum)
dfhumse <- data.frame(region_se_list[[4]])
names(dfhumse) <- region_names
dfhumse$time = 1:100
dfhumse <- pivot_longer(dfhumse,1:8 , names_to = "region", values_to = "y_se")
dfhum_all2 <- full_join(dfhum_all, dfhumse)
gg_hum <- ggplot(dfhum_all2, aes(x_val, y_val, color = region))+geom_line(alpha = 0.7) + geom_ribbon(aes(ymin = y_val - y_se, ymax = y_val + y_se, fill = region), alpha = 0.25) + labs(x = "", title = "% Humidity", y = "") + theme_bw() + theme(text = element_text(size = 14), legend.position = "none")
# plot all 4
p4_region <- plot_grid(gg_home, gg_pm25, gg_temp, gg_hum, ncol = 2)
# extract a legend that is laid out horizontally
legend_b_region <- get_legend(
gg_pm25 +
guides(color = guide_legend(nrow = 2), fill = guide_legend(nrow = 2)) +
theme(legend.position = "bottom")
)
# add legend
plot_grid(p4_region, legend_b_region, ncol = 1,rel_heights = c(1,.1))
source("ReadAllData1.R")
library(mgcv)
library(ggplot2)
#library(tidyverse)
province_names <- unique(italy_25$province_name)
for (province in province_names) {
new_data = subset(italy_25, province_name!=province)
assign( paste(province,"_exclude_gam", sep = ""), gam(log(rt)~
s(home, k = 5)+s(pm25, k = 5)+s(med_temp, k = 4)+s(med_hum, k = 4) , data = new_data, family = gaussian))
plot(get(paste(province, "_exclude_gam", sep = "")),scale=0,se=2, shade=TRUE,pages=1, scheme = 2, main = province)
}
province_data = subset(italy_25, province_name=="Roma")
i = 0
pred_mat <- matrix(NA, length(province_data$rt), length(province_names))
true_mat <- matrix(NA, length(province_data$rt), length(province_names))
pred_cases <- rep(NA, length(province_data$rt))
mse <- rep(NA, length(province_data$rt))
for (province in province_names) {
i = i + 1
province_data = subset(italy_25, province_name==province)
province_gam <- get(paste(province, "_exclude_gam", sep = ""))
pred_cases <- predict(province_gam,province_data, se.fit = TRUE)
mse[i] <- sum((exp(pred_cases$fit) - province_data$rt)^2, na.rm = T)/length(province_data$rt)
#mse2[i] <- pred_cases$se.fit
#plot(y=exp(pred_cases$fit), x=province_data$rt, pch = 16, xlab = "True", ylab = "Predicted",main = province)
df <- data.frame("Predicted" = exp(pred_cases$fit), "True" = province_data$rt, "Date" =province_data$date)
#df_long <- pivot_longer(df, 1:2, names_to = "Method", values_to = "R_e")
plot_province <- suppressWarnings( ggplot(df) + geom_point(aes(x = True, y = Predicted))+ geom_abline(slope =1, intercept = 0, color = "red") + theme_classic()+ggtitle(label = province))
print(plot_province)
}
