Covid19MortalityEurope.Rmd 57.3 KB
Newer Older
1
---
2
title: |
3
4
5
    | ![](logoMoodWideBlue.png){width=2.5in}  
    |
    | Visualizing COVID-19 mortality in Europe in near-real time
6
date: "Last updated: `r date()`"
7
8
9
10
11
12
13
14
15
16
author:
  - Lancelot, R^[renaud.lancelot@cirad.fr]
  - Muñoz, F^[facundo.munoz@cirad.fr]
  - Van Bortel, W^[wvanbortel@itg.be]
  - Wint, WR^[william.wint@gmail.com]
  - Alexander, N^[neil.alexander@zoo.ox.ac.uk]
  - Marsboom, C^[cmarsboom@avia-gis.com]
  - Arsevska, E^[elena.arsevska@cirad.fr]
  - Falala, S^[sylvain.falala@cirad.fr]
  - Van Kleef, E^[evankleef@itg.be]
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
Affiliation:
  - "AVIA-GIS":
    name: "Marsboom, C"
  - "ERGO":
    name: "Alexander, N; Wint, WR"
  - "Institute of Tropical Medicine": 
    name: "Van Kleef, E; Van Bortel, W"
  - "ASTRE":
    name: "Arsevska, E; Falala, S; Lancelot, R; Muñoz, F"
  - "CIRAD":
    name: "Arsevska, E; Lancelot, R; Muñoz, F"
  - "INRAE":
    name: "Falala, S"
address:
  - "Zoersel; Belgium":
    code: "AVIA-GIS"
  - "Institute of Tropical Medicine, Antwerp, Belgium":
    code: "ITM"
  - "Environmental Research Group of Oxford, Oxford, UK": 
    code: "ERGO"
  - "UMR ASTRE, Univ Montpellier (I-MUSE), CIRAD, INRAE, Montpellier, France":
    code: "ASTRE"
  - "UMR ASTRE, Montpellier, France":
    code: "INRAE"
41
fig_caption: yes
42
bibliography: covid.bib
43
44
45
citation_package:
  natbib:
    natbiboptions: "numbers, sort&compress"
46
47
csl: plos-one.csl
fontsize: 12pt
48
geometry: "margin=1in"
49
---
50

51
```{r config, eval = interactive(), cache = FALSE, include = FALSE}
52
53
54
55
56
57
58
59
## This chunk is only executed in interactive mode (i.e. not when compiled)
## and helps to debug the RMarkdown starting from a fresh session.
source("src/packages.R")
source("src/functions.R")
message("Assuming interactive session")
```


60
61
62
```{r settings, include=F}
knit_hooks$set(optipng = hook_optipng)
knit_hooks$set(pngquant = hook_pngquant)
63
knitr::opts_chunk$set(fig.pos = 'H')
64
knitr::knit_hooks$set(crop = knitr::hook_pdfcrop)
65
66

## max width in inches for plots for A4 paper size (21 * 29.7 cm) with
67
## 1-inch margins (top/bottom and left/right)
68
69
W <- 21 / 2.54 - 2

70
71
knitr::opts_chunk$set(
    cache = F,
72
73
74
75
76
77
78
79
80
81
##    dev.args = list(png = list(res = 600)),
    fig.retina = 2,
    fig.align = "center",
    fig.width = W,
    fig.height = W,
    out.width = "100%",
    crop = T,
    include = F,
    echo = F,
    warning = F,
82
    comment = NA)
83

84
pdf.options(family = "Palatino")
85
86
87

```

88
```{r load_targets}
89

90
loadd(
91
    world_map,
92
93
    geo_europe_light,
    roi_europe_light,
94
    xeurope_countries,
95
    sdt,
96
97
    wgs84,
    merc,
98
    eust,
99
    pop_countries_2020, 
100
    covid,
101
102
    lockdown_dates,
    lock
103
)
104
105
106
```

```{r lockdown}
107
## lockdown
108
109
Lock <- do.call(
    "rbind",
110
111
    by(subset(lock, name %in% xeurope_countries$name),
       list(name = subset(lock, name %in% xeurope_countries$name)$name),
112
113
114
115
116
117
118
       function(x){
           smax <- max(x$stringency, na.rm = T)
           x$T90 <- min(x$date[x$stringency >= .9 * smax], na.rm = T)
           return(x)
       })
)

119
deaths_by_ctry_dt <- covid %>%
120
    filter(deaths >= 0) %>% 
121
122
123
124
125
126
    group_by(date, name) %>%
    summarise(deaths = sum(deaths)) %>%
    ungroup() %>%
    mutate(date = as.Date(date))

data_update_date <- format(Sys.Date() - 1,
127
                           "%d %b %Y")
128

129
## clean the data
130
covid$name <- standardise_country_names(covid$name)
131

132
133
134
135
## Variables relative to European countries: 
## group, name, lockdown date, population, cumulated deaths,
## and death rate.
country_data <- 
136
  xeurope_countries %>% 
137
  ## filter(name != "Liechtenstein") %>% 
138
139
140
141
142
  left_join(
    lockdown_dates,
    by = "name"
  ) %>% 
  left_join(
143
    pop_countries_2020,
144
145
    by = "name"
  ) %>% 
146
147
148
149
150
151
152
153
154
  mutate_if(is.character, as.factor) %>% 
  left_join(
    deaths_by_ctry_dt %>% 
      group_by(name) %>% 
      summarise(
        cum_deaths = sum(deaths),
        .groups = "drop"
      ),
    by = "name"
155
  ) %>% 
156
157
158
  mutate(
    ## Cumulated number of deaths per 100 K hab.
    ## Note that pop is expressed in thousands.
159
    death_rate = 1e5 * (cum_deaths / pop)
160
161
  ) %>%
  ## Arrange rows by descending order of death rate
162
    arrange(desc(death_rate))
163

164
165
166
167
168
169
170
171
172
173
174
175
## Deaths by country and date of European countries
## with population and lockdown date added
Dfr <- deaths_by_ctry_dt %>% 
  inner_join(
    country_data,
    by = "name"
  ) %>% 
  mutate(
    # Order name levels in decreasing order of death rate
    name = fct_reorder(name, death_rate, .desc = TRUE)
  ) %>% 
  select(name, date, deaths, pop, lockdown)
176

177
178
## Names of countries with less than 25 cumulated deaths
nam25 <- country_data %>% 
179
  filter(cum_deaths < 25) %>% 
180
  pull(name)
181

182
## character string with the number of deaths by country
183
184
sond <- format(country_data$cum_deaths,
               big.mark = ",", trim = T)
185

186
## character string with the cumulative death rate by country
187
188
socdr <- sprintf(country_data$death_rate, fmt = "%##.1f",
                 digits = 0)
189

190
## overall number of deaths
191
ttnd <- sum(Dfr$deaths)
192

193
194
## character string for overall number of deaths
sttnd <- format(ttnd, big.mark = ",")
195

196
## overall population size (in thousands)
197
ttpop <- sum(country_data$pop)
198

199
## character string for overall pop size
200
sttpop <- format(round(ttpop, 0), big.mark=",")
201

202
## overall cumulative death rate for 100,000 inh.
203
ttcdr <- 1e5 * ttnd / ttpop
204

205
206
## character string for the overall cumulative death rate
sttcdr <- format(ttcdr, digits = 3)
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
```{r spatdata}
## spatial data
covmap <- geo_europe_light
covmap2 <- spTransform(covmap, merc)

## change the extent (units in degrees)
e <- raster::extend(raster::extent(covmap),
                    y = c(-5, 6, -2, 1))
E <- as(e, "SpatialPolygons")
projection(E) <- wgs84

roi <- roi_europe_light
roi2 <- spTransform(roi, merc)

## cumulative mortality
covmap2$cumdr <- NA
## In the series of variables onam_xxx I keep records of the list of
## study countries ordered by decreasing mortality rates
onam <- levels(Dfr$name)

for(i in seq.int(onam)){
    xname <- onam[i]
    xcumdr <- country_data$death_rate[which(country_data$name == xname)]
234
    covmap2@data$cumdr[which(covmap2@data$name == xname)] <- xcumdr
235
}
236

237
238
239
```


240
```{r onam2}
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
## add population size
## we discard countries with < 100 deaths
onam2 <- onam[!(onam %in% nam25)]

## the condition on deaths >= 0 is needed because Spain has declared a
## negative count in early June
Dfr2 <- subset(Dfr, name %in% onam2 & deaths >= 0)
Dfr2$name <- factor(as.character(Dfr2$name),
                    levels = onam2)

Dfr3  <- do.call(
    "rbind",
    by(Dfr2,
       list(name = Dfr2$name),
       function(x){
           x$day = as.numeric(x$date - min(x$date, na.rm = T))
           x
       }))
259

260
261
262
263
264
265
266
267
268
269
270
271
272
273
Dfr3$name <- factor(as.character(Dfr3$name),
                    levels = onam2)

## use function "by" to loop over the country names
L2 <- by(Dfr3,
         INDICES = list(name = Dfr3$name),
         FUN = function(x) trend(x))

## fitted growth rate
fTrend <- do.call("rbind",
                  lapply(L2,
                         function(x)
                             x$fittedTrend))
rownames(fTrend) <- seq(nrow(fTrend))
274

275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
## indicators
fParam <- do.call("rbind",
                  lapply(L2,
                         function(x)
                             x$fittedParam))
rownames(fParam) <- seq(nrow(fParam))

## for plotting
Dfr4 <- fTrend

## assess whether some countries did not peak
list_of_countries_after_trend <- unique(Dfr4$name)
did_not_peak <- setdiff(onam2, list_of_countries_after_trend)
onam3 <- onam2[onam2 %in% list_of_countries_after_trend]

## keek onam order for country names
Dfr4$name <- factor(as.character(Dfr4$name),
                    levels = onam3)

```

296

297
298
```{r prepastats}
## table preparation
299

300
301
302
303
304
zmydat1 <- na.omit(fParam[order(as.character(fParam$name)),
                          match(c("name",
                                  "d2peak", "maxdr",
                                  "gr11", "rgr46"),
                                table = names(fParam)), ])
305
306
307
308
309

mydat1 <- zmydat1[,-1]

rownames(mydat1) <- zmydat1$name

310
311
312
313
314
315
316
zmydat2 <- fParam[order(as.character(fParam$name)),
                  match(c("name", "maxdr", "T10", "T20", "T30",
                          "T40", "T50", "T60", "T70",
                          "T80", "T90", "T95"),
                        table = names(fParam)), ]
mydat2 <- zmydat2[!is.na(zmydat2$maxdr),-(1:2)]
rownames(mydat2) <- zmydat2$name[!is.na(zmydat2$maxdr)]
317
318
319
320
321
322
323
324
325
326
327

nopeak_statement <- character(0)

if(length(did_not_peak) == 0)
    nopeak_statement <- paste(
        "The mortality related to the COVID-19 epidemic",
        "is low in all the countries")

if(length(did_not_peak) > 0){
    nopeak_statement <- paste(
        "The mortality related to the COVID-19 epidemic",
328
        "is low in most countries.",
329
330
331
332
333
334
335
        "However, the epidemiological status is still uncertain in",
        collapse(did_not_peak, cnt=NA))

}

## countries with mortality >= 50% peak mortality
persist <- rownames(mydat2)[which(is.na(mydat2$T50))]
336

337
338
```

339
\newpage
340
341
342
343
344

# Summary {-}

## Goals{-}

345
346
347
348
349
The motivation of this work was to provide a near-real time
visualization tool of the COVID-19 related mortality in Europe,
highlighting the features of COVID-19 transmission dynamic, particular
to help knowing whether a peak in daily mortality rate was reached,
and whether mortality started decreasing.
350

351
352
353
354
In addition, we defined and estimated a set of indicators of the
country-level daily mortality pattern, to enable further studies on
the effect of lockdown measures, and their implementation (e.g.,
consequences of human mobility), on the course of this epidemic.
355
   
356
357
These indicators were then used to identify clusters of countries
sharing similar mortality patterns.
358
359
360

## Main results {-}

361
As per `r data_update_date`, in the `r length(xeurope_countries$name)` countries
362
363
accounted for, the total number of deaths is `r sttnd` for an overall
population size of `r round(ttpop/1e6)` millions representing a
364
cumulative death rate of `r sttcdr` deaths $10^{-5}$ inhabitants.
365

366
**The most affected countries are located in south-western Europe, at
367
368
369
the exception of Sweden, Switzerland, and North Macedonia**. Indeed,
the cumulative death rate is contrasted between these south-western
European countries and the rest od Europe.
370
  
371
372
Low death rates are reported in Finland, the Baltic countries, Central
Europe, and south-eastern Europe. These countries are heterogeneous in
373
374
375
terms of socio-economic features. For COVID-19, factors related to
sub-population connectiveness and vulnerability, might explain these
differences in disease incidence.
376
  
377
378
379
380
381
The highest mortality growth rates were met before the implementation
of lockdown measures, with some fluctuations around the lockdown
date - probably depending on the progressiveness of the preventive
measures decided by the governments, and their actual appropriation by
the populations. The mortality growth rate decreased
382
383
384
385
386
thereafter. Different situations were encountered after the first
mortality peak:

* **a continuous decrease in the mortality growth rate** after the
  peak (Belgium, Italy, France, the Netherlands, Ireland,
387
388
  Switzerland...): the daily mortality rate continuously decreased,
  and now stands at a low or very low level.
389
	
390
391
392
393
* the mortality growth rate increased again after the peak, possibly
  resulting in secondary and further peaks: Spain, North Macedonia,
  Portugal, Austria, Bosnia and Herzegovina, Turkey, Serbia, Kosovo,
  Montenegro, Albania, Czech Republic, and Croatia.
394
395

Daily mortality rate has now reached a low level in most countries
396
(fig. \@ref(fig:fdecay)). However, it is sill at least 50% of the
397
398
399
400
401
402
403
404
405
406
407
408
409
410
national peak mortality rate in `r collapse(persist, NA)`.
  
These observations suggest that while the overall situation has much
improved with respect to the worst days of the outbreak, it is still
heterogeneous and somewhat unstable. The epidemiological situation
still needs to be closely monitored.

Finally, mapping the mortality pattern shows three broad patterns of
decreasing severity: (i) western and southern Europe, (ii) a part of
northern and central Europe (Sweden, Romania...), and (iii) remaining
northern and central Europe (Norway, Finland, Germany,
Austria...). The next step of this work will be to assess whether
these patterns are associated with specific lockdown measures or their
implementation.
411

412
\newpage
413

414
415

## Geographical scope{-}
416

417
The geographical scope of the study is:
418

419
* the European Union and the UK (28 countries): `r collapse(xeurope_countries %>% filter(group == "EU28") %>% pull(name), cnt = NA)`,
420
  
421
* plus member countries of the European Free Trade Association (EFTA):
422
  `r collapse(xeurope_countries %>% filter(group == "EFTA") %>% pull(name), cnt = NA)`;
423
  
424
* plus candidate countries for EU membership: `r collapse(xeurope_countries %>% filter(group == "Candidates") %>% pull(name), cnt = NA)`;
425
  
426
* plus European, non-EU countries: `r collapse(xeurope_countries %>% filter(group == "Non_EU") %>% pull(name), cnt = NA)`.
427

428
Liechtenstein was excluded from the study because of its small population size.
429

430
## Approach{-}
431

432
433
434
* Publicly available data sets (see section \@ref(data)) are used to
  estimate national mortality indicators (see section \@ref(indi))
  starting with daily death rate (daily death count scaled by its
435
  population size [@Porta2008]). These data are heterogeneous in
436
437
  nature, by the fact that different case definitions, testing and
  reporting strategies are implemented in the included countries.
438

439
# Results
440

441
## Cumulative death rate
442

443
In the `r nrow(xeurope_countries)` countries accounted for^[List (number of deaths in braces): `r collapse(nam = levels(Dfr$name), cnt = sond)`.],
444
445
446
the total reported deaths is now `r sttnd` for an overall population size of 
`r round(sum(country_data$pop)/1e3,0)` millions, thus representing a cumulative death rate of 
`r round(100 * sum(Dfr$deaths) / sum(country_data$pop), 1)` deaths $10^{-5}$ inhabitants (inh.).
447
448
449
These figures and ranking must be considered very cautiously because
case definition and reporting strategies were heterogeneous among the
study countries.
450

451
### By country{-}
452

453
```{r cumul, include = T, fig.asp = 1.25, out.width="60%", fig.cap = cap}
454

455

456
cap <- paste("Cumulative apparent COVID-19 death rate (total",
457
458
459
460
             "deaths in braces) from the start of the epidemic to",
             data_update_date,
             "in Europe, as reported by the nationak public-health agencies.",
             "South-western European countries are labelled in red.")
461
462
463
464

red  <- c("Belgium", "France", "Ireland", "Italy", "Luxembourg",
          "Portugal", "Spain", "Netherlands", "United Kingdom")

465
lbl <- paste0("(", sond, ")")
466
cex <- 0.85
467

468
dotplot( ~ rev(setNames(death_rate, name)),
469
470
471
472
473
474
475
476
477
478
479
        data = country_data,
        xlim = extendrange(range(country_data$death_rate, na.rm = T),
                           f = c(.05, .2)),
        scales = list(cex = cex, tck = .5,
                      y = list(labels = rev(country_data$name),
                               col = ifelse(
                                   rev(country_data$name) %in% red,
                                   2,
                                   1)
                               )
                      ),
480
        cex = cex,
481
        col = ifelse(rev(country_data$name) %in% red, 2, 1),
482
        aspect = 1.6,
483
        xlab = list("Deaths / 100,000 inh.",
484
485
486
487
488
489
490
                    cex = cex)) +
    latticeExtra::layer(
      grid::grid.text(
        x = unit(x, "native") + unit(2, "mm"),
        y = unit(y, "native"),
        label = rev(lbl),
        hjust = 0, vjust = .5,
491
492
        gp = gpar(cex = cex,
                  col = ifelse(rev(onam) %in% red, 2, 1))
493
      ),
494
      data = list(lbl = lbl, cex = cex, red = red, onam = country_data$name)
495
    )
496
```
497

498
\newpage
499

500
### Spatial distribution{-}
501

502
```{r covmap, include=T, fig.height=7,fig.cap=cap}
503
cap <- paste0(
504
    "Cumulative apparent mortality rate (deaths / 100,000 inh.) ",
505
506
507
508
509
    "attributed to COVID-19 as reported by the national ",
    "public-health agencies in Europe: situation on ",
    data_update_date,
    ". Countries are listed in the legend, and numbered on ",
    "the map (cumulative reported deaths in braces).")
510

511
num <- sprintf(seq.int(nlevels(Dfr$name)), fmt = "%02d")
512

513
## colors for the map
514
515
cols <- colorRampPalette(rev(brewer.pal(n=9, "RdYlGn"))[-(1:2)])

516
labs <- paste0(num, " ", onam, " (", sond, ")")
517
nch <- nchar(labs)
518
slabs <- labs[which(nch == max(nch))]
519

520
ocovmap2 <- covmap2[covmap2$name %in% onam, ]
521
oXY <- data.frame(coordinates(ocovmap2))
522
names(oXY) <- c("x", "y")
523
oXY$name <- as.character(ocovmap2$name)
524
num <- unlist(sapply(oXY$name, function(x) which(onam == x)))
525
526
            
## parameters of the map
527
cex <- .75
528
e <- raster::extent(roi2)
529
530
H <- stringHeight(slabs)
L <- stringWidth(slabs)
531
n <- length(onam)
532
kex <- 1.1
533

534
535
spplot(ocovmap2, zcol = "cumdr",
       col.regions = cols(50),
536
       xlim = e[1:2], ylim = e[3:4],
537
       scales = list(draw = F),
538
       key = list(space = "right",
539
540
                  background = "transparent",
                  x=.005, y = .005,
541
                  text = list(labs, cex = cex-.1)),
542
       colorkey = list(space = "bottom", height = .7, width = .8),
543
       sub = list("Deaths / 100,000 inh.", cex = cex, font = 1),
544
       lwd = 1/4, col = grey(.5)) +
545
    latticeExtra::layer(sp.polygons(roi2, fill = grey(.9),
546
                      lwd = 1/4, col = grey(.5)),
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
                      data = list(roi2 = roi2), under = T) +
    latticeExtra::layer(grid.rect(
        x = unit(oXY$x, "native"),
        y = unit(oXY$y, "native"),
        width = kex * stringWidth(num),
        height = kex * stringHeight(num),
        gp = gpar(fill = "white", col = NA, alpha = .5)),
        data = list(oXY = oXY, kex = kex, num = num)) +
    latticeExtra::layer(
        grid.text(
            label = num,
            x = unit(oXY$x, "native"),
            y = unit(oXY$y, "native"),
            gp = gpar(cex = cex)),
        data = list(oXY = oXY, num = num, cex = cex))
562
563
564

```

565
* **The most affected countries are located in south-western
566
567
  Europe^[Belgium, France, Ireland, Italy, Luxembourg, Portugal, Spain, The
Netherlands, The United Kingdom.], at the exception of Sweden** (ranking 
568
  `r which(onam ==  "Sweden")` on fig. \@ref(fig:covmap)), where
569
570
571
572
  strict lockdown measures were not adopted. Indeed, the cumulative
  death rate (fig. \@ref(fig:cumul)) is contrasted between these
  south-Western European countries and the others. For instance, low
  death rates were reported in Germany
573
574
  (`r socdr[which(onam == "Germany")]`  $10^{-5}$ inh.), or even lower
  in Poland (`r socdr[which(onam == "Poland")]` $10^{-5}$ inh.).
575
  
576
577
* Low death rates are reported in Finland, the Baltic countries,
  Central Europe, and south-eastern Europe. These countries are
578
579
  heterogeneous in terms of socio-economic features. These differences
  are known to be risk factors for the emergence of infectious
580
581
582
583
584
  diseases, e.g., tick-borne encephalitis (TBE) [@EDEN0205,
  @Randolph2010]. However, TBE was much related to changes in
  agriculture and exposure to tick bites; other socio-economic factors
  are probably important for COVID-19, such as the age-structure o
  populations, their density and connectivity, as well as testing and
585
  reporting strategies. 
586
587
  
* For the rest of this document, we limited the analysis to countries
588
589
590
  totaling more than 25 deaths from the start of the epidemic - defined
  here as the date of the first reported death. Consequently, the
  following countries were discarded from the rest of the study:
591
  `r collapse(nam25, cnt=NA)`. Indeed, the estimation of trends in
592
593
  daily death rates over a time frame of more than
  `r round(as.numeric(Sys.Date() - min(Dfr$date)) / 30.44)` months seemed
594
595
  to be difficult for lower counts.
  
596

597
## Daily death rates
598

599
600
```{r deathRate, include = F}
cap <- paste("Daily death rate (count / 100,000 inh. / d)",
601
602
             "related to COVID-19 infection as reported by the",
             "national public-health agencies in Europe. The",
603
604
605
             "$c_i$",
             "value on the top left corner of each panel of the",
             "right column of plots is the cumulative death rate",
606
             "and the subscript is the timr (d) since the first reported",
607
             "death. The black, vertical line is",
608
609
610
             "drawn at the lockdown date. The blue, vertical lines",
             "are drawn at days 11 and 46 after the lockdown,",
             "corresponding to early and late lag",
611
             "between infection and death. Countries: ")
612
613

## number of rows (countries) / page
614
nbr <- 8
615
616
617
nbp <- ceiling(length(onam2) / nbr)

## generate as many captions as the number of pages. The caption is
618
## made of a common part (defined above) plus the list of countries
619

620
caps <- vector(mode = "character", length = nbp)
621
List <- list()
622
for(i in 1:nbp){
623
624
    b <- if(i < nbp) nbr * i else length(onam2)
    nam <-  onam2[(nbr*(i-1) + 1):b]
625
    List[[i]] <- nam
626
    caps[i] <- paste0(cap, collapse(nam, cnt = NA), ".")
627
}
628

629
630
```

631
```{r dmr1, include=T, fig.height=8, fig.cap=cap}
632

633
cap <- caps[1]
634

635
dmrplot(dmr = fTrend,
636
637
        cex = cex,
        asp = "fill",
638
        nam = List[[1]])
639

640
```
641

642
\newpage
643

644
```{r dmr2, include=T, fig.height=8,fig.cap=cap}
645
cap <- caps[2]
646

647
dmrplot(dmr = fTrend,
648
649
        cex = cex,
        asp = "fill",
650
        nam = List[[2]])
651
```
652

653
\newpage
654

655
```{r dmr3, include=T, fig.height=8, fig.cap=cap}
656

657
cap <- caps[3]
658

659
dmrplot(dmr = fTrend,
660
661
        cex = cex,
        asp = "fill",
662
        nam = List[[3]])
663

664
H <- length(List[[4]]) 
665
666
```

667
\newpage
668

669
```{r dmr4, include=T, fig.height=H, fig.cap=cap}
670

671
cap <- caps[4]
672
dmrplot(dmr = fTrend,
673
674
        cex = cex,
        asp = "fill",
675
        nam = List[[4]])
676
677
```

678

679
680
* Daily death rates (fig. \@ref(fig:dmr1) to \@ref(fig:dmr4), left
  column) show large day-to-day variations for a given country. A part
681
  of these variations is related to delays in data collection, and
682
683
684
  subsequently, to the data updates (corrections) reported by the
  national public-health agencies. Therefore, it is wiser to interpret
  the trends, rather than specific data points.
685
686

* On fig. \@ref(fig:dmr1) to \@ref(fig:dmr4), countries are ordered
687
  from the highest daily death rates on the top plot, left column of
688
689
  the first plot matrix (`r onam2[1]`), to the lowest daily death rate
  on the bottom plot, left column of the fourth plot matrix
690
691
  (`r  onam2[length(onam2)]`). Cumulative death rates are displayed on the
  right column of each plot matrix.
692

693
694
695
* According to the available data, `r onam2[1]` is the most severely
  hit country (apparent cumulative death rate: `r socdr[1]` deaths
  $10^{-5}$ inh.).
696
697

* The 46-d post-lockdown limit has been reached by all countries:
698
699
700
701
702
  fourth (if any) vertical, red, dashed line from the left of each
  plot. This means that - according to the maximum recorded delay
  between the infection and death times (46 days), all the deaths
  reported after this 46-day limit occurred in patients who were
  exposed to the virus after lockdown was implemented.
703
704

* Periodic variations of daily mortality rates are observed in several
705
  countries: see e.g., the pattern for the UK, Sweden, or Netherlands
706
  to be compared with the pattern of other countries on
707
708
709
710
711
712
713
  fig. \@ref(fig:dmr1). The cause of these variations would deserve to
  be investigated. At this stage, we can only assume it is either
  related to periodic virus transmission spikes (e.g., more mobility
  and less social distancing during the week-ends), or more likely, a
  reporting bias (e.g., less reporting during the week-ends). Anyway,
  the consequence is a higher uncertainty in the latest estimated values
  of the trend for these countries.
714

715
* Several countries should deserve close attention:
716
717
  
  * the daily mortality rate has much decreased with respect to the
718
719
720
    daily mortality peak, nut is still higher than at the lockdown:
    the United Kingdom, Sweden (fig. \@ref(fig:dmr1)), and Turkey
    (fig. \@ref(fig:dmr3));
721
722
723
	
  * the daily mortality rate has reached a peak, and decreased
    afterward. However it is increasing again: North Macedonia,
724
725
726
727
    Romania (fig. \@ref(fig:dmr2))), Bosnia and Herzegovina, Kosovo,
    Serbia (\@ref(fig:dmr3)), Poland, Bularia, and Croatia
    (fig. \@ref(fig:dmr4)).
	  
728
  \clearpage
729

730
## Time trends in daily mortality growth rate{-}
731

732
733
```{r rgr1, include=T, fig.height = 9, fig.cap=cap}
cap <- paste0(
734
    "Mortality growth rate related to COVID-19 in European countries ",
735
    "totalling at least 25 deaths. ",
736
737
738
739
740
    "The black line is drawn at the lockdown. ",
    "The pink bans spans over d 11 and 46 ",
    "after the lockdown, corresponding to early ",
    "and late time lag between infection and death. ",
    "The blue line is drawn at the mortality peak.")
741

742
## number of rows /page
743
nbr <- 9
744
nc <- 4
745
asp <- "fill"
746

747
748
749
750
fTrend$name <- reorder(fTrend$name,
                       fTrend$cumdr,
                       function(x)
                           1/max(x, na.rm=T))
751

752
nam <- levels(fTrend$name)
753

754
755
756
757
758
grplot(gr = fTrend,
       nc = nc,
       asp = asp,
       nam = nam,
       dYlim = c(-.15, .4))
759

760
761
```

762
National trends in daily mortality growth rates are shown on
763
fig. \@ref(fig:rgr1). The (first) peak in daily mortality rate was
764
765
observed when the daily mortality growth rate curve crossed the line
of equation $y=0$ (thereafter called "mortality threshold"). Negative
766
values of the growth rate correspond to decreasing mortality rates.
767

768
* In most cases, the daily mortality growth rate had started
769
770
771
  decreasing *before* the lockdown could have ab effect on mortality
  (i.e., before 11 days after the lockdown date), suggesting either
  the populations anticipated the official decisions, or the measures
772
773
774
775
776
  preceding the lockdown were efficient enough to start decreasing the
  daily mortality growth rate. We arbitrarily set the growth-rate
  reference on day 11 after the lockdown because this was the time
  when early death started being reported among patients who got
  infected at the lockdown date.
777
  
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
* The time interval spanning from 11 to 46 days after the lockdown
  generally included the mortality peak, suggesting the virus
  transmission was effectively reduced following the implementation of
  lockdown measures. Among the investigated national data sets,
  Bulgaria is the single country for which the peak was reached after
  the end of this interval. Nonetheless, the variability in peak date
  with respect to the lockdown is rather high, probably depending on
  how lockdown measures were appropriated and actually implemented by
  the populations. With this respect, cultural and social norms
  probably influenced the virus transmission.

* Different situations were met after the first mortality peak:

  + a continuous decrease in the daily growth rate: the daily
    mortality rate is continuously decreasing, and now stands at a low
    or very low level: Belgium, Italy, France, the Netherlands,
    Switzerland...
795
	
796
797
798
799
  + the daily mortality growth rate increased again, and eventually
    crossed the mortality threshold upward, indicating an increasing
    mortality rate: Spain, Portugal, North Macedonia, Romania,
    Austria, Turkey, Bosnia and Herzegovina, Poland, the Czech
800
    Republic, Serbia, Bulgaria, and Croatia.
801
802
803
804
805
806

These observations suggest that while the overall situation has much
improved with respect to the worst days of March and April, it is
still heterogeneous and somewhat unstable.

\clearpage
807
808
809
810
811

## Mortality patterns{-}

We used a panel of mortality indicators to characterize the national
mortality patterns, both for the overall features and growing step of
812
813
the daily mortality pattern (tab. \@ref(tab:tabg)), and for the decay
step (tab. \@ref(tab:tabd)):
814
815
  
* for the overall features and growing-mortality step:  
816

817
  + the estimated maximum daily death rate, i.e., the daily mortality
818
819
820
821
    rate at the first daily mortality peak (`maxdr`), and the time (d)
    when it was reached, with respect to the first reported death
    (`d2peak`);
	  
822
823
  + the daily mortality growth rate 11 days after the lockdown
    (`gr11`) before lockdown was implemented;
824

825
  + the relative change (\%) in fitted mortality growth rate on day 46
826
    after the lockdown, with respect to `gr11`: `rgr46`);
827

828
829
830
831
* for the daily mortality decay step, after the first daily mortality
  peak: the estimated times (d) from the mortality peak to a 10\%,
  20\%, ..., 90\%, and 95\% reduction of the daily mortality rate,
  with regard to the peak `T10`, `T20`,... `T95`.
832

833
\newpage
834

835
836
837
```{r tabg, include=T}

cap1 <- paste("Indicators of the COVID-19 mortality pattern:",
838
839
             "overall features and growing step.")

840
labcol1 <- paste("`d2peak`: time lag (d) from the 1st death to the mortality peak;",
841
                 "`maxdr`: daily mortality rate at the first daily mortality peak",
842
843
                 "(death number / 100,000 ind.),",
                 "`gr11`: mortality growth rate on day 11 after the lockdown;",
844
                 "`rgr46`: relative change (\\%) in fitted mortality",
845
                 "growth rate on d 46 after the lockdown, reported to `gr11`")
846

847
capg <- paste(cap1, collapse(labcol1, NA), "")
848
mydat1$maxdr100k <- 1e5 * mydat1$maxdr
849
fdat1 <- subset(mydat1,
850
851
852
853
854
                select = c("d2peak", "maxdr100k", "gr11", "rgr46"))
fdat1[ , 1] <- formatC(fdat1[,1], format = "f", digits = 0)
fdat1[ , 2] <- formatC(fdat1[,2], format = "f", digits = 2)
fdat1[ , 3] <- formatC(fdat1[,3], format = "f", digits = 2)
fdat1[ , 4] <- formatC(fdat1[,4], format = "f", digits = 2)
855

856
knitr::kable(fdat1,
857
858
             format = "pandoc",
             booktabs = T,
859
860
             row.names = T,
             digits = 0,
861
862
             align = "r",
             caption = capg)
863

864
```
865

866
867
868
869
870
871
\newpage

```{r tabd, include = T}
cap2 <- paste("Indicators of the COVID-19 mortality pattern in Europe:",
             "decay step.")

872
labcol2 <- paste("T10, T20, ..., T95: time lag (d) for a 10\\%, 20\\%, ..., 95\\%",
873
                 "mortality decrease with respect to the peak daily mortality.")
874

875
capd <- paste(cap2, labcol2, "")
876
877
878
879
880

knitr::kable(mydat2,
             format = "pandoc",
             booktabs = T,
             row.names = T,
881
882
             digits = 0,
             caption = capd)
883

884
```
885

886

887
```{r fdecay, include = T, fig.height = 9, fig.cap = cap}
888
889
890
891
cap <- paste("Daily mortality decay after the daily mortality",
             "peak in European countries. The red, dashed, horirontal",
             "line is drawn at daily mortality decay = 50\\%",
             "with respect to the daily mortality peak.")
892

893
onam3 <- onam[onam %in% row.names(mydat1)]
894
decayL <- na.omit(
895
    reshape2::melt(
896
897
                  data = data.frame(name = factor(row.names(mydat2),
                                                  levels = onam3),
898
                                    mydat2),
899
900
                  id.vars = "name"))

901
902
903
decayL$decay <- as.numeric(substr(decayL[,2], 2, 3))
names(decayL) <- c("name", "TXX", "time", "decay")
decayL <- decayL[order(decayL$name, decayL$time), ]
904

905
906
xyplot(decay ~ time | name, data = decayL,
       scales = list(cex = .7, tck=.5),
907
       layout = c(4, 9), 
908
909
910
911
       ylim = c(105, -5),
       as.table = T,
       between = list(x = 1/4, y = 1/4),
       par.strip.text = list(cex = .7),
912
       xlab = list("Days after from the mortality peak", cex = .7),
913
       ylab = list("Mortality decrease (\\%)", cex = .7),
914
915
       panel = function(x, y){
           panel.grid()
916
           panel.abline(h = 50, lty = 4, col = 2)
917
918
919
           panel.xyplot(x, y, type = "l", col = 1)
           lpoints(x, y, pch = 19, col = "white", cex = .9)
           lpoints(x, y, pch = 21, col = 1, fill = "white", cex = .7)
920
921
922
923
       })

```

924

925
926
* Daily mortality rate has decreased to a low level in most countries
  (fig. \@ref(fig:fdecay)). However, it is still at least 50% of the
927
  national peak mortality rate in `r collapse(persist, NA)`. 
928
929


930
```{r pairs, include = T, fig.cap=cap}
931

932
labcol3 <- paste("T50: time lag (d) for a 50\\%",
933
934
935
936
937
938
939
                 "daily mortality decrease with respect to the",
                 "peak in daily mortality.") 

cap <- paste0("Correlation between the considered mortality indicators. ",
              "Above the diagonal: scatter plots augmented with a ",
              "loess-smoothed curve (red line). Under the diagonal: ",
              "absolute value of the linear correlation coefficient. ",
940
              collapse(labcol1, NA), ", ", labcol3, ".")
941
942
943

## put (absolute) correlations on the upper panels,
## with size proportional to the correlations.
944
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
945
946
947
948
949
950
951
952
953
954
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r)
}

955
z1 <- mydat1[, -match("maxdr100k", names(mydat1))]
956
957

z2 <- base::subset(mydat2,
958
                   select =  "T50")
959
960
961
962
963
964
z <- cbind(z1, z2)
miss <- apply(z, 1, function(x) any(is.na(x)))
err <- names(miss)[miss]
onam3 <- onam2[!(onam2 %in% err)]
z <- na.omit(z)

965
par(las = 1)
966

967
pairs(z,
968
      cex.axis = 1,
969
970
      upper.panel = panel.smooth,
      lower.panel = panel.cor,
971
      gap = 0.5, row1attop = FALSE)
972
973

```
974

975
\justify
976

977
We selected variables from the two sets of indicators, (i) based on their
978
completeness regarding the study countries, and (ii) the separate