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

Bayesiaanse Statistiek — door Marc Jacobs — last modified 25-07-2022
Wiskundig raamwerk voor ouderwets leren
R vs Python — door Nick van de Venn — last modified 05-07-2022
Samenwerken is de sleutel
Je fietsroutes eenvoudig in kaart brengen... — door Nick van de Venn — last modified 05-05-2022
Marcel-Jan doet het eenvoudig met behulp van Python.
Het DIKW model — door marco — last modified 24-02-2022
In vier stappen waarde creëren met data
DIKW Academy: Waar theorie en praktijk samen komen — door marco — last modified 21-02-2022
DIKW docenten delen hun expertise uit de praktijk
Laat data voor u renderen! — door marco — last modified 08-02-2022
Bluemine: Analytics as a service
30 jaar intelligence: Nieuwe uitdagingen om met data waarde toe te voegen — door marco — last modified 08-02-2022
Van oude computerterminal naar smartphone
De fascinerende wereld van testen — door marco — last modified 30-12-2021
Verwacht het onverwachte
Verzekeraar creëert meerwaarde met slimme data hub — door marco — last modified 20-12-2021
Klant maakt met betere voorspellingen met data
Op naar een mooi data gedreven 2022! — door marco — last modified 14-12-2021
Data gedreven organisatie dient blijvend te worden gevoed
Machine Learning: De gereedschapskist van de data scientist — door marco — last modified 21-12-2021
Machine Learning algoritmes zijn de gereedschappen voor een data scientist
Met data de wind in de zeilen — door marco — last modified 28-01-2022
Met data management kiest u de juiste koers
Hoe ethisch is Facebook? — door marco — last modified 05-11-2021
Is regulering en wetgeving voor AI nodig?
Boekbespreking: Data Teams van Jesse Anderson — door marco — last modified 02-11-2021
Voor succesvolle big data projecten zijn drie teams nodig
Smells like AI — door marco — last modified 01-11-2021
Artificial Intelligence creëert nieuwe muziek
De waarde van data voor het MKB — door marco — last modified 09-12-2021
Bluemine ontzorgt MKB door data beheer
Data & AI: Kans of bedreiging? — door marco — 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 — 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 — last modified 02-11-2021
Ondersteun uw customer journey met data strategie
Wat is data engineering? — door marco — last modified 03-02-2022
Hoe word je een data engineer?
De fasen om te transformeren naar een data gedreven organisatie — door marco — 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 — last modified 02-11-2021
Transformeren naar een data gedreven organisatie kost tijd
Data gedreven logistiek onderhoud voorkomt uitval — door marco — 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 — last modified 28-01-2022
Duurzame innovatieve logistieke oplossing op basis van data science
Welke sandwich mogen wij voor u bereiden? — door marco — last modified 22-02-2022
Data gedrevenheid is als een goede en juist belegde sandwich
Met data bijdragen aan een betere wereld — door marco — last modified 18-02-2022
DIKW is partner van Sensing Clues
In het verleden behaalde resultaten... — door marco — last modified 02-11-2021
Data zorgt voor betere resultaten in de toekomst
Blokkade Ever Given geeft noodzaak betere data science aan — door marco — last modified 11-02-2022
Containerschip blokkeert Suezkanaal
Van voor naar achteren en van links naar rechts in de logistieke keten — door marco — last modified 24-01-2022
Het verminderen van opslagkosten en verplaatsingen van het aantal containers

Data Science recente blogs

Bayesiaanse Statistiek — door Marc Jacobs — last modified 25-07-2022
Wiskundig raamwerk voor ouderwets leren
R vs Python — door Nick van de Venn — last modified 05-07-2022
Samenwerken is de sleutel
Machine Learning: De gereedschapskist van de data scientist — door marco — last modified 21-12-2021
Machine Learning algoritmes zijn de gereedschappen voor een data scientist

Data Science Nieuws & Evenementen

Data science opleidingen gaan weer van start! — door marco — last modified 08-02-2022
Vanaf 21 september start ons succesnummer weer! Twaalf weken data science in R, we hebben er weer zin in
Aedes data science workshop 2 van 3 — door marco — last modified 07-02-2022
Voor Aedes organiseert DIKW drie workshops data science
DIKW in top 50 beste data science bedrijven — door marco — last modified 22-10-2021
DIKW is één van snelst groeiende bedrijven volgens MKB Data Science rapport