functions.R 24.9 KB
Newer Older
1
2
3
## get lockdown data https://covidtracker.bsg.ox.ac.uk/
## Blavatnik School of Government / Oxford Univ
## Oxford COVID-19 Government Response Tracker
4
5
#stringency <- function(date_range = "2020-01-01/2020-05-31"){
stringency <- function(date_range){
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
    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"
52
53
54
get_world_map <- function(x){
    ans <- readOGR(dsn = x)
    ans$name <- standardise_country_names(ans$ADM0_NAME)
55
56
57
    return(ans)
}

58

59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
#' Simplify a vector layer
geosimplify <- function(x) {
  sx <- gSimplify(x, tol=.025, topologyPreserve=T)
  ans <- as(sx, "SpatialPolygonsDataFrame")
  ans@data <- x@data
  return(ans)
}

#' Extract the region of interest from a context
#'
#' Get a map of the extent of x from a context
region_of_interest <- function(x, context) {
  e <- extend(extent(x), 1)
  E <- as(e, "SpatialPolygons")
  projection(E) <- projection(x)
74
  return(geosimplify(raster::intersect(context, E)))
75
76
77
}


78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
#' ECDC Covid19 data worldwide
#'
#' Returns the published file at a given date.
#'
#' Downloading the file published a month ago is not the same as downloading
#' the latest file and filtering dates up to a mont ago!
#' I'm not sure which is best, but I can only presume that the latest file
#' will contain more accurate figures.
#'
#' @example
#'   covid202004 <- ecdc_covid19("2020-04-30")
#'   covid202005 <- ecdc_covid19("2020-05-20")
#'
#'   covid202004_latest <- covid202005 %>% filter(dateRep < "2020-05-01")
#'   identical(covid202004, covid202004_latest)
#'
#'   ## See differences. There are of all kind.
#'   full_join(
#'     covid202004,
#'     covid202004_latest,
#'     by = c("dateRep", "cases", "deaths", "countriesAndTerritories")
#'   ) %>%
#'   filter_all(any_vars(is.na(.)))
#'
ecdc_covid19 <- function(date) {
  base_url <- "https://www.ecdc.europa.eu/sites/default/files/documents/COVID-19-geographic-disbtribution-worldwide-"
  ecdc_url <- paste0(base_url, date, ".xlsx")
Facundo Muñoz's avatar
Facundo Muñoz committed
105

106
107
108
109
110
111
112
113
114
115
  tf <- tempfile(
    pattern = paste0("COVID-19-geographic-disbtribution-worldwide-", date),
    fileext = ".xlsx"
  )

  download.file(ecdc_url, tf, quiet = T, method = "curl")
  ans <- read_excel(tf)

  ## Cleanup variable and country names
  names(ans) <- c("date", "day", "month", "year", "cases", "deaths",
116
                  "country", "geoid", "code", "pop", "continent", "XXX")
117
118
119
  ans$country <- gsub(pattern = "ô",
                      replacement = "o",
                      x = as.character(ans$country))
120
121
122
123
124
125
  ans$name <- gsub(
    pattern = "_",
    replacement = " ",
    x = as.character(ans$country)
  ) %>%
    standardise_country_names()
Facundo Muñoz's avatar
Facundo Muñoz committed
126
127
128
129

  return(ans)
}

130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
#' Standardise country name variants
#'
#' This function maintains internally a table of country name variants and
#' their corresponding standard names.
#'
#' If a name is not found in the table, it is left as is.
#'
#' @param x Character vector of country names.
#'
#' @return A character vector of the same size with standardised names.
#'
#' @examples
#'   standardise_country_names(c("Spain", "Republic of Moldova", "Czechia", "Czechia", "France"))
standardise_country_names <- function(x) {
  ## Table of country name variants and their standardised name.
  ## Whenever a new country name causes conflict, update this table.
  country_names_table <-
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
      tibble::tribble(
                  ~variant, ~standard_name,
                  "Syrian Arab Republic", "Syria",
                  "State of Palestine", "Occupied Palestinian Territory",
                  "Czechia", "Czech Republic",
                  "Republic of Moldova", "Moldova",
                  "Democratic Republic of the Congo", "Congo DR",
                  "Congo, DRC", "Congo DR",
                  "Congo DRC", "Congo DR",
                  "Congo, DR", "Congo DR",
                  "Congo - Kinshasa", "Congo DR",
                  "Republic of Congo", "Congo",
                  "Congo - Brazzaville", "Congo",
                  "Gambia", "The Gambia",
                  "Cote d'Ivoire", "Cote dIvoire",
                  "Côte dIvoire", "Cote dIvoire",
                  "Côte d'Ivoire", "Cote dIvoire",
                  "Côte d'Ivoire", "Cote dIvoire",
                  "Cabo Verde", "Cape Verde",
                  "Guinea-Bissau", "Guinea Bissau",
                  "Macedonia", "North Macedonia",
168
                  "The former Yugoslav Republic of Macedonia", "North Macedonia",
169
170
171
172
173
                  "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")
174
  std_name <- function(x) {
175
    if ( !is.na(idx <- match(x, country_names_table$variant)) ) {
176
        x <- country_names_table$standard_name[idx]
177
178
179
180
181
182
183
184
185
    }
    return(x)
  }

  return(
    unname(vapply(x, std_name, character(1)))
  )
}

186

187
188
189
190
191
192
193
194
#' Get ISO3 codes for a list of countries
get_iso3 <- function(x) {
  codes <- ccodes()
  codes$ISO3[standardise_country_names(codes$NAME) %in%
               standardise_country_names(x)]
}


195
196
197
198
199
200
201
202
203
204
205
## 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"))
}


206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
## #' Total population by country for a given year
## #'
## #' Population counts in thousands of people.
## #'
## #' Some country names were corrected to match those from ECDC.
## #' Population count for Kosovo was manually added.
## #'
## #' @return Data frame with one row per country. Country name and population count.
## pop_countries <- function(x, year) {

##   ## Filter year and select variables Country and Population count.
##   pop <- base::subset(
##     x,
##     subset = Time == year,
##     select = c("Location", "PopTotal")
##   )
##   names(pop) <- c("name", "pop")
##   rownames(pop) <- NULL

##   ## Standardise country names
##   pop$name <- standardise_country_names(pop$name)

##   ## Manually include missing population for Kosovo
##   ## Source ???
##   pop <- rbind(
##     pop,
##     data.frame(name = "Kosovo", pop = 1845)
##   ) %>%
##     arrange(name)

##   return(pop)
## }
Facundo Muñoz's avatar
Facundo Muñoz committed
238

239
240
241
current_date <- function() format(Sys.Date(), "%d %b %Y")

#' Authors
242
authors <- function(){
243
244
245
  paste(
    c("[N. Alexander (Ergo)](neil.alexander@zoo.ox.ac.uk)",
      "[E. Arsevska (Cirad)](elena.arsevska@cirad.fr)",
246
      "[W. Van Bortel (ITM)](wvanbortel@itg.be)",
247
248
249
250
      "[S. Falala (Inrae)](sylvain.falala@cirad.fr)",
      "[E. Van Kleef (ITM)](evankleef@itg.be)",
      "[R. Lancelot (Cirad)](mailto:renaud.lancelot@cirad.fr)",
      "[C. Marsboom (Avia-GIS)](cmarsboom@avia-gis.com)",
251
      "[F. Mu\u00f1oz (Cirad)](facundo.munoz@cirad.fr)",
252
      "[W. Wint (Ergo)](william.wint@gmail.com)"),
253
254
    collapse=", ")
  }
255

256
257
258
259
260
261
262
263

#### Set locale for time representation. Platform-dependent. ####
if ( .Platform$OS.type == "unix") {
  Sys.setlocale("LC_TIME", "en_GB.UTF-8")
} else {
  ## Windows
  Sys.setlocale("LC_TIME", "English")
}
264

265
266
## utility for text output : takes a vector of names (argument nam)
## and a vector of counts (cnt), and returns a string with the names
267
## and counts within brackets, separated with commas, with "and" for
268
269
## the last one

270
collapse <- function(nam = onam, cnt = nd){
271
    strg <- as.character(nam)
272
    if(!any(is.na(cnt)) & length(cnt) > 1)
273
       if(length(strg) != length(cnt))
274
275
276
277
278
           stop('The lengths of names and counts are different.\n')
    n <- length(strg)
    if(all(!is.na(cnt))){
        strg <- unlist(sapply(strg,
                              function(x){
279
280
                                  o <- which(strg == x)
                                  out <- paste0(strg[o], " (", cnt[o], ")")
281
282
283
284
285
286
287
288
289
290
291
292
                              }))
    }
    if(n > 2){
        butlast <- paste(strg[-n], collapse = ", ")
        out <- paste0(butlast, " and ", strg[n])
        }
    else{
            if(n == 2)
                out <- paste(strg, collapse = " and ")
            else
                out <- strg
        }
293
    return(out)
294
295
    }

296
297
298
299
300
301
302
303
304
305
## ## fn to plot trellis plotd by column
##     packet.panel.bycolumn <- function(layout,
##                                       condlevels,
##                                       page,
##                                       row,
##                                       column,
##                                       skip) {
##         x <- c(layout[1] * (page - 1) + column, row); if (x[1] <= length(condlevels[[1]])) x
##     }

306
## fn to plot the daily mortality rate
307
308
dmrplot <- function(dmr, start = as.Date("2020-03-01"),
                    end = Sys.Date(), cex, nam, asp = 2/3){
309
310
311
312
    z <- subset(dmr, name %in% nam)
    idv <- c("name", "date", "lockdown", "d2peak")
    Z1 <- reshape2::melt(data = z,
                        id.vars = idv,
313
                        measure.vars = c("fitdr", "cumdr"))
314
315
    Z2 <- reshape2::melt(data = z,
                    id.vars = idv,
316
                    measure.vars = "obsdr")
317
    ## incubation 5 d
318
319
    ## early death 6 d after onset, late death 46 d after
    ## median death 14 days after onset
320
321
    ## left side at lockdown + 5 + 6, right side at lockdown + 5 + 46
    ## center at morkdown + 5 + (46 - 6) / 2
322
323

    ## regular xyplot
324
325
326
    ## key
    mykey <- list(
        space = "right", border = 1,
327
328
        text = list(c("Data",
                      "Trend",
329
                      "lockdown",
330
331
332
333
                      "Post-lockdown\nd 11 and 46",
                      "Peak"), cex = cex),
        lines = list(lty  = c(1, 1, 4, 1),
                     col  = c(grey(.8), 2, 1, 2),
334
335
                     size = 3))

336
337
338
339
340
    ## data for 100,000 inh.
    Z1$value100k <- 1e5 * Z1$value
    Z2$value100k <- 1e5 * Z2$value

    g <- xyplot(value100k ~ date | variable + name,
341
342
                xlim = c(as.Date("2020-02-15"), Sys.Date()),
                data = Z1, subscripts = T,
343
344
345
346
347
                scales = list(cex = cex - .1, tck = .5,
                              x = list(rot = 0,
                                       relation = "same"),
                              y = list(rot = 0,
                                       relation = "free")),
348
349
                layout = c(2, length(nam)),
                key = mykey,
350
351
352
353
354
                aspect = asp,
                between = list(x = 1/4, y = 1/4),
                xlab = character(0),
                ylab = list("Deaths / 100,000 inh.",
                            cex = cex),
355
                par.strip.text = list(cex = cex - .1),
356
357
                as.table = T,
                panel = function(x, y, subscripts, ...){
358
359
                    xnam  <- unique(Z1$name[subscripts])
                    xvar  <- unique(Z1$variable[subscripts])
360
361
                    if(xvar == "fitdr"){
                        ## add observed daily death rate
362
                        X <- Z2$date[Z2$name == xnam]
363
                        Y <- Z2$value100k[Z2$name == xnam]
364
                        o <- order(X)
365
366
367
368
                        lpoints(X[o], Y[o],
                                type = "h", col = grey(.8))
                        ## Add trend
                        llines(x, y, col = 2)
369
                    }
370
371
372
373
374
375
376
377
378
379
                    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, ")")
380
381
382
                        grid.text(
                            x = unit(.5, "mm"),
                            y = unit(1, "npc") - unit(.5, "mm"),
383
                            label = eval(parse(text = lbl)),
384
385
                            hjust = 0, vjust = 1,
                            gp = gpar(cex = cex - .1))
386
387
388
389
390
391
392
393
394
                        ## 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))
395
396
                    panel.abline(v = ld + c(0, 11, 46),
                                 col = c(1, 2, 2),
397
398
399
400
401
402
403
404
405
                                 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,
406
                                     col = 2)
407
                    }
408
                })
409
    ## plot with outer strips for better use of space
410
411
412
413
414
    p <- useOuterStrips(
        g,
        strip = strip.custom(
            factor.levels = c("Daily death rate",
                              "Cumulative daily death rate")))
415
    return(p)
416
417
418
}

## fn to compute the trend in mortality growth rate
419
trend <- function(x){
420
    ## use function "by" to loop over the country names
421
422
    xnam <- unique(x$name)
    xlock <- if(is.null(x$lockdown)) NA else unique(x$lockdown)
423
424
    ## initialize variables which might be updaated later or not
    ## depending whether specific conditions are met
425
426
427
428
429
430
431
432
433
434
435
436
437
438
    ## Y fitted values (growth rate)
    ## lo lower confidence bound of fitted valiue
    ## hi upper ""
    ## d2peak: # d from 1st death to the peak
    ## ldpeak number of days between lockdown and the peak
    ## ldmax # d after the lockdowm when the growth rate was max
    ## r19 Relative mortality increase from 19 days after lockdown to the peak
    ## a intercept of the tangent to fitted growth rate at 19 d after lockdown
    ## b slope of the tangent ""

    ## X25, X50, X75, X90: days when mortality decrease was 25%, 50%,
    ## 75%, 90% after the peak dr (vector), dr0 (max) daily death rate

    Y <- lo <- hi <- r19 <- a <- b <- d2peak <- ldpeak <- NA
439
    maxgr <- dr0 <- dr <- NA
440
441
442
443
444
445
446
447
    T95 <- T90 <- T80 <- T70 <- T60 <- T50 <- NA
    T40 <- T30 <- T20 <- T10 <- NA

    cond <- NA

    ## initialize data containers
    fittedTrend <- fittedParam <- data.frame()

448
    k <- length(seq(min(x$date), Sys.Date(), by="1 week"))
449
    ## fits a negative binomial additive model for death coints
450
    fm <- try(gam(deaths ~ s(day, bs = "cr", k = k) + offset(log(pop)),
451
              data = x, family = nb, link = "log",
452
              optimizer = "perf"))
453
    if(any(class(fm) == "try-error"))
454
455
456
457
458
        return(NULL)
    else{
        pred <- predict(fm, se = T)
        Y <- exp(pred$fit)
        ## these are counts we want rates
459
        dmr <- Y / x$pop
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478

        ## 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
479
        gr11 <- dY[day == dlock + 11]
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495

        ## relative change in mortality growth rate, reference value = gr11
        rgr19 <- 100 * dY[day = dlock + 19] / gr11
        rgr46 <- 100 * dY[day = dlock + 46] / gr11

        ## 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
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
            l1  <- 0
            if(!any(day >= dlock & dY < 0))
                return(NULL)
            else{
                d2peak <- min(day[day >= dlock & dY < 0], na.rm = T)
                l1 <- length(d2peak)

                ## lockdown to peak
                ldpeak <- d2peak - dlock

                ## max fitted death rate
                dr0 <- dmr[day == d2peak]

                ## decay
                Y2 <- Y[day >= d2peak]
                Ymax <- Y[day == d2peak]
                d2 <- day[day >= d2peak] - d2peak

                ## collects the results
                fittedTrend <- data.frame(
                    name = x$name,
                    date = x$date,
                    day  = day,
                    date1st = rep(min(x$date), nrow(x)),
520
521
522
                    obsdr = tapply(x$deaths,
                                   x$day,
                                   sum) / x$pop,
523
                    fitdr = dmr,
524
525
526
                    cumdr = cumsum(tapply(x$deaths,
                                          x$day,
                                          sum) / x$pop),
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
                    pop = unique(x$pop),
                    lockdown = rep(xlock, nrow(x)),
                    d2peak = rep(if(l1 > 0) d2peak else NA, nrow(x)),
                    ldpeak = rep(if(l1 > 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

                issue <- l1 == 0
                if(!issue){
                    fittedParam <- data.frame(
                        name     = unique(x$name),
                        lockdown = xlock,
                        d2peak   = d2peak,
                        ldpeak   = ldpeak,
                        maxdr    = dr0,
                        gr11     = gr11,
                        rgr19    = rgr19,
                        rgr46    = rgr46,
                        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))
                }
            }
576
        }
577
578
    }
}
579

580
581
## function to visualize the mortalilty growth rates, as estimated by
## the function trend() described above. The arguments are: gr: the data
582
583
584
585
586
587
588
589
## set containing the growth rates (output of trend()); nc: the number
## of columns in the plot; asp: the aspect of each panel: see the
## argument aspect in xyplot(); cex: character expansion for the
## labesl and other figure annotations; nam: a vector of strings
## containing the names of the countries names to be included in the
## plot. Must be a subset of the country names present in the trend()
## output.

590
591
592
grplot <- function(gr, nc, asp=2/3, cex=.7, dYlim=NA, nam){
    ## legend
    mykey <- list()
593
594
595
596
597
598
    mykey$space <- "right"
    mykey$divide <- 1
    mykey$rep <- F
    mykey$padding.text <- 1.5
    mykey$background <- grey(.95)
    mykey$border <- grey(.5)
599
600
601
602
603
    mykey$text <- list(c(
        "Trend",
        "95% conf. band",
        "Lockdown start",
        "Peak line (y=0)",
604
        "Early / late\nlockdown effect",
605
        "Mortality peak"), cex = cex + .1)
606
    mykey$lines <- list(type = "b",
607
608
609
                        pch = c(NA, 15, NA, NA, 15, 4),
                        col = c(4, 4, 1, 4, "pink", 4),
                        lty = c(1, 0, 4, 2, 0, 4),
610
                        size = 3, lwd = 1.5)
611
    mydata <- subset(gr, name %in% nam)
612
    ## set time = 0 to lockdown
613
614
615
616
    myL <- by(mydata,
              list(name = mydata$name),
              function(x){
                  x$dayL <- as.numeric(x$date - x$lockdown)
617
                  subset(x, dayL >= -10)
618
619
620
621
622
623
624
              })

    mydata <- do.call("rbind", myL)
    mydata$name <- reorder(mydata$name,
                           mydata$cumdr,
                           function(x)
                               1/max(x, na.rm = T))
625
    nr <- ceiling(length(nam) / nc)
626
    extrg <- extendrange(range(mydata$dY, na.rm=T))
627
    ylim <- if(any(is.na(dYlim))) extrg else dYlim
628
629
630
    xyplot(dY ~ dayL | name,
           data = mydata, subscripts = T,
           ylim = ylim,
631
632
           scales = list(cex = cex - .2, tck = .5,
                         x = list(relation = "same"),
633
634
635
636
637
                         y = list(rot = 0,
                                  relation = "same")),
           layout = c(nc, nr), as.table = T,
           between = list(x = 1/4, y = 1/4),
           aspect = asp,
638
           key = mykey,
639
           xlab = list("Days post lockdown\n", cex = cex),
640
641
642
           ylab = list("Mortality growth rate", cex = cex),
           par.strip.text = list(cex = cex),
           panel = function(x, y, subscripts){
643
644
645
646
647
648
649
650
651
652
                   ## lockdown effetc
                   grid.rect(
                       x = unit(11 + (46 - 11) / 2, "native"),
                       y = unit(0.5, "npc"),
                       width = unit(46 - 11, "native"),
                       height = unit(1, "npc"),
                       gp = gpar(fill = "pink",
                                 col = NA,
                                 alpha = .5)
                   )
653
654
               xlock <- unique(
                   mydata$lockdown[subscripts])
655
656
               ldpeak <- unique(
                   mydata$ldpeak[subscripts])
657
               ## draw the confidence band
658
659
660
               o <- order(x); X <- x[o]; Y <- y[o]
               pg <- cbind(
                   c(X, rev(X)),
661
662
                   c(mydata$dYlo[subscripts][o],
                     rev(mydata$dYhi[subscripts][o])))
663
               panel.polygon(pg, border = NA,
664
                             col = "light blue", alpha = .7)
665
666
               ## smoothed mean
               panel.xyplot(X, Y, type = "l", col = 4)
667
               if(!is.na(xlock)){
668
669
670
671
672
673
                   ## lockdown
                   panel.abline(v = 0, col = 1, lty = 4)
                   ## peak
                   panel.abline(v = ldpeak, col = 4, lty = 4)
                   ## death rate threshold
                   panel.abline(h = 0, col = 4, lty = 4)
674
675
676
           }
    })
}