SMMR revision 1

- Updated the database to correspond to the corrected verion of The Lancet paper
- Further tidied up some sections of the code
- Added a "section 3" of the code for the validation

Author: PierreMasselot

.gitignore

@@ -38,3 +38,4 @@ rsconnect/
 # All data
 temp/
 old/
+figures/

00_DownloadData.R

@@ -12,6 +12,7 @@ library(data.table)
 library(readxl)
 library(zen4R)
 library(dplyr)
+library(tidyr)
 library(giscoR)
 library(sf)
 
@@ -27,7 +28,7 @@ origagebreaks <- c(1, seq(5, 100, by = 5))
 # Target age groups (for the analysis)
 agebreaks <- c(0, 45, 65, 75, 85)
 agelabs <- c(paste(sprintf("%02i", agebreaks[-length(agebreaks)]), 
-  agebreaks[-1], sep = ""), sprintf("%i+", agebreaks[length(agebreaks)]))
+  agebreaks[-1] - 1, sep = ""), sprintf("%i+", agebreaks[length(agebreaks)]))
 
 #------------------------
 # Download mortality data
@@ -40,12 +41,15 @@ linkd <- paste0("https://www.istat.it/storage/dati_mortalita",
   "/decessi-comunali-giornalieri_4-21062023.zip")
 
 # Download zip and open csv
-temp <- tempfile()
+temp <- tempfile(fileext = ".zip")
+tdir <- gsub("\\\\[[:alnum:]]*\\.zip", "", temp)
 download.file(linkd, temp)
 ftr <- grep("csv", zip_list(temp)$filename, value = T)
-deathdata <- fread(cmd = sprintf("unzip -p %s %s", temp, ftr), 
+unzip(temp, ftr, exdir = tdir)
+deathdata <- fread(sprintf("%s\\%s", tdir, ftr), 
   encoding = "Latin-1", na.string = "n.d.")
 unlink(temp)
+unlink(sprintf("%s\\%s", tdir, ftr))
 
 # Pad province codes with zeros
 deathdata[, COD_PROVCOM := sprintf("%06d", COD_PROVCOM)]
@@ -80,13 +84,13 @@ years <- sprintf("20%s", substr(totalvars, 3, 4))
 # Reshape to long
 datalong <- melt.data.table(deathdata, id.vars = c("COD_PROVCOM", "REG", 
     "PROV", "CL_ETA", "GE", "CITY_ID", "CITY_NAME"),
-  measure.vars = totalvars, variable.name = "year", value.name = "all")
+  measure.vars = totalvars, variable.name = "year", value.name = "deaths")
 
 # Rename years
 datalong[, year := gsub("T_", "20", year, fixed = T)]
 
 # Aggregate by city, date and age
-datatab <- datalong[, .(all = sum(all, na.rm = T)), 
+datatab <- datalong[, .(deaths = sum(deaths, na.rm = T)), 
   keyby = .(CITY_ID, year, GE, CL_ETA)]
 
 # Create date variables
@@ -95,31 +99,33 @@ datatab[, day := GE - month * 100]
 datatab[, date := as.Date(paste(year, month, day, sep = "-"))]
 
 # Create full list of dates
-fullfactor <- as.data.table(expand.grid(CITY_CODE = unique(datatab$CITY_ID), 
+fullfactor <- as.data.table(expand.grid(city_code = unique(datatab$CITY_ID), 
   date = seq(dstart, dend, by = "day"), agegroup = 0:21))
 
 # Merge with datatab
 fulldata <- merge(fullfactor, datatab, 
-  by.x = c("CITY_CODE", "date", "agegroup"), 
+  by.x = c("city_code", "date", "agegroup"), 
   by.y = c("CITY_ID", "date", "CL_ETA"), all.x = T, all.y = F)
 
 # Fill NAs
-fulldata[, all := nafill(all, fill = 0)]
+fulldata[, deaths := nafill(deaths, fill = 0)]
 
 # Clean date related variables
 fulldata[, ":="(year = NULL, month = NULL, day = NULL, GE = NULL)]
 
 # Aggregate age-groups
 fulldata[, agegroup := cut(c(0, origagebreaks)[agegroup + 1], c(agebreaks, Inf), 
   right = F, labels = agelabs)]
-fulldata <- fulldata[, .(all = sum(all)), by = .(CITY_CODE, date, agegroup)]
+fulldata <- fulldata[, .(deaths = sum(deaths)), 
+  by = .(city_code, date, agegroup)]
 
 # Transform age groups as wide
-fulldata <- dcast.data.table(fulldata, CITY_CODE + date ~ agegroup, 
-  value.var = "all")
+fulldata <- dcast.data.table(fulldata, city_code + date ~ agegroup, 
+  value.var = "deaths")
+setnames(fulldata, agelabs, sprintf("deaths_%s", agelabs))
 
-# Rename
-setnames(fulldata, agelabs, sprintf("all_%s", agelabs))
+# Arrange
+setorder(fulldata, city_code, date)
 
 #----- Save
 
@@ -137,18 +143,21 @@ if (!file.exists(fname)){
 
 #----- Metadata
 
-# Download metadata used in The Lancet Planetary Health paper from Zenodo
-download_zenodo("10.5281/zenodo.7672108", path = "data",
+# Download data used in The Lancet Planetary Health paper from Zenodo
+# Also include the metadata
+# download_zenodo("10.5281/zenodo.7672108", path = tdir,
+#   files = list("metadata.csv"))
+download_zenodo("10.5281/zenodo.10288665", path = tdir,
   files = list("metadata.csv"))
 
 # Read metadata
-metadata <- read.csv("data/metadata.csv")
-unlink("data/metadata.csv")
+metadata <- read.csv(sprintf("%s/metadata.csv", tdir))
+unlink(sprintf("%s/metadata.csv", tdir))
 
 # Select only Italy
 metadata <- subset(metadata, CNTR_CODE == "IT", -c(CNTR_CODE, cntr_name, region, 
-  nmiss, mcc_code, cityname, country, inmcc)) |>
-  rename(CITY_CODE = "URAU_CODE", CITY_NAME = "URAU_NAME")
+    nmiss, mcc_code, cityname, country, inmcc, LABEL)) |>
+  rename(city_code = "URAU_CODE", city_name = "URAU_NAME")
 
 #----- Info about municipalities
 
@@ -157,12 +166,15 @@ linkg <- paste0("https://www.istat.it/storage/codici-unita-amministrative",
   "/Elenco-codici-statistici-e-denominazioni-delle-unita-territoriali.zip")
 
 # Download and load into session
-temp <- tempfile()
+temp <- tempfile(fileext = ".zip")
+tdir <- gsub("\\\\[[:alnum:]]*\\.zip", "", temp)
 download.file(linkg, temp)
-ftr <- grep("csv", zip_list(temp)$filename, value = T)
-geoinfo <- fread(cmd = sprintf("unzip -p %s %s", temp, ftr), 
+ftr <- grep("csv", zip_list(temp)$filename, value = T, useBytes = T)
+unzip(temp, ftr, exdir = tdir, junkpaths = T)
+geoinfo <- fread(sprintf("%s\\%s", tdir, gsub("^.*/", "", ftr, useBytes = T)), 
   encoding = "Latin-1", na.string = "n.d.")
 unlink(temp)
+unlink(sprintf("%s\\%s", tdir, gsub("^.*/", "", ftr, useBytes = T)))
 
 # Link to cities, select and rename
 geoinfo <- mutate(geoinfo, `LAU CODE` = sprintf("%06d", 
@@ -173,12 +185,29 @@ geoinfo <- mutate(geoinfo, `LAU CODE` = sprintf("%06d",
   unique()
 
 # Add to metadata
-metadata <- merge(metadata, geoinfo, by.x = "CITY_CODE", by.y = "CITY_ID")
+metadata <- merge(metadata, geoinfo, by.x = "city_code", by.y = "CITY_ID")
+
+# Separate spatial and age-related variables
+agevars <- grep("[[:alpha:]]+_[[:digit:]]+$", names(metadata), value = T)
+metadata_spatial <- metadata[, !names(metadata) %in% agevars]
+
+# Age-related variables into long and characterise age groups
+metadata_age <- metadata[, c("city_code", agevars)]
+names(metadata_age) <- gsub("([^_])[[:digit:]]{2}$", "\\1", names(metadata_age))
+metadata_age <- pivot_longer(metadata_age, cols = !city_code, 
+    names_to = c(".value", "age"), names_sep = "_") |>
+  arrange(city_code, age) |>
+  mutate(age = as.numeric(age), agehigh = lead(age) - 1, .by = city_code) |>
+  select(city_code, age, agehigh, prop, deathrate, lifexp)
+
 
 # Write metadata
-fname <- "data/metadata.csv.gz"
+fname <- "data/metadata_spatial.csv.gz"
 if (!file.exists(fname)){
-  fwrite(metadata, fname, quote = F, compress = "gzip")
+  fwrite(metadata_spatial, "data/metadata_spatial.csv.gz", 
+    quote = F, compress = "gzip")
+  fwrite(metadata_age, "data/metadata_age.csv.gz", 
+    quote = F, compress = "gzip")
 }
 
 #----- Geographical data

11_Packages.R

@@ -41,3 +41,4 @@ library(ggdist)
 library(fields)
 library(tidyterra)
 library(ggnewscale) # To have several fill scales on the same plot
+library(ggforce) # For facet paginate

12_Parameters.R

@@ -29,7 +29,7 @@ lagknots <- logknots(maxlag, 3)
 
 # List of metapredictors in composite indices
 metaprednames <- c(Population = "pop", 'Population above 65' = "prop_65p", 
-  'Population density' = "popdens", 'Life expectancy' = "lifexp_00", 
+  'Population density' = "popdens", 
   'Social isolation' = "isol", GPD = "gdp", 'Unemployment rate' = "unempl", 
   'Educational level' = "educ", 'Deprivation rate' = "depriv", 
   'Hospital bed rate' = "bedrates", Imperviousness = "imperv", 
@@ -40,15 +40,15 @@ metaprednames <- c(Population = "pop", 'Population above 65' = "prop_65p",
   'Temperature range' = "trange")
 
 # Number of composite vulnerability indices
-npc <- 5
+npc <- 6
 
 # Number of simulations for eCI
 nsim <- 1000
 
 #----- Results
 
 # Temperature percentile grid for ERFs
-predper <- c(seq(0.1,1,0.1), 2:98, seq(99,99.9,0.1))
+predper <- c(seq(0,1,0.1), 2:98, seq(99,100,0.1))
 
 # Acceptable MMP range
 mmprange <- c(25, 99)
@@ -58,3 +58,6 @@ axisper <- c(1, 25, 50, 75, 99)
 
 # Denominator for death rates
 byrate <- 10^5
+
+# Ages to which to display results
+agebreaks <- c(45, 65, 75, 85)

13_PrepData.R

@@ -15,44 +15,46 @@ mortdata <- fread("data/mortality.csv.gz")
 
 # Read temperature series and merge
 tempdata <- fread("data/tmean.csv.gz")
-fulldata <- merge(mortdata, tempdata, by.x = c("CITY_CODE", "date"), 
-  by.y = c("URAU_CODE", "date"))
+tsdata <- merge(mortdata, tempdata, all.x = T)
 
 # Order
-setkey(fulldata, CITY_CODE, date)
+setkey(tsdata, city_code, date)
 
 # Create date-related variables
-fulldata[, ":="(year = year(date), dow = weekdays(date))]
-
-# Split
-dlist <- split(fulldata[,!"CITY_CODE"], fulldata$CITY_CODE)
-dlist <- lapply(dlist, as.data.frame)
+tsdata[, ":="(year = year(date), dow = weekdays(date))]
 
 #---------------------------
 # Read Metadata
 #---------------------------
 
-# Load data from EUcityTRM
-metadata <- read.csv("data/metadata.csv.gz")
+# Create all city-age combinations
+agelabs <- grep("deaths_[[:digit:]]", names(tsdata), value = T) |> 
+  gsub(pattern = "deaths_", replacement = "")
+metadf <- expand.grid(agegroup = agelabs, 
+  city_code = unique(tsdata$city_code))
+
+# Load data from EUcityTRM and merge
+metadata_spatial <- read.csv("data/metadata_spatial.csv.gz")
+metadf <- merge(metadf, metadata_spatial) 
+
+# Load age-specific demographic data
+metadata_age <- read.csv("data/metadata_age.csv.gz")
+
+
+#---------------------------
+# Separate observed and predicted cities
+#---------------------------
 
 # Draw observed cities for the first stage
 set.seed(1)
-obs <- sample.int(nrow(metadata), nobs)
-metadata$obs <- seq_len(nrow(metadata)) %in% obs
-
-# Create all city-age combinations
-agelabs <- grep("all_[[:digit:]]", names(fulldata), value = T) |> 
-  gsub(pattern = "all_", replacement = "")
-stage2df <- expand.grid(agegroup = agelabs, city = metadata$CITY_CODE)
+obs <- sample(unique(tsdata$city_code), nobs)
 
-# Add name, and geographical information
-stage2df <- merge(stage2df, 
-  metadata[, c("CITY_CODE" ,"CITY_NAME", "geozone", "lon", "lat", "obs")], 
-  by.x = "city", by.y = "CITY_CODE")
+# Separate time series data
+tspred <- tsdata[!city_code %in% obs]
+tsdata <- tsdata[city_code %in% obs]
 
-# Separate data.frame
-stage2_obs <- subset(stage2df, obs)
-stage2_pred <- subset(stage2df, !obs)
+# Include info in metadata
+metadf <- mutate(metadf, obs = city_code %in% obs)
 
 #---------------------------
 # Additional data for results
@@ -64,12 +66,12 @@ italymap <- st_read("data/italymap.shp")
 #----- Common basis to represent curves
 
 # Estimate an overall empirical distribution of temperature
-tmeandist <- t(sapply(dlist, 
-  function(x) quantile(x$tmean, predper / 100, na.rm = T)))
-ovper <- colMeans(tmeandist)
+tmeandist <- tempdata[, .(per = predper, tper = quantile(tmean, predper / 100)), 
+  by = city_code]
+ovper <- tmeandist[, .(tmean = mean(tper)), by = per]$tmean
 
 # Create average basis
-ovknots <- ovper[sprintf("%i.0%%",varper)]
+ovknots <- ovper[predper %in% varper]
 ovbasis <- onebasis(ovper, fun = varfun, degree = vardegree, knots = ovknots)
 
 # Axis locations for plots
@@ -82,13 +84,11 @@ mmtper <- ovper[between(predper, mmprange[1], mmprange[2])]
 #----- Useful objects
 
 # Number of spline coefficients
-nc <- length(varper) + 2
+nc <- ncol(ovbasis)
 
 # Number of cities
-n <- nrow(metadata)
+n <- nrow(metadata_spatial)
 
 # Number of city-ages
-na <- nrow(stage2df)
+na <- nrow(metadf)
 
-# Length of time series
-nts <- nrow(dlist[[1]])