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
-
Differential Privacy — door Nick van de Venn — last modified 18-09-2023
- Gevoelige gegevens verwerken zonder dat de gevoelige informatie kan uitlekken
-
Datagedreven werken Deel 2 — door Nick van de Venn — last modified 07-07-2023
- Welke data heeft eigenlijk waarde voor uw organisatie?
-
ChatGPT for Business Intelligence — door Nick van de Venn — last modified 18-09-2023
- Chatten met je datawarehouse, utopie of werkelijkheid?
-
Intelligence Factory — door Nick van de Venn — last modified 05-07-2023
- Agile design thinking met een ML-ops sausje
-
Data gedreven werken — door Nick van de Venn — last modified 26-05-2023
- Bij datagedreven werken staan mens én data centraal
-
Hoe gebruik je ChatGPT om je data pipeline te bouwen? — door Nick van de Venn — last modified 11-05-2023
- Marcel-Jan zocht het uit voor zijn grote hobby: astronomie.
-
Lead consultant en manager Business Analytics Patrick Meulstee — door Nick van de Venn — last modified 06-04-2023
- Over remote werken op Bonaire
-
Hoogwaardige data opleidingen bij DIKW Academy — door nick van de venn — last modified 28-02-2023
- Waarom wij geloven in de praktische opzet van onze cursussen
-
Lead data engineer Remy Lamberty — door Nick van de Venn — last modified 27-03-2023
- Over remote werken vanuit India
-
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
Data Science recente blogs
-
ChatGPT for Business Intelligence — door Nick van de Venn — last modified 18-09-2023
- Chatten met je datawarehouse, utopie of werkelijkheid?
-
Intelligence Factory — door Nick van de Venn — last modified 05-07-2023
- Agile design thinking met een ML-ops sausje
-
Bayesiaanse Statistiek — door Marc Jacobs — last modified 25-07-2022
- Wiskundig raamwerk voor ouderwets leren
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