Commit 8a3b19c7 authored by Facundo Muñoz's avatar Facundo Muñoz ®️
Browse files

Covid response report.

parent 1a055e3d
---
title: |
|
| Available data to assess the effect of lockdown on COVID-19 transmission in Europe and its neighborhood
subtitle: A [MOOD](https://mood-h2020.eu/) output, released on `r sdt`
author: "`r au`"
date: "Last updated: `r date()`"
tables: yes
bibliography: d:/biblio/covid.bib
csl: d:/biblio/csl/plos-one.csl
output:
rmdformats::readthedown:
toc_depth: 3
mathjax: "default"
lightbox: true
gallery: true
use_bookdown: true
number_sections: true
---
```{js logo-js, echo=FALSE}
$(document).ready(function() {
$('#header').parent().prepend('<div id=\"logo\"><img src=\"logoMood_vert.png\" style=\"position:absolute; top:0; left:85; padding:0px; height:60px\"</div>');
$('#header').css('margin-leftt', '85px')
});
```
```{r settings, include = F}
knit_hooks$set(optipng = hook_optipng)
knit_hooks$set(pngquant = hook_pngquant)
knitr::knit_hooks$set(crop = knitr::hook_pdfcrop)
knitr::opts_chunk$set(cache = F, dev = "png", type = "cairo", dpi = 300, fig.retina = 2,
crop = T, include = F, echo = F, warning = F, message = F,
comment = NA, results = "markup")
Sys.setlocale("LC_TIME", "English")
sdt <- format(Sys.Date(), "%d %b %Y")
wgs84 <- "+proj=longlat +ellps=WGS84 +datum=WGS84"
merc <- "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
library(RODBC)
## install.packages('countrycode')
library(countrycode)
library(rjson)
library(grid)
library(lattice)
library(ade4)
##help(package="adegraphics")
library(ggplot2)
## install.packages("ggthemes")
library(ggthemes)
## install.packages("ggrepel")
library(ggrepel)
library(latticeExtra)
library(sp)
library(raster)
library(cshapes)
library(RColorBrewer)
library(INLA)
library(ade4)
eust <- paste(
"MOOD has received funding from the",
"European Union's Horizon 2020 research",
"and innovation programme under grant",
"agreement No 874850")
eu28 <- c("Austria", "Belgium", "Bulgaria", "Croatia",
"Cyprus", "Czech Republic", "Denmark", "Estonia",
"Finland", "France", "Germany", "Greece", "Hungary",
"Ireland", "Italy", "Latvia", "Lithuania", "Luxembourg",
"Malta", "Netherlands", "Poland", "Portugal",
"Romania", "Slovakia", "Slovenia", "Spain", "Sweden",
"United Kingdom")
efta <- c("Iceland", "Liechtenstein", "Norway", "Switzerland")
candi <- c("Albania", "Montenegro", "North Macedonia",
"Serbia", "Turkey")
noneu <- c("Bosnia and Herzegovina", "Kosovo")
enp <- c("Algeria", "Armenia", "Azerbaijan", "Belarus", "Egypt",
"Georgia", "Israel", "Jordan", "Lebanon", "Libya",
"Moldova", "Morocco", "Syria", "Tunisia", "Ukraine")
ctr <- sort(c(eu28, efta, candi, noneu, enp))
```
# Oxford COVID-19 Government Response Tracker
## Goal{-}
Quoting the [web site](https://covidtracker.bsg.ox.ac.uk/),
"_We aim to track and compare worldwide government responses to the
coronavirus rigorously and consistently. Systematic information on
which measures governments take, and when, can help us understand
the responses in a consistent way, aiding efforts to fight the
pandemic._" [@Hale2020]
## Data{-}
"_We collate publicly available information on a number of
indicators of government response. The first seven indicators
(S1-S7), taking policies such as school closures, travel bans, etc,
are recorded on an ordinal scale; the remainder are financial
indicators such as fiscal or monetary measures._"
"_The COVID-19 Government Response Stringency Index is an additive
score of the seven indicators (S1-S7) measured on an
ordinal scale, rescaled to vary from 0 to 100. This measure is for
comparative purposes only._"
```{r getdata1, include=F}
levXX <- c("Belgium", "Spain", "Italy", "United Kingdom",
"France", "Sweden", "Netherlands", "Ireland",
"Switzerland", "Luxembourg", "Portugal", "Denmark",
"Germany", "Austria")
## donload the data
jsn <- "https://covidtrackerapi.bsg.ox.ac.uk/api/v2/stringency/date-range/2020-01-01/2020-05-30"
## get the data into lists
confi <- fromJSON(file = jsn)
## function to get the names from th code
foo <- function(x){
ifelse(x == "RKS",
"Kosovo",
countrycode(sourcevar = x,
origin = "iso3c",
destination = "country.name"))
}
## USE IT
cgrtName <- data.frame(
code = confi$countries,
name = sapply(confi$countries, function(x) foo(x)),
stringsAsFactors = F)
## clean the names
idx <- grep(pattern = "Czech", x = cgrtName$name)
cgrtName$name[idx] <- "Czech Republic"
## filter tarfet countrues
code <- subset(cgrtName, name %in% ctr)
## build condition to check data at the lowest level of nesting:
## needed because of some mistakes in the json file
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 &",
"x$country_code %in% code$code")
## get the data from the complicated structure of the data: list of
## lists of lists
cgrtData <- confi$data
Liste <- lapply(cgrtData,
function(x)
do.call("rbind",
lapply(x,
function(x){
if(eval(parse(text = cond)))
data.frame(x)
else
NULL
})))
Confi <- do.call("rbind", Liste)
Confi$date <- as.Date(substr(rownames(Confi), 1, 10))
rownames(Confi) <- seq(nrow(Confi))
## get the country names in the data.frzme
lock <- merge(Confi, code, by.x = "country_code", by.y = "code")
## sort by country and date
lock <- lock[order(lock$name, lock$date), ]
```
## Visualization{-}
Preliminary examination revealed a very strong correlation of the four
different items.We only used `stringency_actual`. The dataset is shown
on fig. \@ref(fig:show).
```{r show, include = T, fig.width=7, fig.height = 9, fig.cap = cap}
cap <- paste("Changes in stringency index used to monitor",
"government response along the 2020",
"COVID-19 outbreak in Europe.")
cex <- .7
lock$name <- reorder(lock$name, lock$stringency)
xyplot(stringency_actual ~ date | name,
data = lock, subset = name %in% levXX,
scales = list(cex = cex, tck = .5),
xlab = list("Time (d)", cex = cex),
ylab = list("Stringency index", cex = cex),
layout = c(2, length(levXX)/2),
between = list(x = 1/4, y = 1/4),
par.strip.tec = list(cex = cex),
type = "l") +
layer(panel.grid(), under = T)
```
## Usage{-}
We propose to categorize the stringency $s_i$ of government's response
in a given country $i$ ($i \in {1, \ldots, I}$) according to different
features of the response pattern:
* Earliness $E_i$, and relative earliness $RE_i$:
+ Find the date $t_0$ of the first rise in stringency index
$s_{i, t_0}$ across the $I$ countries taken as the
reference,
+ For each country $i$, calculate the time lag between the time of
first increase in stringency index $t_{0, i}$ and $t_0$:
earliness $E_i = t_{0, i} - t_0$, and relative earliness
$$ RE_i = \frac{E_i - min (E_i)}{max(E_i) - max(E_i)} $$
* Lateness $L_i$, and relative lateness $RL_i$::
+ Find the date $t_{i, max}$ when the maximum stringency index
$s_{i, t_{max, i}}$ was reached across the $I$ countries,
+ Compute the lateness $L_i = t_{max, i} - min(t_{max, i})$, and the
relative lateness:
$$ RL_i = \frac{L_i}{max(t_{max, i}) - min(t_{max, i})} $$
* Progressiveness $P_{i}$ is the rate of change in stringency from
$t_{0, i}$ to $t_{max, i}$:
$$ P_{i} = \frac{s_{i, t_{max, i}} - s_{i, t_{0, i}}}{t_{max,
i} - t_{0, i}} $$
* Final stringency index $FS_i = s_{i, t_{max, i}}$, and relative
final stringency index:
$$ RFS_i = \frac{s_{i, t_{max, i}} - min(s_{i, t_{max, i}})}{max(s_{i, t_{max, i}}) - min(s_{i, t_{max, i}})} $$
To find the clusters of countries with similar response patterns, we
run a principal components analysis on the four indicators defined
above, followed by a hierarchical ascending classification on the
scoring factors corresponding to the principal components.
## Implementation{-}
```{r indic, include = F}
z <- subset(lock,
subset = name %in% levXX,
select = c("name", "date", "stringency_actual"))
names(z)[3] <- "strg"
MyList <- by(
data = z,
INDICES = list(name = z$name),
FUN = function(x){
## earliness
t0 <- min(z$date[z$strg > 0],
na.rm = T)
t0i <- min(x$date[x$strg > 0], na.rm = T)
Ei <- t0i - t0
## lateness
tmaxi <- min(x$date[which(x$strg == max(x$strg, na.rm = T))], na.rm = T)
## Progressiveness
Deltat <- as.numeric(tmaxi - t0i)
Pi <- (x$strg[x$date == tmaxi] - x$strg[x$date == t0i]) / Deltat
data.frame(name = unique(x$name),
t0 = t0,
Ei = as.numeric(Ei),
tmaxi = tmaxi,
Pi = Pi,
Smaxi = x$strg[x$date == tmaxi])
})
tmp3 <- do.call("rbind", MyList)
tmp3$REi <- (tmp3$Ei - min(tmp3$Ei)) / as.numeric(diff(range(tmp3$Ei)))
tmp3$Li <- tmp3$tmaxi - min(tmp3$tmaxi)
tmp3$RLi <- as.numeric(tmp3$Li) / as.numeric(diff(range(tmp3$Li)))
tmp3$RFSi <- with(tmp3, (Smaxi - min(Smaxi)) / (max(Smaxi) - min(Smaxi)))
resu <- subset(tmp3,
name != "Sweden",
select = c("name", "REi", "RLi", "Pi", "RFSi"))
```
```{r pca, include = T, fig.width = 7, fig.height = 7, fig.cap = cap}
cap <- paste("Categorization of 13 European countries according to",
"the stringency indicators of their government'e response",
"along the 2020 COVID-19 outbreak. A: eigen values (%) of",
"the principal component analysis (PCA) of the indicator",
"table. B: correlation circle of the indicators included",
"in the PCA. C: dendrogram of the hierarchical ascending",
"classification on the PCA scores.")
acp <- dudi.pca(resu[,-1], scannf = F, nf = 4)
sli <- acp$li
d <- dist(sli)
h <- hclust(d, method = "ward.D")
layout(mat = matrix(c(1, 2, 3, 3), ncol =2, byrow = T))
## aigen values
barplot(100 * acp$eig / sum(acp$eig), las = 1,
ylab = "Variance (%)")
title(main = "A", adj = 0)
## correlation circle
s.corcircle(acp$co)
title(main = "B", adj = 0)
## dendrogram
plot(h, xlab="", hang = .1, main = "", sub = "", las = 1)
title(main = "C", adj = 0)
k <- 3
resu$k <- cutree(h, k = k)
cols <- brewer.pal(n = k, "Dark2")
Cols <- resu$k
tab <- table(resu$k, h$order)
idx <- apply(tab, 1, function(x) paste(x, collapse=""))
tab <- tab[order(idx, decreasing = T), ]
ord <- as.numeric(rownames(tab))
rect.hclust(h, k = k, which = ord, border = cols)
legend("topright", bg= "white", cex = .7,
legend = paste("class", 1:3),
fill = "white", border = cols[c(3,1,2)])
out <- cbind(resu, acp$li)
write.table(x = out,
file = "PCAscoresOCGRT.csv",
sep = ";",
row.names = T)
```
* The principal component analysis (PCA) highlights the preponderance
of the first two eigen values, representing 90% of the overall
variance (fig. \@ref(fig:pca)A).
* The correlation circle (fig. \@ref(fig:pca)B) shows relative
lateness and final stringency are positively correlated, as well as
relative earliness, and progressiveness. However, these two pairs of
variables are negatively correlated, and this dipole makes the first
PCA axis, with high values of final stringency on the left.
* The second axis discriminates countries accorging to the intensity of
these variables, with strong values on the bottom.
* Following country classification using the PCA scores
(fig. \@ref(fig:pca)C), we selected a partition with three classes.
```{r plotAcp, include = T, fig.width=7, fig.height=7, fig.cap = cap}
cap <- paste("Projection of country categories",
"on the plane defined by the first and",
"second axes of the principal component",
"analysis of the stringency indicators of",
"government'e response along the 2020",
"COVID-19 outbreak in Europe.")
xl <- extendrange(range(acp$li[,1]))
yl <- extendrange(range(acp$li[,2]))
acp$li$Class <- factor(resu$k)
acp$li$name <- rownames(acp$li)
Groups <- cols[resu$k]
ggplot(acp$li, aes(Axis1, Axis2)) +
coord_fixed() +
geom_point(aes(shape = "19", colour = Groups), show.legend = F) +
scale_colour_manual(labels = as.character(1:k),
values = cols) +
theme_bw() +
geom_hline(lty = 4, yintercept = 0, show.legend = F) +
geom_vline(lty = 4, xintercept = 0, show.legend = F) +
xlim(xl[1], xl[2]) +
ylim(yl[1], yl[2]) +
ggtitle(label = paste("Projection of countries on the plane defined",
"by 1st and 2nd PCA axes")) +
labs(y = "Country scores on 1st axis",
x = "Country scores on 2nd axis") +
geom_text_repel(aes(label = name, colour = Groups))
```
The scatter plot of country scores (fig. \@ref(fig:plotAcp)) shows an
ordering of these countries along the first PCA axis. According to this
metric, France had the highes stringency, and Luxembourg the lowest.
# References{-}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment