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.
Wat 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.
De 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!
De 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…
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.
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.
Wat goed is aan het dashboard is dat deze bandbreedte benoemd wordt en beschikbaar is. Wat jammer is dat deze “bandbreedte” niet expliciet gemaakt wordt: voor zover het uit de beschikbaar gestelde data te halen is, 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.
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?
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.
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.
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.
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 datein_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)
#df.ZKH.R <- res_uncertain_si$Rnms <- 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 plotplot2 <- 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 plotsn=2s <- paste0('000',i)s <-substr(s, nchar(s)-n, nchar(s))# sample_posterior_Rsp <- 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 greencolors = 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))}
DIKW Intelligence
Wattbaan 1
3439 ML Nieuwegein
