Dashboard coronavirus

Eerste observaties van een datascientist

Het dashboard van Hugo de Jonge waarin de overheid rapporteert over de stand van zaken met betrekking tot de COVID-19 epidemie is op vrijdag 5 juni  live gegaan. Na het bestuderen van het dashboard vallen een aantal zaken op. Daar gaan we in dit blog wat uitgebreider op in.

blog corona dashboard obeservaties van een datascientist

Open data

Wat mij erg aanspreekt is het feit dat er veel referenties staan naar data, grafieken, code en wetenschappelijke artikelen waarop de cijfers die worden gepresenteerd zijn gebaseerd. Dat is echt positief en een mooi voorbeeld van openheid die je van een overheid verwacht. 

Presentatievorm

dashboard Intensive care opnames per dagDe gekozen presentatievorm is eenvoudig en strak en richt zich op het nu, de tijd ontbreekt dus de trend is onzichtbaar. Wel wordt in ieder figuur een “signaalwaarde” opgenomen, ook wel een Key Performance Indicator genoemd. Het gepresenteerde cijfer is een voortschrijdend gemiddelde over de afgelopen drie dagen, vandaar de decimaal. De stoplicht kleuren geven een indicatie wat veilig en wat risico aangeeft, heel vertrouwd.  


Als we meer uitleg laten zien komt er een kale tijdreeks te voorschijn, thats it ?! De benchmark/signaalwaarde is niet zichtbaar. Ook is er geen trend indicator die aangeeft of we omhoog of omlaag aan het bewegen zijn, hier is ruimte voor verbetering!

verloop intensivecare opnames per dag iover de tijdDe meeste cijfers die worden gepresenteerd zijn tellingen. Het proces van verzamelen van deze data is redelijk overzichtelijk en een driedaags gemiddelde zorgt voor een demping van de variantie in het dashboard. Vanwege dit driedaags gemiddelde zien we op 5 juni ook de cijfers tot en met drie dagen daarvoor (3 juni in dit geval).

Voor het reproductiegetal R hebben we het over een heel ander beest…





Reproductiegetal R

Het reproductiegetal R is een schatting op basis van een heel aantal aannames en modellen van weer andere input parameters (zoals infectiegraad, ziekenhuisopnames, verspreidingskansen etc., etc.). Er zijn een aantal interessante artikelen verschenen diewat meer duiding geven, een goed artikel in de Correspondent van Sanne Blauw en een artikel uit NRC van de hand van Wouter van Loon over de beperkte waarde van R geven achtergrond.dashboard reproductiegetal R

Dit getal is het resultaat van een zogenaamde Monte Carlo simulatie: door heel vaak een bepaalde model berekening uit te voeren waarvan de model parameters die in dat model worden gebruikt zelf niet 1 vaste waarde hebben maar uit een kansverdeling worden getrokken (stochastische variabelen) ontstaat een beeld van de onzekerheden die er in de model uitkomst spelen. De uitkomst van een dergelijk model is een verwachtingswaarde, zeg maar een meest waarschijnlijk gemiddelde van alle mogelijke waardes die voor kunnen komen, gegeven wat we nu weten. Een belangrijk gegeven bij een dergelijke simulatie is de spreiding van die verwachtingswaarde, door het RIVM de bandbreedte genoemd.

De “Bandbreedte”

bandbreedte van reproductiegetal RWat ik goed vind aan het dashboard is dat deze bandbreedte benoemd wordt en beschikbaar is. Wat jammer is dat deze “bandbreedte” niet expliciet gemaakt wordt: voor zover ik het uit de beschikbaar gestelde data kan halen is het een betrouwbaarheidsinterval van 90%. Wat zoveel wil zeggen als, de werkelijke waarde ligt in 9 van de 10 gevallen tussen deze grenzen. Waar ik wat meer moeite mee heb is dat deze bandbreedte in dit geval heel groot is en dat de implicaties daarvan naar mijn mening onderbelicht blijven.

Reproductiegetal R over de loop van de tijdWat zien we hier?

Dit is wederom een tijdreeks, een dikke blauwe lijn en een grijs oppervlak, de zogenaamde bandbreedte. We zien dat de blauw lijn veel eerder ophoudt, de laatste schatting van het getal R is van 15 mei. Dit komt onder andere door de onzekerheid over de incubatietijd, het ontwikkelen van ernstige symptomen, en de rapportage vertraging. Kijken we naar de “bandbreedte” dan zien we dat de R van 15 Mei op 5 juni wordt geschat met een bandbreedte tussen 0.44 en 1.37.De bandbreedte loopt wel door, zou R na 15 mei alweer boven de 1 kunnen liggen? Of is de bandbreedte nu zo groot dat je er eigenlijk geen conclusies meer aan kunt verbinden (de aantallen nieuwe ziekenhuisopnames zijn zo laag dat de bandbreedte enorm toeneemt)? Wat zegt dit over de waarde die we aan het reproductiegetal R moeten hangen? 

Conclusie 

Het dashboard is live, dat is een mooie stap! Er is veel open data en science beschikbaar, een voorbeeld voor overheidsinstanties wat mij betreft. De gerapporteerde tellingen over ziekenhuisopnames en aantallen patiënten zijn eenvoudig maar zonder trend ook beperkt informatief. Mijn grootste probleem zit hem in de presentatie van het reproductiegetal R, deze wordt op dezelfde manier gepresenteerd terwijl er veel meer haken en ogen aan zitten dan je op basis van die presentatie zou vermoeden. In de uitleg staat wel een soort van kanttekening maar daar zal 80% van de Kamerleden geen notie van nemen.

Welke indicator is dan een beter alternatief? 

Kunnen we een alternatief aanbieden aan de minister voor het dashboard? Als we toch iets over dat reproductiegetal willen rapporteren dan opteer ik voor de “Kans dat R groter is dan 1”, dat is naar mijn bescheiden mening een dashboard waardig instrument. Hierbij kun je ook goed de onbetrouwbaarheid weergeven door gewoon de complete feitelijke uitkomst van de simulatie te laten zien, de hele posterior kansverdeling dus. Zo kan iedereen zien hoe die verdeling veranderd in de tijd en zich ook een goed beeld vormen van de spreiding en de onzekerheid van het achterliggende rekenmodel.

DISCLAIMER  : Alle hierna getoonde data is voor demonstratie toepassing van de inhoud van deze blog. Ik ben geen viroloog of anderszins deskundig op de medische inhoud.

De getoonde animatie geeft de posterior kansverdeling weer van de mogelijke waarde die het reproductiegetal R op enig moment kan hebben. De kansverdeling wordt gemaakt door 1000x een trekking te doen uit de geschatte kansverdeling die het reproductiegetal kan hebben gegeven de onzekerheden die er in het spel zijn.

In het begin van de epidemie zien we een heel brede spreiding ver boven de 1, later zakt de R onder de 1 en neemt de spreiding af. Vanaf half mei zien we weer een iets bredere spreiding en neemt de onzekerheid over de waarde van het reproductiegetal weer iets toe.

Animatie reproductiegetal R

De code

In het kader van open source beleid ook hier natuurlijk een stukje R code waarmee je de bovenstaande animatie kunt maken.

# epiestim package demo
# https://cran.r-project.org/web/packages/EpiEstim/vignettes/demo.html

library(EpiEstim)
library(ggplot2)
# op basis van ziekenhuis opnames
# https://www.rivm.nl/coronavirus-covid-19/grafieken
# download het csv bestand van de ziekenhuisopnames

library(readr)

in_ziekenhuis_opgenomen_patiënten <- read_delim("data/in-ziekenhuis-opgenomen-patiënten.csv",
";", escape_double = FALSE, trim_ws = TRUE)
names(in_ziekenhuis_opgenomen_patiënten) <- c('datum','nieuw','bekend')
# convert to date
in_ziekenhuis_opgenomen_patiënten$datum <- as.Date(paste(in_ziekenhuis_opgenomen_patiënten$datum,'2020',' '),format='%d %b %Y')

library(incidence)

df.ZKH.incd <- as.incidence(in_ziekenhuis_opgenomen_patiënten$bekend, dates = in_ziekenhuis_opgenomen_patiënten$datum)

plot(df.ZKH.incd)

## UNCERTAIN SI

## waardes afkomstig uit
# https://royalsocietypublishing.org/doi/full/10.1098/rspb.2006.3754

config <- make_config(list(mean_si = 4.8, std_mean_si = 0.6,
min_mean_si = 3.8, max_mean_si = 6.1,
std_si = 2.3, std_std_si = 0.5,
min_std_si = 1.6, max_std_si = 3.5))

res_uncertain_si <- estimate_R(df.ZKH.incd,
method = "uncertain_si",
config = config)
#> Default config will estimate R on weekly sliding windows.
#> To change this change the t_start and t_end arguments.
plot(res_uncertain_si, legend = FALSE)

resultaat epiestim fit reproductiegetal R op Nederlandse ziekenhuisopnames

#
df.ZKH.R <- res_uncertain_si$R
nms <- c("t_start", "t_end", "Mean_R", "Std_R", "Q0_025","Q0_05","Q0_25","Median_R","Q0_75","Q0_95","Q0_975","dates")
df.ZKH.R$dates <- res_uncertain_si$dates[8:100]
names(df.ZKH.R) <- nms

ggplot(df.ZKH.R, aes(x=dates, y=Median_R)) +
geom_line(col='blue',size=1.2) +
geom_ribbon(data=df.ZKH.R,aes(ymin=Q0_025,ymax=Q0_975),alpha=0.3) +
scale_y_continuous(limits = c(-0, 6)) +
geom_hline(yintercept=1, linetype="dashed", color = "red")

# a way to get separate plots for each plot
plot2 <- function(theplot, name, ...) {
name <- paste0(name, ".png")
png(filename=name, width = 1200, height = 627)
print(theplot)
dev.off()
} #plotting function

for (i in 1:length(res_uncertain_si$dates)) {
# maak voor iedere datum een plot van de posterior van R

# gerommel met de naam van de plots
n=2
s <- paste0('000',i)
s <-substr(s, nchar(s)-n, nchar(s))

# sample_posterior_R
sp <- sample_posterior_R(res_uncertain_si, n=1000, win=i)
df.post <- as.data.frame(sp)
names(df.post) <- c('posterior.R')

# color if true then red else green
colors = c('TRUE'='red', 'FALSE'='green')

p <- ggplot(df.post, aes(x=posterior.R)) + geom_histogram(binwidth=.01, aes(fill = posterior.R > 1.0)) +
geom_vline(aes(xintercept=mean(posterior.R)),
color="blue", linetype="dashed") +
geom_vline(aes(xintercept=1.0),
color="red", linetype="dashed") +
scale_fill_manual(values = colors) +
annotate(geom="text", x=3.5, y=15, label=paste("R verw : ", round(mean(df.post$posterior.R),2))) +
annotate(geom="text", x=3.5, y=5, label=paste("Datum : ",res_uncertain_si$dates[i])) +
guides(fill=FALSE) +
xlim(c(0,4)) +
ylim(c(0,200)) +
ggtitle('Reproductiegetal R', subtitle = '(Schatting o.b.v.Posterior Kansverdeling n=1000)') +
theme_bw() +
theme(plot.title = element_text(size = 40, face = "bold"),
plot.subtitle = element_text(size = 20))

plot2(p, name=paste0("TEST", s))

}


Voor de duidelijkheid, deze analyse probeert de analyses van de RIVM na te doen op basis van de verstrekte informatie. Deze analyse is niet de werkelijke analyse zoals door het RIVM gedaan en is ook niet als zodanig bedoeld. We willen hier laten zien hoe de onzekerheid met betrekking tot het reproductiegetal tot stand komt en hoe je, door die onzekerheid duidelijker te laten zien, een alternatieve weergave kunt maken voor in het dashboard van de overheid.

Natuurlijk is alle commentaar welkom, met de kracht van open source data analyse kunnen we een samen een aantal mooie suggesties voor het coronavirus  dashboard ontwikkelen!