Plots credible intervals and median for the observed data under the posterior predictive distribution, and for a specific observation type. The user can control the interval levels (i.e. 30%, 50% etc.) and the plotted group(s). This is a generic function.
plot_obs(object, ...)
# S3 method for class 'epimodel'
plot_obs(
object,
type,
groups = NULL,
dates = NULL,
date_breaks = "2 weeks",
date_format = "%Y-%m-%d",
cumulative = FALSE,
by_100k = FALSE,
bar = TRUE,
levels = c(30, 60, 90),
log = FALSE,
...
)
spaghetti_obs(
object,
type,
draws = min(500, posterior_sample_size(object)),
alpha = 1/sqrt(draws),
groups = NULL,
dates = NULL,
date_breaks = "2 weeks",
date_format = "%Y-%m-%d",
cumulative = FALSE,
by_100k = FALSE,
bar = TRUE,
log = FALSE,
smooth = 1,
...
)A fitted model object returned by epim.
See epimodel-objects.
Additional arguments for
posterior_predict.epimodel. Examples include
newdata, which allows
predictions or counterfactuals.
A string specifying the name of the observations to plot. This should match one
of the names of the response variables in the obs argument used int the call to epim.
Either NULL or a character vector specifying the groups
to plot for. Default is NULL, which plots all modeled groups.
A length 2 vector of Date objects. This defines the
start and end dates of the date-range to be plotted. Must be coercible to
Date if not NA. If an element of the vector is NA then
the default lower/upper limit is used. See examples.
A string giving the distance between date tick labels.
Default is "2 weeks". This is passed as the date_breaks argument to
scale_x_date. Please see here for details.
This function attempts to coerce the dates argument
to a vector of Date objects. date_format is passed as the format
argument to as.Date. Default is "%Y-%m-%d".
If TRUE then cumulative observations are plotted rather than daily. Defaults to FALSE.
If TRUE, all quantities are plotted per 100k of population. Only possible if the model used a population adjustment.
If TRUE, observations are plotted as a bar plot. Otherwise, a scatterplot is used. Defaults to TRUE.
A numeric vector defining the levels of the plotted credible intervals.
If TRUE, plot quantities on a log10-scale. This argument
must be logical, and defaults to FALSE.
The number of sample paths to plot.
Sets transparency of sample paths.
An integer specifying the window used to smooth the reproduction rates. The
default is 1, which corresponds to no smoothing.
A ggplot object which can be further modified.
# Fitting requires an external CmdStan install (cmdstanr::install_cmdstan()).
# \donttest{
data("EuropeCovid2")
data <- EuropeCovid2$data
data <- dplyr::filter(data, date > date[which(cumsum(deaths) > 10)[1] - 30])
data <- dplyr::filter(data, date < as.Date("2020-05-05"))
rt <- epirt(
formula = R(country, date) ~ 0 + (1 + public_events + schools_universities +
self_isolating_if_ill + social_distancing_encouraged + lockdown || country) +
public_events + schools_universities + self_isolating_if_ill +
social_distancing_encouraged + lockdown,
prior = shifted_gamma(shape=1/6, scale = 1, shift = log(1.05)/6),
prior_covariance = decov(shape = c(2, rep(0.5, 5)),scale=0.25),
link = scaled_logit(6.5)
)
inf <- epiinf(gen = EuropeCovid$si, seed_days = 6)
deaths <- epiobs(
formula = deaths ~ 1,
i2o = EuropeCovid2$inf2death,
prior_intercept = normal(0,0.2),
link = scaled_logit(0.02)
)
args <- list(rt=rt, inf=inf, obs=deaths, data=data, seed=12345)
args$group_subset <- c("Italy", "Austria", "Germany")
args$algorithm <- "fullrank"
args$iter <- 1e4
args$tol_rel_obj <- 1e-3
fm <- do.call(epim, args)
#> ------------------------------------------------------------
#> EXPERIMENTAL ALGORITHM:
#> This procedure has not been thoroughly tested and may be unstable
#> or buggy. The interface is subject to change.
#> ------------------------------------------------------------
#> Gradient evaluation took 0.000767 seconds
#> 1000 transitions using 10 leapfrog steps per transition would take 7.67 seconds.
#> Adjust your expectations accordingly!
#> Begin eta adaptation.
#> Iteration: 1 / 250 [ 0%] (Adaptation)
#> Iteration: 50 / 250 [ 20%] (Adaptation)
#> Iteration: 100 / 250 [ 40%] (Adaptation)
#> Iteration: 150 / 250 [ 60%] (Adaptation)
#> Iteration: 200 / 250 [ 80%] (Adaptation)
#> Iteration: 250 / 250 [100%] (Adaptation)
#> Success! Found best value [eta = 0.1].
#> Begin stochastic gradient ascent.
#> iter ELBO delta_ELBO_mean delta_ELBO_med notes
#> 100 -12564.656 1.000 1.000
#> 200 -7616.279 0.825 1.000
#> 300 -4418.754 0.791 0.724
#> 400 -3482.389 0.661 0.724
#> 500 -3286.507 0.540 0.650
#> 600 -2748.405 0.483 0.650
#> 700 -2429.817 0.433 0.269
#> 800 -2466.154 0.380 0.269
#> 900 -1858.797 0.374 0.269
#> 1000 -1783.876 0.341 0.269
#> 1100 -1606.415 0.252 0.196
#> 1200 -1608.129 0.187 0.131
#> 1300 -1412.967 0.129 0.131
#> 1400 -1418.648 0.102 0.110
#> 1500 -1239.419 0.111 0.131
#> 1600 -1187.959 0.096 0.110
#> 1700 -1188.581 0.083 0.043
#> 1800 -1136.147 0.086 0.046
#> 1900 -1161.259 0.055 0.043
#> 2000 -1086.384 0.058 0.046
#> 2100 -1102.434 0.048 0.043
#> 2200 -1028.861 0.055 0.046
#> 2300 -1029.647 0.042 0.043
#> 2400 -1015.375 0.043 0.043
#> 2500 -991.992 0.030 0.024
#> 2600 -987.169 0.027 0.022
#> 2700 -985.940 0.027 0.022
#> 2800 -998.150 0.023 0.015
#> 2900 -959.803 0.025 0.015
#> 3000 -984.480 0.021 0.015
#> 3100 -972.671 0.021 0.014
#> 3200 -958.460 0.015 0.014
#> 3300 -961.169 0.015 0.014
#> 3400 -947.345 0.015 0.015
#> 3500 -956.334 0.014 0.012
#> 3600 -956.851 0.013 0.012
#> 3700 -946.955 0.014 0.012
#> 3800 -944.959 0.013 0.012
#> 3900 -950.063 0.010 0.010
#> 4000 -951.692 0.007 0.009
#> 4100 -945.356 0.007 0.007
#> 4200 -935.736 0.006 0.007
#> 4300 -946.324 0.007 0.009
#> 4400 -936.581 0.007 0.009
#> 4500 -929.176 0.007 0.008
#> 4600 -935.461 0.007 0.008
#> 4700 -939.771 0.007 0.007
#> 4800 -930.044 0.008 0.008
#> 4900 -929.907 0.007 0.008
#> 5000 -925.010 0.007 0.008
#> 5100 -930.120 0.007 0.008
#> 5200 -934.518 0.007 0.007
#> 5300 -928.864 0.006 0.006
#> 5400 -926.138 0.005 0.005
#> 5500 -924.682 0.005 0.005
#> 5600 -923.552 0.004 0.005
#> 5700 -925.678 0.004 0.005
#> 5800 -930.174 0.003 0.005
#> 5900 -920.703 0.004 0.005
#> 6000 -923.863 0.004 0.005
#> 6100 -915.439 0.005 0.005
#> 6200 -917.973 0.004 0.003
#> 6300 -925.978 0.005 0.003
#> 6400 -917.889 0.005 0.005
#> 6500 -921.228 0.006 0.005
#> 6600 -916.181 0.006 0.006
#> 6700 -915.070 0.006 0.006
#> 6800 -919.995 0.006 0.006
#> 6900 -919.765 0.005 0.005
#> 7000 -914.784 0.005 0.005
#> 7100 -922.327 0.005 0.005
#> 7200 -917.651 0.005 0.005
#> 7300 -911.132 0.005 0.005
#> 7400 -909.527 0.004 0.005
#> 7500 -912.507 0.004 0.005
#> 7600 -912.139 0.004 0.005
#> 7700 -912.779 0.004 0.005
#> 7800 -912.633 0.003 0.003
#> 7900 -908.169 0.004 0.005
#> 8000 -917.951 0.004 0.005
#> 8100 -910.985 0.004 0.005
#> 8200 -913.086 0.004 0.003
#> 8300 -907.917 0.004 0.003
#> 8400 -909.226 0.004 0.003
#> 8500 -912.256 0.004 0.003
#> 8600 -909.347 0.004 0.003
#> 8700 -909.429 0.004 0.003
#> 8800 -914.688 0.005 0.005
#> 8900 -907.193 0.005 0.006
#> 9000 -912.019 0.004 0.005
#> 9100 -909.071 0.004 0.003
#> 9200 -910.826 0.004 0.003
#> 9300 -906.600 0.004 0.003
#> 9400 -912.381 0.004 0.005
#> 9500 -907.631 0.004 0.005
#> 9600 -906.583 0.004 0.005
#> 9700 -905.301 0.004 0.005
#> 9800 -905.782 0.004 0.005
#> 9900 -903.136 0.003 0.003
#> 10000 -904.283 0.003 0.003
#> Informational Message: The maximum number of iterations is reached! The algorithm may not have converged.
#> This variational approximation is not guaranteed to be meaningful.
#> Drawing a sample of size 1000 from the approximate posterior...
#> COMPLETED.
#> Finished in 3.0 seconds.
# different ways of using plot_rt
p <- plot_rt(fm) # default, plots all groups and dates
#> Running standalone generated quantities after 1 MCMC chain...
#>
#> Chain 1 finished in 0.0 seconds.
p <- plot_rt(fm, dates=c("2020-03-21", NA)) # plot 21 March 2020 onwards
#> Running standalone generated quantities after 1 MCMC chain...
#>
#> Chain 1 finished in 0.0 seconds.
p <- plot_rt(fm, dates=c(NA, "2020-03-20")) # plot up to 20 March 2020
#> Running standalone generated quantities after 1 MCMC chain...
#>
#> Chain 1 finished in 0.0 seconds.
p <- plot_rt(fm, dates=c("2020-03-20", "2020-04-20"))
#> Running standalone generated quantities after 1 MCMC chain...
#>
#> Chain 1 finished in 0.0 seconds.
p <- plot_rt(fm,
dates=c("2020-03-20", "2020-04-20"),
date_breaks="1 day") # ticks every day
#> Running standalone generated quantities after 1 MCMC chain...
#>
#> Chain 1 finished in 0.0 seconds.
p <- plot_rt(fm,
dates=c("2020-20-03", "2020-20-04"),
date_format="%Y-%d-%m") # (different date format)
#> Running standalone generated quantities after 1 MCMC chain...
#>
#> Chain 1 finished in 0.0 seconds.
# other plotting functions
p <- plot_obs(fm, type = "deaths")
#> Running standalone generated quantities after 1 MCMC chain...
#>
#> Chain 1 finished in 0.0 seconds.
p <- plot_infections(fm)
#> Running standalone generated quantities after 1 MCMC chain...
#>
#> Chain 1 finished in 0.0 seconds.
p <- plot_infectious(fm)
#> Running standalone generated quantities after 1 MCMC chain...
#>
#> Chain 1 finished in 0.0 seconds.
# }