COVID-19 Weersverwachting
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!
Meer weten over wat data voor uw organisatie kan betekenen?
Neem contact op met Hugo Koopmans
Telefoon: +31643106780
Email:hugo.koopmans@dikw.com
Blogs
-
Van pizzakoerier naar business intelligence consultant — door Nick van de Venn — last modified 15-12-2022
- Het verhaal hoe Brian een vaste waarde werd binnen DIKW
-
Data management — door Nick van de Venn — last modified 18-11-2022
- Met vertrouwen bouwen op data
-
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
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