functions.R 24.9 KB
Newer Older
1
2
3
4
5
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
## 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"
51
52
53
get_world_map <- function(x){
    ans <- readOGR(dsn = x)
    ans$name <- standardise_country_names(ans$ADM0_NAME)
54
55
56
    return(ans)
}

57

58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#' 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)
73
  return(geosimplify(raster::intersect(context, E)))
74
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
#' 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
104

105
106
107
108
109
110
111
112
113
114
  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",
115
                  "country", "geoid", "code", "pop", "continent", "XXX")
116
117
118
  ans$country <- gsub(pattern = "ô",
                      replacement = "o",
                      x = as.character(ans$country))
119
120
121
122
123
124
  ans$name <- gsub(
    pattern = "_",
    replacement = " ",
    x = as.character(ans$country)
  ) %>%
    standardise_country_names()
Facundo Muñoz's avatar
Facundo Muñoz committed
125
126
127
128

  return(ans)
}

129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
#' 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 <-
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
      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",
167
                  "The former Yugoslav Republic of Macedonia", "North Macedonia",
168
169
170
171
172
                  "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")
173
  std_name <- function(x) {
174
    if ( !is.na(idx <- match(x, country_names_table$variant)) ) {
175
        x <- country_names_table$standard_name[idx]
176
177
178
179
180
181
182
183
184
    }
    return(x)
  }

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

185

186
187
188
189
190
191
192
193
#' 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)]
}


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


205
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
## #' 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
237

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

#' Authors
241
authors <- function(){
242
243
244
  paste(
    c("[N. Alexander (Ergo)](neil.alexander@zoo.ox.ac.uk)",
      "[E. Arsevska (Cirad)](elena.arsevska@cirad.fr)",
245
      "[W. Van Bortel (ITM)](wvanbortel@itg.be)",
246
247
248
249
      "[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)",
250
      "[F. Mu\u00f1oz (Cirad)](facundo.munoz@cirad.fr)",
251
      "[W. Wint (Ergo)](william.wint@gmail.com)"),
252
253
    collapse=", ")
  }
254

255
256
257
258
259
260
261
262

#### 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")
}
263

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

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

295
296
297
298
299
300
301
302
303
304
## ## 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
##     }

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

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

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

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

## fn to compute the trend in mortality growth rate
418
trend <- function(x){
419
    ## use function "by" to loop over the country names
420
421
    xnam <- unique(x$name)
    xlock <- if(is.null(x$lockdown)) NA else unique(x$lockdown)
422
423
    ## initialize variables which might be updaated later or not
    ## depending whether specific conditions are met
424
425
426
427
428
429
430
431
432
433
434
435
436
437
    ## 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
438
    maxgr <- dr0 <- dr <- NA
439
440
441
442
443
444
445
446
    T95 <- T90 <- T80 <- T70 <- T60 <- T50 <- NA
    T40 <- T30 <- T20 <- T10 <- NA

    cond <- NA

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

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

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

        ## 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
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
            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)),
519
520
521
                    obsdr = tapply(x$deaths,
                                   x$day,
                                   sum) / x$pop,
522
                    fitdr = dmr,
523
524
525
                    cumdr = cumsum(tapply(x$deaths,
                                          x$day,
                                          sum) / x$pop),
526
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
                    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))
                }
            }
575
        }
576
577
    }
}
578

579
580
## function to visualize the mortalilty growth rates, as estimated by
## the function trend() described above. The arguments are: gr: the data
581
582
583
584
585
586
587
588
## 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.

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

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