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:

reproductie getal covid-19 op basis van simulatie ziekenhuisopnames

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...)

Covid-19 weersverwachting

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!

Meer weten over wat data voor uw organisatie kan betekenen?

Neem contact op met Hugo Koopmans

Telefoon: +31643106780

Email:hugo.koopmans@dikw.com

CONTACT

BLOG OVERZICHT

Blogs

Data gedrevenheid is proces van lange adem — door Marco van den Doel — last modified 16-09-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 03-09-2021: Ondersteun uw customer journey met data strategie
Wat is data engineering? — door Marco van den Doel — last modified 03-09-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-09-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-09-2021: Transformeren naar een data gedreven organisatie kost tijd
Data gedreven logistiek onderhoud voorkomt uitval — door Marco van den Doel — last modified 06-09-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 03-09-2021: Duurzame innovatieve logistieke oplossing op basis van data science
Welke sandwich mogen wij voor u bereiden? — door Marco van den Doel — last modified 03-09-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 03-09-2021: DIKW is partner van Sensing Clues
In het verleden behaalde resultaten... — door Marco van den Doel — last modified 02-09-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 02-09-2021: Containerschip blokkeert Suezkanaal
Van voor naar achteren en van links naar rechts in de logistieke keten — door Marco van den Doel — last modified 02-09-2021: Het verminderen van opslagkosten en verplaatsingen van het aantal containers
De zeven pilaren van DataOps — door Marco van den Doel — last modified 03-09-2021: DataOps wordt gedefinieerd door zeven hoofdkenmerken
Kijk verder dan je dashboard — door Marco van den Doel — last modified 02-09-2021: Met de DIKW Analytical Roadmap kijk je verder!
Data discovery tools — door Marco van den Doel — last modified 03-09-2021: 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 03-09-2021: 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 03-09-2021: Ter ondersteuning van het corona dashboard van de rijksoverheid
Koning TOTO : Sjaak vs Bayes — door Marco van den Doel — last modified 03-09-2021: Definitieve uitslag voetbal eredivisie op basis van een wiskundig model, Baysiaanse statistiek, en een kleine Monte Carlo simulatie
COVID-19 Weersverwachting — door Marco van den Doel — last modified 03-09-2021: Ter ondersteuning van het corona dashboard van de rijksoverheid
Dashboard coronavirus — door Marco van den Doel — last modified 03-09-2021: Eerste observaties van een datascientist
Van BICC naar DACoE deel 1 — door Marco van den Doel — last modified 03-09-2021: Van BICC naar Data & Analytics Center of Excellence: waarom je moet veranderen om relevant te blijven
SERVAL Open Ears AI machine listening — door Marco van den Doel — last modified 03-09-2021: Building artificial ears for (urban) jungle applications
Granuliet WOB documenten — door Marco van den Doel — last modified 03-09-2021: Textmining LDA Topic Models toegepast op 2 GB aan WOB documenten over granuliet
BERT en Transformer Learners — door Marco van den Doel — last modified 02-09-2021: Ontwikkelingen op het gebied van Natural Language Processing
Textmining vs NLP — door Marco van den Doel — last modified 03-09-2021: De verschillen en toepassingen van textmining en Natural Language Processing
Wat is nieuw in IBM Cognos Analytics 11 — door Marco van den Doel — last modified 03-09-2021: Cognos Analytics ontwikkelt zich snel en voortvarend als een betrouwbaar self-service platform voor data analyse
Predicting butter prices — door Marco van den Doel — last modified 03-09-2021: Case study assignment for Certified Data Science Proffesional course DIKW Academy
Forecast in R — door Marco van den Doel — last modified 03-09-2021: Tutorial Forecast in R
Tijdreeksanalyse in R — door Marco van den Doel — last modified 03-09-2021: Tijdreeksanalyse ARIMA in R, Handleiding modelselectie in R

Data Science recente blogs

Data gedreven logistiek onderhoud voorkomt uitval — door Marco van den Doel — last modified 06-09-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 03-09-2021: Duurzame innovatieve logistieke oplossing op basis van data science
Met data bijdragen aan een betere wereld — door Marco van den Doel — last modified 03-09-2021: DIKW is partner van Sensing Clues

Data Science Nieuws & Evenementen

DIKW in top 50 beste data science bedrijven — door Marco van den Doel — last modified 07-09-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 06-09-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 06-09-2021: De innovatie boost van AEDES is binnen gehaald door de werkgroep Big data.