COVID-19 Weersverwachting

Ter ondersteuning van het corona dashboard van de rijksoverheid

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. Mijn eerdere blog gaat in op het feit dat het op zijn minst vreemd is dat het reproductiegetal op twee cijfers achter de komma wordt gepresenteerd naast rechte tellingen terwijl de (on)betrouwbaarheid van dat getal enorm groot is, misschien wel te groot om daar beleid op te bepalen. Om je een indruk te geven van die onzekerheid  heb ik dat eerste blog geschreven 

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.

Rapportage onzekerheid

Een van de zaken die aandacht verdient is de rapportage onzekerheid, de rapportage van de ziekenhuisopnames ijlt nogal na. Er druppelen nog dagen later opnames binnen en zeker als er weinig nieuwe opnames zijn heeft dat een enorme invloed op het reproductiegetal. Wat opvalt in het dashboard van de overheid is dat de betrouwbaarheidsgrenzen (ook wel bandbreedte genoemd door het RIVM) wel wordt gepresenteerd.

Hoe komen ze daar aan als ze de R zelf niet hebben?

Simulatie

Een mogelijke aanpak kan zijn het aantal ziekenhuis opnames te simuleren. Je zou kunnen zeggen dat het exacte aantal ziekenhuis opnames dat we daadwerkelijk hebben gezien voor een deel toevallig is. Jan of Piet die al een paar dagen klachten heeft kan denken " Weet je wat ik ga die COVID test morgen wel even doen", dus voor je het weet het je 1 of 2 ziekenhuisopnames meer of minder. Dat kunnen we natuurlijk prima in de computer naspelen, een simulatie.

Als je de ziekenhuisopnames van de afgelopen 14 dagen als uitgangspunt neemt kun je daar eenvoudig een negatief binomiale distributie op fitten. Dit type distributie kun je gebruiken om voor 1 dag het aantal nieuwe ziekenhuis opnames te schatten, als je dat voor 14 dagen doet heb je dus een mogelijke versie van de werkelijkheid/toekomst. Die werkelijkheid kun je achter de cijfers plakken die we al gezien hebben en door het RIVM worden vrijgegeven en dan kun je het reproductie getal sommetje weer maken.

Maar dan heb je nog maar 1x het mogelijke verloop van het aantal ziekenhuisopnames gebruikt, we kunnen (en moeten) dat veel vaker doen om een beeld te krijgen van de onzekerheid die er in de toekomst verscholen ligt. Laten we zeggen dat we dit 25 x doen. Dan ziet de ontwikkeling van het reproductiegetal er op basis van die 25 simulatieruns er zo uit:


Het is op het moment van het maken van deze figuur 15 juni 2020, de grijze verticale stippellijn geeft vandaag aan. Voor iedere simulatie van de ziekenhuisopnames rekenen we het verloop van R uit (en haar bandbreedte), op vandaag zetten we een dikke rode punt, we doen dat 25 x, de "25 random runs".

Deze 25 random runs kunnen we dan verder gebruiken om een schatting te doen van de posterior R (zoals ook al in mijn eerdere blog beschreven. Als we de posteriors van alle simulaties samen nemen voor iedere dag kunnen we voor iedere dag een distributie laten zien van mogelijke waarde die R op die dag zou kunnen hebben.

De 14 daagse COVID-19 weersverwachting

Dus op basis van de gesimuleerde data uit de voorgaande runs komen we dan tot de volgend weersverwachting (inclusief de enorme bandbreedtes die nu eenmaal bij een weersverwachting horen...)

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.

Dus we hebben naar aanloop van het komende weekeinde kans op een onstuimige R groter dan 1, na het weekeinde weer iets rustiger weer en de R zakt weer richting de 1.

(Omdat het aantal ziekenhuisopnames de afgelopen week zeer klein was en de R een verhoudingsgetal voor het aantal ziekenhuisopnames nu in vergelijking met een week gelden is,is  het niet gek dat met slecht een aantal extra ziekenhuisopnames de R weer door de 1 schiet.)

De code

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

# tweede blog over simualtie ziekenhuiopnames

library(EpiEstim)
library(ggplot2)
library(incidence)
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')
# orginele data behouden ter controle
df.izop <- in_ziekenhuis_opgenomen_patiënten
df.izop$alles <- df.izop$nieuw + df.izop$bekend
# incidence object heeft epistim nodig
df.ZKH.incd <- as.incidence(df.izop$bekend, dates = df.izop$datum)

# config settings voor epistim
## 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_0 <- estimate_R(df.ZKH.incd,
method = "uncertain_si",
config = config)

r0 <- dim(res_uncertain_si_0$R)[1]

plot(res_uncertain_si_0, 'R')


#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

# generate new cases from distribution
# fit kans op nieuwe ziekenhuis opnames

# neem laatste 14 dagen
r <- nrow(in_ziekenhuis_opgenomen_patiënten)
df.14 <- in_ziekenhuis_opgenomen_patiënten[(r-14):r,]
df.14$alles <- df.14$nieuw + df.14$bekend
hist(df.14$alles, breaks=150)

nu <- in_ziekenhuis_opgenomen_patiënten[r,]$datum

# fitdistplus
library(fitdistrplus)

# fit the negative binomial distribution op de laatste 14 dagen van de ziekenhuisopnames
fit <- fitdist(df.14$alles, "nbinom")
plot(fit)

# plaatje checken
f <- rnbinom(n=1000, size=fit$estimate[1], mu=fit$estimate[2])
hist(f, breaks=25)

# simulatie 2 weken vooruit
n <- 14 #dagen

# vanaf laatste datum n x doortellen
d <- seq(max(df.izop$datum+1), by='day', length.out = n)

# aantal random runs
rep = 25

# lijst met simulatie resultaten
sim = list()
sim_complete = list()

# Start simulatie runs
for (i in 1:rep) {
# sample from fitted dist
s <- rnbinom(n=n, size=fit$estimate[1], mu=fit$estimate[2])

# sample 14 dagen rapportages
df.s <- data.frame(d,rep(0,n),s,s)
names(df.s) <- names(df.izop)

df.izop.aug <- rbind(df.izop,df.s)
df.izop.aug.incd <- as.incidence(df.izop.aug$alles, dates = df.izop.aug$datum)

res_uncertain_si <- estimate_R(df.izop.aug.incd,
method = "uncertain_si",
config = config)

#plot(res_uncertain_si, legend = FALSE)
#plot(res_uncertain_si, 'R')

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")
l <- length(res_uncertain_si$dates)
df.ZKH.R$dates <- res_uncertain_si$dates[1:(l-7)] # laatste week geen schatting
names(df.ZKH.R) <- nms

# append result to sim list
sim[[i]] <- df.ZKH.R
sim_complete[[i]] <- res_uncertain_si
}

df.ZKH.R <- sim[[1]]

# opzet van de plot
p <- ggplot(df.ZKH.R, aes(x=dates, y=Mean_R)) +
#geom_line() +
geom_hline(yintercept=1, color="red", linetype="dashed") +
geom_vline(xintercept=as.Date(nu), color="grey", linetype="dashed") +
annotate(geom="text", x=as.Date('2020-05-01'), y=1.5, label=paste(" ",nu), colour = 'grey80', size=20, alpha = 0.5) +
geom_line(data = df.ZKH.R[1:r0,], aes(x = dates, y = Mean_R), color='blue') +
geom_point(data = df.ZKH.R[r0,], aes(x=dates,y=Mean_R),size=2,color='blue') +
geom_ribbon(data=df.ZKH.R[1:r0,],aes(ymin=Q0_05,ymax=Q0_95),alpha=0.2) +
scale_fill_manual(values = colors) + ylim(c(0,3.5)) +
guides(fill=FALSE) +
ggtitle('Reproductiegetal R COVID-19',
subtitle = 'Op basis van simulatie ziekenhuisopnames') +
#geom_label(label=paste("Random run ",i)) +
theme_bw() +
theme(plot.title = element_text(size = 35, face = "bold"),
plot.subtitle = element_text(size = 20))

# invullen run na run
for (i in 1:rep) {
temp <- sim[[i]]

p <- p +
geom_line(data = temp[r0:(r0+n),], aes(x = dates, y = Mean_R), color='red') +
geom_ribbon(data=temp[r0:(r0+n),],aes(ymin=Q0_05,ymax=Q0_95),alpha=0.1) +
geom_point(data = temp, aes(x=dates[r],y=Mean_R[r]),size=2,color='red') +
annotate(geom="label", x=as.Date('2020-05-01'), y=2.5, label=paste("Random run ",i),
colour = 'grey', size=20, alpha = 0.8, fill = 'white')
# plaatje maken voor animatie
plot2(p, name=paste0("TEST_RAP", i))
# rode lijn overschrijven met grijze
p <- p +
geom_line(data = temp[r0:(r0+n),], aes(x = dates, y = Mean_R), color='grey80')

}

# schatting Reproductiegetal van 14 dagen forecast
for (w in (r-7):(r+7)) {

res_unc <- sim_complete[[1]]
# gerommel met de naam van de plots
n=2
s <- paste0('000',w)
s <-substr(s, nchar(s)-n, nchar(s))

# sample_posterior_R
sp <- sample_posterior_R(res_unc, n=1000, win=w)

for (i in 2:rep) {
res_unc <- sim_complete[[i]]
# sample_posterior_R
sp_i <- sample_posterior_R(res_unc, n=1000, win=w)
sp <- c(sp,sp_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')

p2 <- 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") +
annotate(geom="text", x=2, y=600, label=paste("R = ", round(mean(df.post$posterior.R),2)), colour = 'grey', size=40, alpha = 0.4) +
annotate(geom="text", x=2, y=300, label=paste(" ",res_uncertain_si$dates[w]), colour = 'grey', size=40, alpha = 0.4) +
scale_fill_manual(values = colors) +
guides(fill=FALSE) +
xlim(c(0,4)) +
ylim(c(0,1000)) +
ggtitle('Reproductiegetal R COVID-19',
subtitle = '(Schatting o.b.v. Nederlandse ziekenhuisopnames, Posterior Kansverdeling n=1000)') +
theme_bw() +
theme(plot.title = element_text(size = 40, face = "bold"),
plot.subtitle = element_text(size = 20))

plot2(p2, name=paste0("FORE", 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!

BLOG OVERZICHT

Blogs

30 jaar intelligence: Nieuwe uitdagingen om met data waarde toe te voegen — door Marco van den Doel — last modified 25-01-2022
Van oude computerterminal naar smartphone
De fascinerende wereld van testen — door Marco van den Doel — last modified 30-12-2021
Verwacht het onverwachte
Verzekeraar creëert meerwaarde met slimme data hub — door Marco van den Doel — last modified 20-12-2021
Klant maakt met betere voorspellingen met data
Op naar een mooi data gedreven 2022! — door Marco van den Doel — last modified 14-12-2021
Data gedreven organisatie dient blijvend te worden gevoed
Machine Learning: De gereedschapskist van de data scientist — door Marco van den Doel — last modified 21-12-2021
Machine Learning algoritmes zijn de gereedschappen voor een data scientist
Met data de wind in de zeilen — door Marco van den Doel — last modified 28-01-2022
Met data management kiest u de juiste koers
Hoe ethisch is Facebook? — door Marco van den Doel — last modified 05-11-2021
Is regulering en wetgeving voor AI nodig?
Boekbespreking: Data Teams van Jesse Anderson — door Marco van den Doel — last modified 02-11-2021
Voor succesvolle big data projecten zijn drie teams nodig
Smells like AI — door Marco van den Doel — last modified 01-11-2021
Artificial Intelligence creëert nieuwe muziek
De waarde van data voor het MKB — door Marco van den Doel — last modified 09-12-2021
Bluemine ontzorgt MKB door data beheer
Data & AI: Kans of bedreiging? — door Marco van den Doel — last modified 21-12-2021
Waarde creëren met data en AI zorgt voor nieuwe business mogelijkheden
Data gedrevenheid is proces van lange adem — door Marco van den Doel — last modified 02-11-2021
Data is een ingrediënt dat zorgt voor meerwaarde op lange termijn
Hoe data leidt tot de optimalisatie van de customer journey — door Marco van den Doel — last modified 02-11-2021
Ondersteun uw customer journey met data strategie
Wat is data engineering? — door Marco van den Doel — last modified 02-11-2021
Hoe word je een data engineer?
De fasen om te transformeren naar een data gedreven organisatie — door Marco van den Doel — last modified 02-11-2021
Welke vier fasen doorloopt een organisatie naar data gedrevenheid?
Data gedreven organisaties hebben grotere kans om te overleven — door Marco van den Doel — last modified 02-11-2021
Transformeren naar een data gedreven organisatie kost tijd
Data gedreven logistiek onderhoud voorkomt uitval — door Marco van den Doel — last modified 20-12-2021
Operationele en logistieke kosten lager door gebruik van data
Er zijn meer logistieke wegen die naar Rome leiden — door Marco van den Doel — last modified 28-01-2022
Duurzame innovatieve logistieke oplossing op basis van data science
Welke sandwich mogen wij voor u bereiden? — door Marco van den Doel — last modified 02-11-2021
Data gedrevenheid is als een goede en juist belegde sandwich
Met data bijdragen aan een betere wereld — door Marco van den Doel — last modified 02-11-2021
DIKW is partner van Sensing Clues
In het verleden behaalde resultaten... — door Marco van den Doel — last modified 02-11-2021
Data zorgt voor betere resultaten in de toekomst
Blokkade Ever Given geeft noodzaak betere data science aan — door Marco van den Doel — last modified 20-12-2021
Containerschip blokkeert Suezkanaal
Van voor naar achteren en van links naar rechts in de logistieke keten — door Marco van den Doel — last modified 24-01-2022
Het verminderen van opslagkosten en verplaatsingen van het aantal containers
De zeven pilaren van DataOps — door Marco van den Doel — last modified 02-11-2021
DataOps wordt gedefinieerd door zeven hoofdkenmerken
Kijk verder dan je dashboard — door Marco van den Doel — last modified 23-11-2021
Met de DIKW Analytical Roadmap kijk je verder!
Data discovery tools — door Marco van den Doel — last modified 28-01-2022
Hoe zorgen data gedreven organisaties er voor dat data snel gevonden wordt en dat nieuwe medewerkers snel productief zijn?
Van BICC naar DACoE deel 3 — door Marco van den Doel — last modified 18-01-2022
Van BICC naar Data & Analytics Center of Excellence deel 3:  Van ambitie naar realiteit
COVID-19 Oktober forecast : Het kan vriezen het kan dooien — door Marco van den Doel — last modified 23-11-2021
Ter ondersteuning van het corona dashboard van de rijksoverheid
Koning TOTO: Sjaak vs Bayes — door Marco van den Doel — last modified 25-01-2022
Definitieve uitslag voetbal eredivisie op basis van een wiskundig model, Bayesiaanse statistiek, en een kleine Monte Carlo simulatie

Data Science recente blogs

Machine Learning: De gereedschapskist van de data scientist — door Marco van den Doel — last modified 21-12-2021
Machine Learning algoritmes zijn de gereedschappen voor een data scientist
Hoe ethisch is Facebook? — door Marco van den Doel — last modified 05-11-2021
Is regulering en wetgeving voor AI nodig?
Smells like AI — door Marco van den Doel — last modified 01-11-2021
Artificial Intelligence creëert nieuwe muziek

Data Science Nieuws & Evenementen

DIKW in top 50 beste data science bedrijven — door Marco van den Doel — last modified 22-10-2021
DIKW is één van snelst groeiende bedrijven volgens MKB Data Science rapport
AI Hub Midden Nederland gelanceerd! — door Marco van den Doel — last modified 11-11-2021
DIKW is partner van de AI Hub Midden Nederland en ondersteund en helpt het MKB in de regio
AEDES innovatie boost datascience powered by DIKW — door Marco van den Doel — last modified 11-11-2021
De innovatie boost van AEDES is binnen gehaald door de werkgroep Big data