Commit f85acb8e authored by Renaud Lancelot's avatar Renaud Lancelot 🌍
Browse files

Changes in functions to handle vector layers of asministrative data, motivated...

Changes in functions to handle vector layers of asministrative data, motivated by the use of the GAUL data set of administrative country borders. Also, changes in the trend() function to account for special cases, and in the dmrplot function to improve data rendering.
parent cc8414c1
Pipeline #15359 passed with stage
in 4 minutes and 48 seconds
## get lockdown data https://covidtracker.bsg.ox.ac.uk/
## Blavatnik School of Government / Oxford Univ
## Oxford COVID-19 Government Response Tracker
stringency <- function(date_range = "2020-01-01/2020-05-31"){
web <- "https://covidtrackerapi.bsg.ox.ac.uk/"
api <- "api/v2/stringency/"
fjs <- paste(web,
api,
"date-range/",
date_range,
sep = "")
dat <- fromJSON(file = fjs)
## control to check for the existence of data
cond <- paste("length(x$date_value) > 0 &",
"length(x$country_code) > 0 &",
"length(x$deaths) > 0 &",
"length(x$stringency_actual) > 0 &",
"length(x$stringency) > 0 &",
"length(x$stringency_legacy) > 0 &",
"length(x$stringency_legacy_disp) > 0")
## get data when control os TRUE
Liste <- lapply(dat$data,
function(x)
do.call(
what = "rbind",
args = lapply(x,
function(x){
if(eval(parse(text = cond)))
data.frame(x)
else
NULL
})))
## stack
stringency <- do.call("rbind", Liste)
stringency$date <- as.Date(
substr(rownames(stringency), 1, 10))
rownames(stringency) <- seq(nrow(stringency))
## add country name
nam <- sapply(dat$countries,
function(x)
countryname(x))
code <- data.frame(name = standardise_country_names(nam),
country_code = dat$countries)
resu <- merge(stringency, code, by = "country_code")
return(resu[order(resu$name,
resu$date), ])
}
## get country borders for the whole world and get rid of the
## apostrophe in "Cote d'Ivoire"
get_world_map <- function(){
ans <- readOGR(dsn = "./data",
layer = "GAULCOUNTRY2013")
xname <- as.character(ans$ADM0_NAME)
pos <- grep(pattern = "Ivoire", x = xname)
xname[pos] <- "Cote dIvoire"
ans$name <- standardise_country_names(xname)
return(ans)
}
#' Simplify a vector layer #' Simplify a vector layer
geosimplify <- function(x) { geosimplify <- function(x) {
sx <- gSimplify(x, tol=.025, topologyPreserve=T) sx <- gSimplify(x, tol=.025, topologyPreserve=T)
...@@ -13,11 +73,10 @@ region_of_interest <- function(x, context) { ...@@ -13,11 +73,10 @@ region_of_interest <- function(x, context) {
e <- extend(extent(x), 1) e <- extend(extent(x), 1)
E <- as(e, "SpatialPolygons") E <- as(e, "SpatialPolygons")
projection(E) <- projection(x) projection(E) <- projection(x)
return(raster::intersect(context, E)) return(geosimplify(raster::intersect(context, E)))
} }
#' ECDC Covid19 data worldwide #' ECDC Covid19 data worldwide
#' #'
#' Returns the published file at a given date. #' Returns the published file at a given date.
...@@ -56,7 +115,7 @@ ecdc_covid19 <- function(date) { ...@@ -56,7 +115,7 @@ ecdc_covid19 <- function(date) {
## Cleanup variable and country names ## Cleanup variable and country names
names(ans) <- c("date", "day", "month", "year", "cases", "deaths", names(ans) <- c("date", "day", "month", "year", "cases", "deaths",
"country", "geoid", "code", "pop", "continent") "country", "geoid", "code", "pop", "continent", "XXX")
ans$country <- gsub(pattern = "ô", ans$country <- gsub(pattern = "ô",
replacement = "o", replacement = "o",
x = as.character(ans$country)) x = as.character(ans$country))
...@@ -87,31 +146,35 @@ standardise_country_names <- function(x) { ...@@ -87,31 +146,35 @@ standardise_country_names <- function(x) {
## Table of country name variants and their standardised name. ## Table of country name variants and their standardised name.
## Whenever a new country name causes conflict, update this table. ## Whenever a new country name causes conflict, update this table.
country_names_table <- country_names_table <-
tribble( tibble::tribble(
~variant, ~standard_name, ~variant, ~standard_name,
"Syrian Arab Republic", "Syria", "Syrian Arab Republic", "Syria",
"State of Palestine", "Occupied Palestinian Territory", "State of Palestine", "Occupied Palestinian Territory",
"Czechia", "Czech Republic", "Czechia", "Czech Republic",
"Republic of Moldova", "Moldova", "Republic of Moldova", "Moldova",
"Democratic Republic of the Congo", "Congo DR", "Democratic Republic of the Congo", "Congo DR",
"Congo, DRC", "Congo DR", "Congo, DRC", "Congo DR",
"Congo DRC", "Congo DR", "Congo DRC", "Congo DR",
"Congo, DR", "Congo DR", "Congo, DR", "Congo DR",
"Republic of Congo", "Congo", "Congo - Kinshasa", "Congo DR",
"Gambia", "The Gambia", "Republic of Congo", "Congo",
"Cote dIvoire", "Cote d'Ivoire", "Congo - Brazzaville", "Congo",
"Côte d'Ivoire", "Cote d'Ivoire", "Gambia", "The Gambia",
"Côte dIvoire", "Cote d'Ivoire", "Cote d'Ivoire", "Cote dIvoire",
"Cabo Verde", "Cape Verde", "Côte dIvoire", "Cote dIvoire",
"Guinea-Bissau", "Guinea Bissau", "Côte d'Ivoire", "Cote dIvoire",
"Guinea Bissau", "Guinea Bissau", "Côte d'Ivoire", "Cote dIvoire",
"Macedonia", "North Macedonia" "Cabo Verde", "Cape Verde",
) "Guinea-Bissau", "Guinea Bissau",
"Macedonia", "North Macedonia",
"U.K. of Great Britain and Northern Ireland", "United Kingdom",
"United kingdom", "united Kingdom",
"São Tomé and Príncipe", "Sao Tome and Príncipe",
"Sao Tome & Príncipe", "Sao Tome and Príncipe",
"São Tomé & Príncipe", "Sao Tome and Príncipe")
std_name <- function(x) { std_name <- function(x) {
if ( !is.na(idx <- match(x, country_names_table$variant)) ) { if ( !is.na(idx <- match(x, country_names_table$variant)) ) {
x <- country_names_table$standard_name[idx] x <- country_names_table$standard_name[idx]
x <- iconv(x, from="UTF-8", to="ASCII//TRANSLIT")
} }
return(x) return(x)
} }
...@@ -121,6 +184,7 @@ standardise_country_names <- function(x) { ...@@ -121,6 +184,7 @@ standardise_country_names <- function(x) {
) )
} }
#' Get ISO3 codes for a list of countries #' Get ISO3 codes for a list of countries
get_iso3 <- function(x) { get_iso3 <- function(x) {
codes <- ccodes() codes <- ccodes()
...@@ -129,6 +193,17 @@ get_iso3 <- function(x) { ...@@ -129,6 +193,17 @@ get_iso3 <- function(x) {
} }
## get countrt names from the code (to manage the case of Kosovo=)
countryname <- function(x){
ifelse(x == "RKS",
"Kosovo",
countrycode::countrycode(
sourcevar = x,
origin = "iso3c",
destination = "country.name"))
}
#' Total population by country for a given year #' Total population by country for a given year
#' #'
#' Population counts in thousands of people. #' Population counts in thousands of people.
...@@ -236,13 +311,13 @@ dmrplot <- function(dmr, start = as.Date("2020-03-01"), ...@@ -236,13 +311,13 @@ dmrplot <- function(dmr, start = as.Date("2020-03-01"),
idv <- c("name", "date", "lockdown", "d2peak") idv <- c("name", "date", "lockdown", "d2peak")
Z1 <- reshape2::melt(data = z, Z1 <- reshape2::melt(data = z,
id.vars = idv, id.vars = idv,
measure.vars = c("obsdr", "cumdr")) measure.vars = c("fitdr", "cumdr"))
Z2 <- reshape2::melt(data = z, Z2 <- reshape2::melt(data = z,
id.vars = idv, id.vars = idv,
measure.vars = "fitdr") measure.vars = "obsdr")
## incubation 5 d ## incubation 5 d
## early deat 6 d after onset, late eath 46 d after ## early death 6 d after onset, late death 46 d after
## median death 14 days onset ## median death 14 days after onset
## left side at lockdown + 5 + 6, right side at lockdown + 5 + 46 ## left side at lockdown + 5 + 6, right side at lockdown + 5 + 46
## center at morkdown + 5 + (46 - 6) / 2 ## center at morkdown + 5 + (46 - 6) / 2
...@@ -279,15 +354,15 @@ dmrplot <- function(dmr, start = as.Date("2020-03-01"), ...@@ -279,15 +354,15 @@ dmrplot <- function(dmr, start = as.Date("2020-03-01"),
panel = function(x, y, subscripts, ...){ panel = function(x, y, subscripts, ...){
xnam <- unique(Z1$name[subscripts]) xnam <- unique(Z1$name[subscripts])
xvar <- unique(Z1$variable[subscripts]) xvar <- unique(Z1$variable[subscripts])
if(xvar == "obsdr"){ if(xvar == "fitdr"){
## Add daily mortality rate ## add observed daily death rate
lpoints(x, y,
type = "h", col = grey(.8))
## add trend line
X <- Z2$date[Z2$name == xnam] X <- Z2$date[Z2$name == xnam]
Y <- Z2$value[Z2$name == xnam] Y <- Z2$value[Z2$name == xnam]
o <- order(X) o <- order(X)
llines(X[o], Y[o], col=2) lpoints(X[o], Y[o],
type = "h", col = grey(.8))
## Add trend
llines(x, y, col = 2)
} }
if(xvar == "cumdr"){ if(xvar == "cumdr"){
## annotation with mx(umdr) ## annotation with mx(umdr)
...@@ -339,6 +414,7 @@ dmrplot <- function(dmr, start = as.Date("2020-03-01"), ...@@ -339,6 +414,7 @@ dmrplot <- function(dmr, start = as.Date("2020-03-01"),
## fn to compute the trend in mortality growth rate ## fn to compute the trend in mortality growth rate
trend <- function(x){ trend <- function(x){
browser()
## use function "by" to loop over the country names ## use function "by" to loop over the country names
xnam <- unique(x$name) xnam <- unique(x$name)
xlock <- if(is.null(x$lockdown)) NA else unique(x$lockdown) xlock <- if(is.null(x$lockdown)) NA else unique(x$lockdown)
...@@ -369,10 +445,10 @@ trend <- function(x){ ...@@ -369,10 +445,10 @@ trend <- function(x){
k <- length(seq(min(x$date), Sys.Date(), by="1 week")) k <- length(seq(min(x$date), Sys.Date(), by="1 week"))
## fits a negative binomial additive model for death coints ## fits a negative binomial additive model for death coints
fm <- try(gam(deaths ~ s(day, bs="cr", k = k) + offset(log(1000 * pop)), fm <- try(gam(deaths ~ s(day, bs = "cr", k = k) + offset(log(1000 * pop)),
data = x, family = nb, link = "log", data = x, family = nb, link = "log",
optimizer = "perf")) optimizer = "perf"))
if(class(fm) == "try-error") if(any(class(fm) == "try-error"))
return(NULL) return(NULL)
else{ else{
pred <- predict(fm, se = T) pred <- predict(fm, se = T)
...@@ -398,14 +474,11 @@ trend <- function(x){ ...@@ -398,14 +474,11 @@ trend <- function(x){
day <- as.numeric(x$date - min(x$date)) day <- as.numeric(x$date - min(x$date))
## daily mortality growth rate for the early stage of the epidemic ## daily mortality growth rate for the early stage of the epidemic
gr11 <- (dY[day == dlock + 11]) gr11 <- dY[day == dlock + 11]
l0 <- length(gr11)
## relative change in mortality growth rate, reference value = gr11 ## relative change in mortality growth rate, reference value = gr11
rgr19 <- 100 * dY[day = dlock + 19] / gr11 rgr19 <- 100 * dY[day = dlock + 19] / gr11
l19 <- length(rgr19)
rgr46 <- 100 * dY[day = dlock + 46] / gr11 rgr46 <- 100 * dY[day = dlock + 46] / gr11
l46 <- length(rgr46)
## negative growth rate ## negative growth rate
## ldpeak = day when growth rate became < 0 = mortality peak ## ldpeak = day when growth rate became < 0 = mortality peak
...@@ -418,83 +491,87 @@ trend <- function(x){ ...@@ -418,83 +491,87 @@ trend <- function(x){
## for regional consistency ## for regional consistency
## First death to the peak after lockdown ## First death to the peak after lockdown
d2peak <- min(day[day >= dlock & dY < 0], na.rm = T) l1 <- 0
l1 <- length(d2peak) if(!any(day >= dlock & dY < 0))
return(NULL)
## lockdown to peak else{
ldpeak <- d2peak - dlock d2peak <- min(day[day >= dlock & dY < 0], na.rm = T)
l2 <- length(ldpeak) l1 <- length(d2peak)
## max fitted death rate ## lockdown to peak
dr0 <- dmr[day == d2peak] ldpeak <- d2peak - dlock
l3 <- length(dr0)
## max fitted death rate
## decay dr0 <- dmr[day == d2peak]
Y2 <- Y[day >= d2peak]
Ymax <- Y[day == d2peak] ## decay
l4 <- length(Y2) Y2 <- Y[day >= d2peak]
d2 <- day[day >= d2peak] - d2peak Ymax <- Y[day == d2peak]
l5 <- length(d2) d2 <- day[day >= d2peak] - d2peak
## collects the results ## collects the results
fittedTrend <- data.frame( fittedTrend <- data.frame(
name = x$name, name = x$name,
date = x$date, date = x$date,
day = day, day = day,
date1st = rep(min(x$date), nrow(x)), date1st = rep(min(x$date), nrow(x)),
obsdr = 100 * tapply(x$deaths, obsdr = 100 * tapply(x$deaths,
x$day, x$day,
sum) / x$pop, sum) / x$pop,
fitdr = dmr, fitdr = dmr,
cumdr = cumsum(100 * tapply(x$deaths, cumdr = cumsum(100 * tapply(x$deaths,
x$day, x$day,
sum) / x$pop), sum) / x$pop),
pop = unique(x$pop), pop = unique(x$pop),
lockdown = rep(xlock, nrow(x)), lockdown = rep(xlock, nrow(x)),
d2peak = rep(if(l1 > 0) d2peak else NA, nrow(x)), d2peak = rep(if(l1 > 0) d2peak else NA, nrow(x)),
ldpeak = rep(if(l2 > 0) ldpeak else NA, nrow(x)), ldpeak = rep(if(l1 > 0) ldpeak else NA, nrow(x)),
Y = Y, Y = Y,
Ylo = lo, Ylo = lo,
Yhi = hi, Yhi = hi,
dY = dY, dY = dY,
dYlo = dlo, dYlo = dlo,
dYhi = dhi) dYhi = dhi)
## mortality decay ## mortality decay
size <- unique(x$pop) size <- unique(x$pop)
cond10 <- Y2 <= .90 * Ymax cond10 <- Y2 <= .90 * Ymax
cond20 <- Y2 <= .80 * Ymax cond20 <- Y2 <= .80 * Ymax
cond30 <- Y2 <= .70 * Ymax cond30 <- Y2 <= .70 * Ymax
cond40 <- Y2 <= .60 * Ymax cond40 <- Y2 <= .60 * Ymax
cond50 <- Y2 <= .50 * Ymax cond50 <- Y2 <= .50 * Ymax
cond60 <- Y2 <= .40 * Ymax cond60 <- Y2 <= .40 * Ymax
cond70 <- Y2 <= .30 * Ymax cond70 <- Y2 <= .30 * Ymax
cond80 <- Y2 <= .20 * Ymax cond80 <- Y2 <= .20 * Ymax
cond90 <- Y2 <= .10 * Ymax cond90 <- Y2 <= .10 * Ymax
cond95 <- Y2 <= .05 * Ymax cond95 <- Y2 <= .05 * Ymax
fittedParam <- data.frame( issue <- l1 == 0
name = unique(x$name), if(!issue){
lockdown = xlock, fittedParam <- data.frame(
d2peak = if(l1 > 0) d2peak else NA, name = unique(x$name),
ldpeak = if(l2 > 0) ldpeak else NA, lockdown = xlock,
maxdr = if(l3 > 0) dr0 else NA, d2peak = d2peak,
gr11 = if(l0 > 0) gr11 else NA, ldpeak = ldpeak,
rgr19 = if(l19 > 0) rgr19 else NA, maxdr = dr0,
rgr46 = if(l46 > 0) rgr46 else NA, gr11 = gr11,
T10 = if(any(cond10)) min(d2[cond10]) else NA, rgr19 = rgr19,
T20 = if(any(cond20)) min(d2[cond20]) else NA, rgr46 = rgr46,
T30 = if(any(cond30)) min(d2[cond30]) else NA, T10 = if(any(cond10)) min(d2[cond10]) else NA,
T40 = if(any(cond40)) min(d2[cond40]) else NA, T20 = if(any(cond20)) min(d2[cond20]) else NA,
T50 = if(any(cond50)) min(d2[cond50]) else NA, T30 = if(any(cond30)) min(d2[cond30]) else NA,
T60 = if(any(cond70)) min(d2[cond60]) else NA, T40 = if(any(cond40)) min(d2[cond40]) else NA,
T70 = if(any(cond80)) min(d2[cond70]) else NA, T50 = if(any(cond50)) min(d2[cond50]) else NA,
T80 = if(any(cond90)) min(d2[cond80]) else NA, T60 = if(any(cond70)) min(d2[cond60]) else NA,
T90 = if(any(cond90)) min(d2[cond90]) else NA, T70 = if(any(cond80)) min(d2[cond70]) else NA,
T95 = if(any(cond95)) min(d2[cond95]) else NA) T80 = if(any(cond90)) min(d2[cond80]) else NA,
T90 = if(any(cond90)) min(d2[cond90]) else NA,
T95 = if(any(cond95)) min(d2[cond95]) else NA)
return(list(fittedTrend = fittedTrend,
fittedParam = fittedParam))
}
}
} }
return(list(fittedTrend = fittedTrend,
fittedParam = fittedParam))
} }
} }
...@@ -545,10 +622,7 @@ grplot <- function(gr, nc, asp=2/3, cex=.7, dYlim=NA, nam){ ...@@ -545,10 +622,7 @@ grplot <- function(gr, nc, asp=2/3, cex=.7, dYlim=NA, nam){
1/max(x, na.rm = T)) 1/max(x, na.rm = T))
nr <- ceiling(length(nam) / nc) nr <- ceiling(length(nam) / nc)
extrg <- extendrange(range(mydata$dY, na.rm=T)) extrg <- extendrange(range(mydata$dY, na.rm=T))
if(is.na(dYlim)) ylim <- if(any(is.na(dYlim))) extrg else dYlim
ylim <- extrg
else
ylim <- dYlim
xyplot(dY ~ dayL | name, xyplot(dY ~ dayL | name,
data = mydata, subscripts = T, data = mydata, subscripts = T,
ylim = ylim, ylim = ylim,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment