Commit 4b9f1c74 authored by Renaud Lancelot's avatar Renaud Lancelot 🌍
Browse files

Major revisions of the code of graphics functions dmrplot and grplot.

parent 865baafa
......@@ -203,107 +203,109 @@ collapse <- function(nam = onam, cnt = nd){
## fn to plot the daily mortality rate
dmrplot <- function(dmr, start = as.Date("2020-03-01"),
end = Sys.Date(), cex, nam, asp = 2/3){
no.ld <- is.null(dmr$lockdown)
if(no.ld)
dmr$lockdown <- NA
nam <- as.character(nam)
## cumulative mortality in countries listed in nam
zzz1 <- base::subset(dmr, as.character(name) %in% nam)
zzz1$name <- as.character(zzz1$name)
## data set with cumulative mortality
L <- by(zzz1,
list(name = zzz1$name),
function(x){
o <- order(x$date)
numo <- 100 * cumsum(x$deaths[o])
deno <- x$pop
data.frame(
name = x$name,
date = x$date[o],
pop = x$pop,
lockdown = if(no.ld) rep(NA, nrow(x)) else x$lockdown,
value = numo / deno)
})
## stack list elements
zzz2 <- do.call("rbind", args = L)
## rates in deaths / 100.000 inh.
zzz1$value <- 100 * with(zzz1, deaths / pop)
## stack
Z <- rbind(
base::subset(zzz1,
select = c("name", "date", "pop", "lockdown", "value")),
base::subset(zzz2,
select = c("name", "date", "pop", "lockdown", "value")))
Z$variable <- factor(c(rep("Daily death rate", nrow(zzz1)),
rep("Cumulative death rate", nrow(zzz2))),
levels = c("Daily death rate",
"Cumulative death rate"))
Z$name <- reorder(Z$name,
Z$value,
FUN = function(x) 1 / max(x))
z <- subset(dmr, name %in% nam)
idv <- c("name", "date", "lockdown", "d2peak")
Z1 <- reshape2::melt(data = z,
id.vars = idv,
measure.vars = c("obsdr", "cumdr"))
Z2 <- reshape2::melt(data = z,
id.vars = idv,
measure.vars = "fitdr")
## incubation 5 d
## early deat 6 d after onset, late eath 46 d after
## median death 14 days onset
## left side at lockdown + 5 + 6, right side at lockdown + 5 + 46
## center at morkdown + 5 + (46 - 6) / 2
z2name <- unique(zzz2$name)
cdr <- tapply(zzz2$value, zzz2$name, function(x) max(x, na.rm=T))
scdr <- sprintf(cdr, fmt = "%##.1f")
t_end <- tapply(zzz1$date, zzz1$name, function(x) max(x, na.rm=T))
t_beg <- tapply(zzz1$date, zzz1$name, function(x) min(x, na.rm=T))
t_fin <- as.numeric(t_end - t_beg)
lbl <- paste0("expression(italic(c)[", t_fin, "] ==", scdr, ")")
## regular xyplot
g <- xyplot(value ~ date | variable * name,
data = Z, subscripts = T,
xlim = c(start, end),
## key
mykey <- list(
space = "right", border = 1,
text = list(c("data",
"trend",
"lockdown",
"post-lockdown\nd 11, 19, 46",
"peak"), cex = cex),
lines = list(lty = c(1, 1, 4, 4, 1),
col = c(grey(.8), 2, 1, 2, 2),
size = 3))
g <- xyplot(value ~ date | variable + name,
xlim = c(as.Date("2020-02-15"), Sys.Date()),
data = Z1, subscripts = T,
scales = list(cex = cex - .1, tck = .5,
x = list(rot = 0,
relation = "same"),
y = list(rot = 0,
relation = "free")),
layout = c(2, length(unique(zzz2$name))),
layout = c(2, length(nam)),
key = mykey,
aspect = asp,
between = list(x = 1/4, y = 1/4),
xlab = character(0),
ylab = list("Deaths / 100,000 inh.",
cex = cex),
par.strip.text = list(cex = cex - .2),
par.strip.text = list(cex = cex - .1),
as.table = T,
panel = function(x, y, subscripts, ...){
panel.grid(h = -1)
w <- unique(Z$variable[subscripts])
n <- unique(Z$name[subscripts])
no.ld <- all(is.na(Z$lockdown[subscripts]))
if(!no.ld){
lddt <- Z$lockdown[subscripts]
ld <- unique(as.numeric(lddt))
panel.abline(v = ld + c(0, 11, 19, 46),
col = c(1, 2, 2, 2),
lty = 4,
subscripts)
xnam <- unique(Z1$name[subscripts])
xvar <- unique(Z1$variable[subscripts])
if(xvar == "obsdr"){
## Add daily mortality rate
lpoints(x, y,
type = "h", col = grey(.8))
## add trend line
X <- Z2$date[Z2$name == xnam]
Y <- Z2$value[Z2$name == xnam]
o <- order(X)
llines(X[o], Y[o], col=2)
}
if(w == "Cumulative death rate"){
lbln <- lbl[which(z2name == n)]
if(xvar == "cumdr"){
## annotation with mx(umdr)
my <- max(y, na.rm = T)
smy <- sprintf(fmt = "%##.1f", my)
hi <- max(x, na.rm = T)
lo <- min(x, na.rm = T)
i <- as.numeric(hi - lo)
lbl <- paste0(
"expression(italic(c[",
i, "]) == ", smy, ")")
grid.text(
x = unit(.5, "mm"),
y = unit(1, "npc") - unit(.5, "mm"),
label = eval(parse(text = lbln)),
label = eval(parse(text = lbl)),
hjust = 0, vjust = 1,
gp = gpar(cex = cex - .1))
## add data (cumdr)
lpoints(x, y,
col = grey(.8), type = "h")
}
## add grid
panel.grid()
## add lockdown date and pivotal dates
lddt <- Z1$lockdown[subscripts]
ld <- unique(as.numeric(lddt))
panel.abline(v = ld + c(0, 11, 19, 46),
col = c(1, 2, 2, 2),
lty = 4,
subscripts)
## add peak date
xd2peak <- unique(Z1$d2peak[subscripts])
xdate1st <- min(Z1$date[subscripts],
na.rm = T)
if(!is.na(xd2peak)){
xdatepeak <- xdate1st + xd2peak
panel.abline(v = xdatepeak,
col=2, lwd=1, lty=1)
}
}) +
latticeExtra::layer(lpoints(x, y, type = "h", col = 4),
columns = 1) +
latticeExtra::layer(lpoints(x, y, type = "l", col = 4),
columns = 2)
})
## plot with outer strips for better use of space
p <- useOuterStrips(g)
## p <- update(p, layout = c(2, XX))
p <- useOuterStrips(
g,
strip = strip.custom(
factor.levels = c("Daily death rate",
"Cumulative daily death rate")))
return(p)
}
## fn to compute the trend in mortality growth rate
......@@ -336,127 +338,135 @@ trend <- function(x){
## initialize data containers
fittedTrend <- fittedParam <- data.frame()
k <- length(seq(min(x$date), Sys.Date(), by="1 week"))
## fits a negative binomial additive model for death coints
fm <- gam(deaths ~ s(day) + offset(log(1000 * pop)),
fm <- try(gam(deaths ~ s(day, bs="cr", k = k) + offset(log(1000 * pop)),
data = x, family = nb, link = "log",
optimizer = "perf")
pred <- predict(fm, se = T)
Y <- exp(pred$fit)
## these are counts we want rates
dmr <- 100 * Y / x$pop
## 95% CI
z <- qnorm(p = .975, lower.tail = T)
lo <- exp(pred$fit - z * pred$se.fit)
hi <- exp(pred$fit + z * pred$se.fit)
## compute the growth rates (1st derivative)
dfm <- derivatives(fm, n = nrow(x))
dY <- dfm$derivative
dhi <- dY + z * dfm$se
dlo <- dY - z * dfm$se
## time from 1st death to lockown
dlock <- as.numeric(xlock - min(x$date))
## time from the first reported death
day <- as.numeric(x$date - min(x$date))
## daily mortality growth rate for the early stage of the epidemic
gr11 <- (dY[day == dlock + 11])
l0 <- length(gr11)
## relative change in mortality growth rate, reference value = gr11
rgr19 <- 100 * dY[day = dlock + 19] / gr11
l19 <- length(rgr19)
rgr46 <- 100 * dY[day = dlock + 46] / gr11
l46 <- length(rgr46)
## negative growth rate
## ldpeak = day when growth rate became < 0 = mortality peak
## ldmax = day when max growth rate was observed
if(any(dY < 0)){
## we look for a peak
## NB for Sweden e have set a lockdown date similar to the
## one adopted by neighbor countries for regional
## consistency
## First death to the peak after lockdown
d2peak <- min(day[day >= dlock & dY < 0], na.rm = T)
l1 <- length(d2peak)
## lockdown to peak
ldpeak <- d2peak - dlock
l2 <- length(ldpeak)
## max fitted death rate
dr0 <- dmr[day == d2peak]
l3 <- length(dr0)
## decay
Y2 <- Y[day >= d2peak]
Ymax <- Y[day == d2peak]
l4 <- length(Y2)
d2 <- day[day >= d2peak] - d2peak
l5 <- length(d2)
## collects the results
fittedTrend <- data.frame(
name = x$name,
date = x$date,
day = day,
date1st = rep(min(x$date), nrow(x)),
obsdr = 100 * tapply(x$deaths,
x$day,
sum) / x$pop,
pop = unique(x$pop),
lockdown = rep(xlock, nrow(x)),
d2peak = rep(if(l1 > 0) d2peak else NA, nrow(x)),
ldpeak = rep(if(l2 > 0) ldpeak else NA, nrow(x)),
Y = Y,
Ylo = lo,
Yhi = hi,
dY = dY,
dYlo = dlo,
dYhi = dhi)
## mortality decay
size <- unique(x$pop)
cond10 <- Y2 <= .90 * Ymax
cond20 <- Y2 <= .80 * Ymax
cond30 <- Y2 <= .70 * Ymax
cond40 <- Y2 <= .60 * Ymax
cond50 <- Y2 <= .50 * Ymax
cond60 <- Y2 <= .40 * Ymax
cond70 <- Y2 <= .30 * Ymax
cond80 <- Y2 <= .20 * Ymax
cond90 <- Y2 <= .10 * Ymax
cond95 <- Y2 <= .05 * Ymax
fittedParam <- data.frame(
name = unique(x$name),
lockdown = xlock,
d2peak = if(l1 > 0) d2peak else NA,
ldpeak = if(l2 > 0) ldpeak else NA,
maxdr = if(l3 > 0) dr0 else NA,
gr11 = if(l0 > 0) gr11 else NA,
rgr19 = if(l19 > 0) rgr19 else NA,
rgr46 = if(l46 > 0) rgr46 else NA,
T10 = if(any(cond10)) min(d2[cond10]) else NA,
T20 = if(any(cond20)) min(d2[cond20]) else NA,
T30 = if(any(cond30)) min(d2[cond30]) else NA,
T40 = if(any(cond40)) min(d2[cond40]) else NA,
T50 = if(any(cond50)) min(d2[cond50]) else NA,
T60 = if(any(cond70)) min(d2[cond60]) else NA,
T70 = if(any(cond80)) min(d2[cond70]) 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)
optimizer = "perf"))
if(class(fm) == "try-error")
return(NULL)
else{
pred <- predict(fm, se = T)
Y <- exp(pred$fit)
## these are counts we want rates
dmr <- 100 * Y / x$pop
## 95% CI
z <- qnorm(p = .975, lower.tail = T)
lo <- exp(pred$fit - z * pred$se.fit)
hi <- exp(pred$fit + z * pred$se.fit)
## compute the growth rates (1st derivative) and 95% CI%
dfm <- derivatives(fm, n = nrow(x))
dY <- dfm$derivative
dhi <- dY + z * dfm$se
dlo <- dY - z * dfm$se
## time from 1st death to lockown
dlock <- as.numeric(xlock - min(x$date))
## time from the first reported death
day <- as.numeric(x$date - min(x$date))
## daily mortality growth rate for the early stage of the epidemic
gr11 <- (dY[day == dlock + 11])
l0 <- length(gr11)
## relative change in mortality growth rate, reference value = gr11
rgr19 <- 100 * dY[day = dlock + 19] / gr11
l19 <- length(rgr19)
rgr46 <- 100 * dY[day = dlock + 46] / gr11
l46 <- length(rgr46)
## negative growth rate
## ldpeak = day when growth rate became < 0 = mortality peak
## ldmax = day when max growth rate was observed
if(any(dY < 0)){
## we look for a peak NB for Sweden - which did not
## implement à strict lockdown, we have set a lockdown
## date similar to the one adopted by neighbor countries
## for regional consistency
## First death to the peak after lockdown
d2peak <- min(day[day >= dlock & dY < 0], na.rm = T)
l1 <- length(d2peak)
## lockdown to peak
ldpeak <- d2peak - dlock
l2 <- length(ldpeak)
## max fitted death rate
dr0 <- dmr[day == d2peak]
l3 <- length(dr0)
## decay
Y2 <- Y[day >= d2peak]
Ymax <- Y[day == d2peak]
l4 <- length(Y2)
d2 <- day[day >= d2peak] - d2peak
l5 <- length(d2)
## collects the results
fittedTrend <- data.frame(
name = x$name,
date = x$date,
day = day,
date1st = rep(min(x$date), nrow(x)),
obsdr = 100 * tapply(x$deaths,
x$day,
sum) / x$pop,
fitdr = dmr,
cumdr = cumsum(100 * tapply(x$deaths,
x$day,
sum) / x$pop),
pop = unique(x$pop),
lockdown = rep(xlock, nrow(x)),
d2peak = rep(if(l1 > 0) d2peak else NA, nrow(x)),
ldpeak = rep(if(l2 > 0) ldpeak else NA, nrow(x)),
Y = Y,
Ylo = lo,
Yhi = hi,
dY = dY,
dYlo = dlo,
dYhi = dhi)
## mortality decay
size <- unique(x$pop)
cond10 <- Y2 <= .90 * Ymax
cond20 <- Y2 <= .80 * Ymax
cond30 <- Y2 <= .70 * Ymax
cond40 <- Y2 <= .60 * Ymax
cond50 <- Y2 <= .50 * Ymax
cond60 <- Y2 <= .40 * Ymax
cond70 <- Y2 <= .30 * Ymax
cond80 <- Y2 <= .20 * Ymax
cond90 <- Y2 <= .10 * Ymax
cond95 <- Y2 <= .05 * Ymax
fittedParam <- data.frame(
name = unique(x$name),
lockdown = xlock,
d2peak = if(l1 > 0) d2peak else NA,
ldpeak = if(l2 > 0) ldpeak else NA,
maxdr = if(l3 > 0) dr0 else NA,
gr11 = if(l0 > 0) gr11 else NA,
rgr19 = if(l19 > 0) rgr19 else NA,
rgr46 = if(l46 > 0) rgr46 else NA,
T10 = if(any(cond10)) min(d2[cond10]) else NA,
T20 = if(any(cond20)) min(d2[cond20]) else NA,
T30 = if(any(cond30)) min(d2[cond30]) else NA,
T40 = if(any(cond40)) min(d2[cond40]) else NA,
T50 = if(any(cond50)) min(d2[cond50]) else NA,
T60 = if(any(cond70)) min(d2[cond60]) else NA,
T70 = if(any(cond80)) min(d2[cond70]) 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))
}
## function to visualize the mortalilty growth rates, as estimated by
......@@ -469,8 +479,9 @@ trend <- function(x){
## plot. Must be a subset of the country names present in the trend()
## output.
grplot <- function(gr, nc, asp=2/3, cex=.7, nam){
mykey <- list()
grplot <- function(gr, nc, asp=2/3, cex=.7, dYlim=NA, nam){
## legend
mykey <- list()
mykey$space <- "right"
mykey$divide <- 1
mykey$rep <- F
......@@ -489,26 +500,47 @@ grplot <- function(gr, nc, asp=2/3, cex=.7, nam){
col = c(4, 4, 1, 4, 2, 4),
lty = c(1, 0, 4, 2, 4, 4),
size = 3, lwd = 1.5)
mydata <- subset(gr, name %in% nam)
## set time = 0 to lockdown - 10
myL <- by(mydata,
list(name = mydata$name),
function(x){
x$dayL <- as.numeric(x$date - x$lockdown)
subset(x, dayL >= 0)
})
mydata <- do.call("rbind", myL)
mydata$name <- reorder(mydata$name,
mydata$cumdr,
function(x)
1/max(x, na.rm = T))
nr <- ceiling(length(nam) / nc)
xyplot(dY ~ day | name, data = mydata, subscripts = T,
extrg <- extendrange(range(mydata$dY, na.rm=T))
if(is.na(dYlim))
ylim <- extrg
else
ylim <- dYlim
xyplot(dY ~ dayL | name,
data = mydata, subscripts = T,
ylim = ylim,
scales = list(cex = cex - .2, tck = .5,
x = list(relation = "same"),
y = list(rot = 0, relation = "same")),
layout = c(nc, nr), between = list(x = 1/4, y = 1/4),
aspect = asp, as.table = T,
y = list(rot = 0,
relation = "same")),
layout = c(nc, nr), as.table = T,
between = list(x = 1/4, y = 1/4),
aspect = asp,
key = mykey,
xlab = list("Days post lockdown\n", cex = cex),
ylab = list("Mortality growth rate", cex = cex),
par.strip.text = list(cex = cex),
panel = function(x, y, subscripts){
xdate1st <- unique(mydata$date1st[subscripts])
xlock <- unique(mydata$lockdown[subscripts])
d2peak <- unique(mydata$d2peak[subscripts])
ld <- as.numeric(xlock - xdate1st)
xlock <- unique(
mydata$lockdown[subscripts])
t0 <- as.numeric(xlock - xdate1st) - 10
d2peak <- unique(mydata$d2peak[subscripts]) - t0
ld <- -10
## draw the confidence band
o <- order(x); X <- x[o]; Y <- y[o]
pg <- cbind(
......
Supports Markdown
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