key: cord-0944073-3q9wcbgg
authors: Congdon, Peter
title: Mid-Epidemic Forecasts of COVID-19 Cases and Deaths: A Bivariate Model Applied to the UK
date: 2021-02-12
journal: Interdiscip Perspect Infect Dis
DOI: 10.1155/2021/8847116
sha: 9619040c71cfbf0e30a2ab142e1ce5bd28d8bde6
doc_id: 944073
cord_uid: 3q9wcbgg

BACKGROUND: The evolution of the COVID-19 epidemic has been accompanied by efforts to provide comparable international data on new cases and deaths. There is also accumulating evidence on the epidemiological parameters underlying COVID-19. Hence, there is potential for epidemic models providing mid-term forecasts of the epidemic trajectory using such information. The effectiveness of lockdown or lockdown relaxation can also be assessed by modelling later epidemic stages, possibly using a multiphase epidemic model. METHODS: Commonly applied methods to analyse epidemic trajectories or make forecasts include phenomenological growth models (e.g., the Richards family of densities) and variants of the susceptible-infected-recovered (SIR) compartment model. Here, we focus on a practical forecasting approach, applied to interim UK COVID data, using a bivariate Reynolds model (for cases and deaths), with implementation based on Bayesian inference. We show the utility of informative priors in developing and estimating the model and compare error densities (Poisson-gamma, Poisson-lognormal, and Poisson-log-Student) for overdispersed data on new cases and deaths. We use cross validation to assess medium-term forecasts. We also consider the longer-term postlockdown epidemic profile to assess epidemic containment, using a two-phase model. RESULTS: Fit to interim mid-epidemic data show better fit to training data and better cross-validation performance for a Poisson-log-Student model. Estimation of longer-term epidemic data after lockdown relaxation, characterised by protracted slow downturn and then upturn in cases, casts doubt on effective containment. CONCLUSIONS: Many applications of phenomenological models have been to complete epidemics. However, evaluation of such models based simply on their fit to observed data may give only a partial picture, and cross validation against actual trends is also valuable. Similarly, it may be preferable to model incidence rather than cumulative data, although this raises questions about suitable error densities for modelling often erratic fluctuations. Hence, there may be utility in evaluating alternative error assumptions.

Epidemic forecasts have been an essential element in policy decisions regarding the COVID-19 epidemic, such as lockdown imposition and relaxation. Forecasting has been assisted by well-organized efforts to provide international data on new cases and deaths. ese include the daily updated comparative data provided by the European Centre for Disease Prevention and Control (ECDC) (https://www. ecdc.europa.eu/en/publications-data), and monitoring profiles provided by John Hopkins University [1] . ere is also a growing literature providing evidence on the parameters of the COVID-19 infection (for example, case fatality ratios and serial intervals). Hence, the potential occurs for epidemic models that are applicable to routinely collected data, which make use of accumulated evidence, and can provide forecasts for epidemics observed at mid-stage. Policy decisions in many countries (imposition of lockdowns, and later, relaxation) have been made based on trends in observed numbers of cases and deaths, while admitting these may be subject to measurement error, for example, identified cases may understate total numbers infected; there are fluctuations in daily new cases due to variation in daily testing, and there may be COVID-19 diagnostic errors.

Here, we focus on a practical forecasting approach using routinely available data on new cases and deaths (from ECDC) to estimate parameters in a bivariate version of the Reynolds phenomenological model. Implementation is based on Bayesian inference principles, incorporating accumulated evidence on relevant parameters via informative priors. e operation and utility of this approach is demonstrated using data on new cases and deaths in the United Kingdom (UK), with a focus on predictive accuracy for 20day-ahead forecasts of cases and deaths, based on mid-epidemic data. Other policy relevant parameters such as the effective reproduction ratio and case fatality ratio are also estimated in analysis of longer-term epidemic data. e use of a bivariate approach provides originality compared to existing research, which, in the case of phenomenological models, is limited to analysing incidence only. Some studies have mentioned how the mortality curve parallels the epidemic curve [2] , and we formalize these linkages under a bivariate approach. e benefits of a bivariate approach include the ability to monitor and forecast severity measures such as the case fatality ratio. We also consider issues in modelling daily incidence and new deaths, as opposed to cumulative incidence and mortality. ere are methodological issues in analysing cumulative outcomes, discussed below, but also questions (not so far considered in the literature) on how best to represent the overdispersion present in uncumulated outcomes. e analysis below provides new evidence on the relative fit and forecasting performance of different ways of representing Poisson overdispersion and shows that the usually adopted negative binomial performs less well than other options. e implications of the research are that effective medium-term forecasts of COVID-19 incidence and mortality can be provided by the proposed methodology. Such forecasts are useful in planning healthcare provision and assessing closeness to full capacity in hospital bed occupancy [3, 4] . At the time of writing, daily COVID-19 hospitalisation data were not available for the UK. However, extension of the bivariate approach to include incidence, mortality, and hospitalisations provides additional scope for forecasting of indicators relevant to severity assessment [5, 6] . Another possible outcome in a multivariate model is recovered cases (cf. [4] ), with a focus then on the ratio of predicted recovered to predicted confirmed cases as an indicator of care need and effectiveness of interventions. e methodology also presents a way to monitor longer-term infection numbers leading to early detection of incipient upturns in infection numbers, via continuous monitoring and forecasting of the effective reproduction ratio.

In the following sections, we review relevant literature and research gaps (Section 2) and set out the methodology (Section 3). We then consider aspects of the UK case study application (Section 4), present the results (Section 5), and discuss the implications of the study's findings and methodology in the context of broader research (Section 6). Section 7 provides concluding remarks.

Commonly applied approaches to quantitative modelling of aggregate epidemic data differ in the data inputs they require, their assumptions, their estimability, and their scope for practical application to making forecasts. Commonly applied methods include phenomenological growth models [7] , such as the Richards family of densities [8] , and variants of the susceptible-infected-recovered (SIR) compartment model [9] . Phenomenological models are parameterised in terms of epidemic trajectories and provide estimates of crucial epidemiological parameters [10, 11] , while avoiding the complexity of more formal mechanistic models of disease transmission, which can be difficult to estimate and provide forecasts, and may not be realistic approximations to real epidemic dynamics [2, 9, 12] .

As mentioned by Chowell et al. [13] , phenomenological models are "particularly suitable when significant uncertainty clouds the epidemiology of an infectious disease." By contrast, as noted in [14] , compartmental transmission models may be based on untested assumptions such as random mixing between all individuals in a given population, may be sensitive to starting assumptions, and may provide estimates that differ considerably between models. Such models often rely on preset parameters, which may mean prediction uncertainty is understated. ey may also be complex to specify when an epidemic has more than one phase, whereas multiphase phenomenological models [15] are available.

Regarding forecasts, the study by Zhao et al. [10] exemplifies application of phenomenological models to forecasts of the Zika epidemic in 2015. Autoregressive modelling of new cases, with potential for short-term forecasting, is illustrated (for foot and mouth disease) by the first-order autoregressive model of Lawson et al. [16] , while (for multiple spatial units) the model of Bracher and Held [17] specifies a first-order autoregression based on the mean incidence in adjacent areas.

Regarding the COVID-19 epidemic in particular, studies with differing methodologies have been made, some of which forecast different aspects of the COVID epidemic or related health care need. In its impact on policy making in the UK and US, perhaps most influential has been the Imperial College model [18] . is is based on microsimulation with transmission through contacts between susceptible and infectious individuals in various settings, or randomly in the community, depending on spatial distance between contacts. A number of epidemic parameters (e.g., incubation periods and basic reproduction numbers) are preset. Forecasts are provided for deaths and hospital beds. Also, providing forecasts across countries is the model of the Institute for Health Metrics and Evaluation [14] . is has no underlying representation of epidemic dynamics, but is based on fitting a hierarchical parametric model for observed cumulative death rates in different countries, and then projecting these forward.

Various types of time series forecast of the COVID-19 epidemic have also been made, using ARIMA models [4, 19, 20] , exponential smoothing [21] , or autoregression in cases; for example, the study by Johndron et al [22] postulates daily deaths as a lagged function of earlier new cases. Applications of phenomenological models to COVID-19 incidence forecasts include the studies in [23, 24] .

One may identify some research gaps in the existing literature. us, existing applications (e.g., [23] in the case of COVID-19) have most commonly been to incidence data and have not considered interplay between outcomes in terms of a bivariate model. However, examples of potentially interlinked bivariate outcomes from COVID-19, and other epidemics, include incidence and mortality, and incidence and hospitalisations. e study in [25] converts COVID-19 incidence, mortality, and case fatality into beta variables and applies univariate beta regression models to each outcome. e study in [4] applies separate ARIMA models to incidence and recovered cases. However, separate univariate regressions or time series models do not reflect potential interlinkages between the processes that facilitate the setting of model assumptions or, in the case of Bayesian analysis, facilitate the setting of priors on model parameters. Regarding forecasting, wider experience of time series modelling (with nonepidemic applications) shows the benefit of borrowing strength over outcomes [26] .

Furthermore, many existing applications, such as studies in [11, 25] in the case of COVID-19, study in [27] in the case of H1N1 and Ebola, and study in [10] in the case of Zika, have been to cumulative incidence. However, the drawbacks of studying cumulative incidence have been pointed out [28, 29] . Cumulative incidence data have serially correlated measurement error, leading to understatement of parameter uncertainty. As stated in [29] , "independence of sequential measurement errors, . . .is clearly violated when observations are accumulated through time." However, estimation using new cases or deaths (uncumulated) puts a much greater focus on how to deal with stochastic variation in the data. For daily data, fluctuations in new events may be considerable (including daily "spikes"), whereas cumulative cases and deaths are usually smooth functions.

ere are choices in how to model the overdispersion in new events, based on Poisson mixtures [30] , but in the epidemic literature, use of the negative binomial is standard, and no evaluations of its relative performance are available in the literature so far.

Basic phenomenological models for epidemic trajectories include the logistic, Gompertz, and Rosenzweig, which have been the basis for a range of generalisations [7, 8] . Application of the logistic model to COVID-19 is exemplified by the studies of Batista [31] and Shen [32] . For time t, the logistic model for new cases C ′ (t) and cumulative cases C(t) is

where K is the maximum number of cases (final epidemic size), r > 0 measures the intensity of exponential growth in cases in the early epidemic phase, and τ L is the inflection point where new cases are highest. e Richards model [33] modifies the logistic incidence function to

with solution

e parameter a > 0 modifies the incidence decline phase of the logistic, that is, measures the extent of deviation from the standard logistic curve. e turning point τ, when incidence peaks, is obtained when C(t) equals K(1 + a) − 1/a [34] . e peak incidence is important for the healthcare planning, for example, aligning the forecast peak with hospital bed capacity [35] .

Other commonly used models are the Gompertz model [36] with

while the Rosenzweig model [27] has

e incidence function represented by C ′ (t) can be used to define mean incidence in statistical likelihoods for new cases data. us, time series of incidence counts can often be satisfactorily modelled as a Poisson, with means defined by C ′ (t) functions [37, 38] . Similarly, the cumulative cases function C(t) can be used to define mean epidemic size in models for cumulative case counts [27, 39] .

While for smaller epidemics, a Poisson density for mean incidence may be applied, for larger epidemics such as COVID-19, a negative binomial density is often preferred, both because of large incidence counts and to represent often erratic incidence fluctuations that lead to overdispersion relative to the Poisson [29] . However, the literature does not contain any evaluations of the negative binomial to represent overdispersion. ere are a number of other overdispersed versions of the Poisson that can be achieved by mixing the Poisson with a suitable density (e.g., a lognormal density) [30, 40] , and this may be beneficial in detecting unusual observations.

Overdispersion. Consider the Richards model parameters. Let c t and C t denote incidence and cumulative incidence counts at times t � 1, . . . , T (days in the case of COVID data from ECDC). We condition on the first case or cases (i.e., the first observation) and take incidence at time t as a function of cumulative cases at t − 1, so that for a Poisson model for cases (with subscript c for parameters) we have

In practice, many epidemic datasets are overdispersed relative to the Poisson, and epidemic studies generally adopt a negative binomial model instead. We can specify an overdispersed model (including the negative binomial) by introducing multiplicative random effects ( [40] , Equation ( 2)), such that the Poisson means for incidence are specified by

where ϵ ct are positive random effects. For the Poissongamma model (equivalent to a negative binomial), ϵ ct are gamma distributed with mean 1, namely,

where 1/λ c is the overdispersion parameter mentioned by [29] . Note that the assumed parameterisation of the gamma density with random variable x is

Other options in (7) are to take u ct � log(ϵ ct ) as normally distributed [40] :

or Student t distributed [41] ,

where ] c is a degrees-of-freedom parameter. ese two options define the Poisson-lognormal (PLN) and Poissonlog-Student (PLS) options, respectively [30] . e PLN and PLS representations may provide a more robust alternative to the Poisson-gamma [40] [41] [42] [43] [44] , as their tails are heavier than for the gamma distribution, and have been found to be better at accommodating outliers (such as daily "spikes" in an epidemic application).

In the analysis below, we apply a bivariate estimation with both new cases and deaths modelled using the Richards specification. us, denote d t and D t as new and cumulative deaths at time t. e joint likelihood for an overdispersed Poisson model for both outcomes then specifies

For a Bayesian application, we need to specify prior densities, or priors for short, for the parameters. For the epidemic size parameter K c , a diffuse prior confined to positive values, such as a diffuse gamma density, for example, Gamma(1, ε) or Gamma(ε, ε), with ε small, was found to lead to convergence problems. As noted in [31] , ". . .in the early stage, the logistic curve follows an exponential growth curve, so the estimation of K is practically impossible." is difficulty persists when an epidemic is past its peak but early in a downturn.

However, Batista [31] mentions a relationship (for the logistic model) between successive cumulative case counts that may assist in providing a prior for K c . Specifically, for three points spaced m time units apart, one may obtain the relationship for a point estimator of K c , namely,

One may use this point estimator to define a prior mean for K c in the Richards model (which is a generalisation of the logistic). Specifically, one may take a lognormal density prior for K c , with log(K e c ) as mean, and a suitable variance, such that the prior is still relatively diffuse. For example, suppose K e c is 250,000, and the variance in the lognormal is set at 1. en, the 97.5 percentile for the lognormal prior is 1.77 million.

In the bivariate specification (for new cases and new deaths jointly), we seek to share prior information between outcomes. One option for the prior on K d (the final death total) is as a function of K c , namely,

where Φ is a form of case fatality ratio (CFR). An informative prior for Φ could be based on the COVID experience in similar countries, or on experience of epidemics of similar diseases. Considering the first option, and an appropriate prior for analysing UK data, an informative prior for Φ could be provided by the case fatality ratio across the European Union (the UK being no longer an EU member). International information on case fatality is provided at https://ourworldindata.org/mortality-risk-covid#thecurrent-case-fatality-rate-of-covid-19. Alternatively, one may link K c and K d using both the point estimator K e c and a case fatality ratio; namely, K e d � ΦK e c . en, a lognormal density prior can be taken for K d , with log(K e d ) � log(K e c ) + log(Φ) as the mean, and a suitable variance such that the prior is still relatively diffuse.

Another possible prior to link ultimate cases and deaths would involve a time series in time-specific case fatality ratios, such as an autoregression ϕ t ∼ N(ρϕ t−1 , σ 2 ϕ ), with ϕ t estimated from cumulative data on deaths and cases. Some analyses of epidemics show that the CFR early in an epidemic may underestimate later values [45] , in which case the prior on Φ may be constrained to exceed the final ϕ t based on observed data. With regard to COVID-19, this pattern may not necessarily apply, with the US (for example) showing a decline in CFRs at later epidemic stages. ere is also evidence that the mortality to infection ratio (a more precise measure than the CFR) has fallen [46] . To allow for such a scenario, the prior for Φ could be centred at the last observed ϕ t , rather than constrained to exceed it.

Joint priors on other parameters could be considered, for example, bivariate normal priors on the logs of r c and r d , or on the logs of a c and a d . In the empirical analysis below, we focus on priors linking the final epidemic and death total parameters, K c and K d , as these are an important influence on forecasts.

Many applications of phenomenological models are to historic data on epidemics, where the epidemic has run to its full extent. Here, we consider applications to incomplete epidemics (e.g., epidemics observed to their mid-point or early in the downturn), and to forecasts using such data. Forecasts at an intermediate point within the observation span are of interest in themselves for policy purposes. However, they can also be used in comparative model evaluation by using cross validation, with only some data used for estimating the model, and some held out for validation.

us, suppose the training sample is formed by ob- where I(A) � 1 if condition A is true, and I(A) � 0 otherwise. us, at each iteration, we compare average new cases (modelled) during the validation period with actual average new cases.

Probabilities of overprediction for new cases, ω c , are estimated as S s�1 (O (c) s /S). Probabilities of underprediction can be obtained as 1 − ω c . A satisfactory prediction would have 0.05 < ω c < 0.95, with ω c over 0.95 indicating a high probability of overprediction, while ω c under 0.05 indicates a high probability of underprediction. Underprediction means underforecasting of future cases and may lead to incorrect inferences regarding epidemic control, as it implies a lessening in incidence earlier than actually occurred.

Reproduction Ratios, and Case Fatality Ratios. Strategic decisions regarding containment of the COVID-19 epidemic, in the UK and other countries, have depended on trends in new infections and deaths, but also on the effective reproduction rate. us, in the UK, the choice on whether or not to relax the initial COVID lockdown restrictions was based on five criteria, with two being numeric: first, "a sustained and consistent fall in daily death rates" and second that the "rate of infection is decreasing to manageable levels," meaning that the effective reproduction ratio is demonstrably below 1. e reproduction rate may also become especially relevant at later epidemic stages (postlockdown), after a downturn from the initial peak and after lockdown measures have been relaxed. Here, the concern is to prevent a resurgence of infection, indicated by an upturn in R t .

In the case of a protracted downturn, but with new cases still occurring, the concern is especially that there may be a substantial resurgence in cases, and possibly also deaths, at some point. is scenario is colloquially known as a "second wave," and in most European countries, there have been pronounced second waves in the COVID-19 epidemic during 2020, albeit at different times. Such a resurgence indicates use of a multiphase model [15] , with a second phenomenological model applied to data after a latent switch-point between epidemic regimes. Note that at the time of writing (August 2020), a fully developed second wave had not yet happened in the UK, although the signs were of an upturn in R t , as the analysis below confirms.

In planning for hospital care, longer-term trends in disease severity may be relevant. For example, an upturn in cases may reflect more cases among younger people at lower mortality risk. Hence, trends in, and forecasts of, the case fatality rate are an important aspect of strategic management. In a bivariate model of cases and deaths, as here, we can trace the modelled CFR through time, where the modelled CFR at day t is given by the ratio of predicted total deaths to predicted total cases D new,t /C new,t . In a natural extension of the bivariate model to include hospitalisations, trends in, and forecasts of, hospitalisation rates (ratios of hospitalised cases to total cases) can also be estimated.

Assume there is some evidence from new cases data of an upturn in cases, even before a second wave epidemic is fully established. Consider a model for new cases only. en, a two-phase model to reflect the upturn would involve two phenomenological models, before and after a latent switchpoint, with each model having distinct parameters. Denote the single switch-point as κ c , such that for a two-phase Richards model,

where I(A) � 1 when condition A is true, and 0 otherwise. Parameters, such as (r c1 , r c2 ), are differentiated by outcome and by wave. e parameter κ c can be assigned a uniform prior (on a positive interval) or a positive valued prior, such as an exponential density. If there is a second wave upturn in deaths also, then a bivariate model can be used, with switchpoints κ d assumed later for deaths than cases, due to possible delays in mortality upturns following incidence upturns. A three-phase model for cases would have two switch-points, with mean

An estimator of the effective reproduction rate R t at time t is based on predictions from this phenomenological model and from an estimate of the serial interval density. e serial interval is the time between symptom onset in an infected subject and symptom onset in the infectee. e serial interval density can be discretised in the form of weights ρ j , applied to serial interval lengths (in days) up to a maximum J. ese can be used to estimate effective reproduction ratios R t within a phenomenological model to analyse new cases; see the papers [37, 48, 49] . us, R t can be estimated as

where c new,t are predicted new case data from the phenomenological model. By virtue of the MCMC sampling strategy used below, we can readily obtain 95% credible intervals for R t , and the probability that R t < 1, which is important for assessing epidemic containment strategy.

We consider the application of the above methods to UK data on new cases and deaths from 1st February 2020 (when the first two cases of COVID-19 in the UK were reported according to ECDC). Observations are assigned dates as in the ECDC data, with times t in days from 1st February 2020. Figures 1(a) and 1(b) show daily trends in these outcomes up to 8th August 2020, with erratic fluctuations apparent in both outcomes. However, there is a broad downward trend in both outcomes from days 70 to 80, although with a more protracted decline as opposed to the steep initial upturn. Figures 2(a) and 2(b) show the relatively smooth evolution of cumulative cases and deaths. 

where ] c is the degree of freedom parameter and u ct � log(ϵ ct ), we take

e indicators δ ct have average 1 but are significantly lower than 1 for more outlying observations, such as daily spikes in cases.

Cross-validation estimations are made at points M < T. us, we consider twenty-day-ahead forecasts at three different stages of the UK COVID-19 epidemic. For the first cross validation, estimations are based on training data up to day M � 80 with F � 20 (i.e., the cross-validation period consists of days 81 to 100). Cross-validation estimations with F � 20 are also made for M � 100, and for M � 120, with forecast accuracy based on comparing forecasts with holdout data for days M + 1, . . . , M + F.

For prior densities on the unknowns in the medium-term forecasts, exponential priors with mean 1 are assumed on r c , r d , a c , and a d . For the precisions 1/σ 2 uc and 1/σ 2 ud in the PLN and PLS options, and for the parameters λ c and λ d in the Poisson-gamma model, a gamma prior Gamma (1, 0.001) is assumed. For the PLS option, we take ] c � ] d � 4 as a preset option. e degrees of freedom can be difficult to estimate for relatively small datasets, and the option of the preset value ] � 4 is a robust option [52, 53] . For the maximum cases (epidemic size) parameter K c , we assume as lognormal prior, centred at log(K e c ) and with variance 1, as discussed in Section 3.2. For the maximum deaths parameter K d , we assume K d � ΦK c , where the prior for Φ is beta distributed with mean defined by the EU-wide case fatality ϕ EU . e beta has total prior count C set at 5, C � 5. us, the prior on Φ is Beta(Cϕ EU , C(1 − ϕ EU )). For example, on 20-04-2020 (at day t � 80), the EU-wide case fatality rate was 0.101, and the prior is set at Φ ∼ Beta (0.505, 4.495). is illustrates the use of relevant external information, rather than diffuse priors.

Fit to the observed training data is assessed using the Watanabe-Akaike information criterion (WAIC) [54] .

is is a fit measure which takes account of fit but also penalises model complexity. Crossvalidation fit (fit out of the observed sample) is assessed by the probabilities of forecasting overprediction, ω c and ω d , and by predictive coverage: whether the credible intervals for predicted cases and deaths (averages over each validation period) include the observed averages. Estimation uses the BUGS program [55] , with posterior estimates based on the second halves of two chain runs of 100,000 iterations, and convergence assessed using Brooks-Gelman-Rubin criteria [56] .

To provide a longerterm perspective on epidemic containment, we apply the best performing model from analysis of the first wave to the full set of observations as at 8th August 2020 (T �190), when the first epidemic peak had passed. So this application is to a situation where the epidemic first wave has passed, but there are still nonnegligible numbers of new cases, and a potential for possible upturns and further waves.

New cases and deaths (as 7-day moving averages) had reached maxima of 4850 and 950, respectively, when the UK epidemic peaked in April 2020. As a result of lockdown measures imposed in late March, daily new cases averaged just over 500 daily by July. However, lockdown relaxations from July were accompanied by the risk of resurgence. In that regard, an upturn was apparent with new cases in late July and early August averaging over 800 daily (see Figure 3 ), although deaths continued to fall, averaging around 50 per day in early August.

To model the full time series, and since there was no upturn in deaths apparent in early August 2020, we focus on new cases only. To reflect the evidence of an upturn in cases, a two-phase model is applied. For this model, the shift parameter κ c is assigned an exponential prior density with mean 150. Priors on the parameters of the second-phase Richards model are as before for the exponential ascent and logistic modifier parameters. For K c2 (the forecast total cases under wave 2), we assume K c2 � η c K c1 with η c assigned an exponential prior with mean 1. Of policy interest here is epidemic containment, as summarised by the effective reproduction ratio: specifically the question of policy relevance Interdisciplinary Perspectives on Infectious Diseases is whether this ratio consistently is below 1, and whether its 95% credible interval also is entirely below 1. If the ratio is not below 1, this suggests a significant upturn. To provide estimates of the effective reproduction ratio, we use accumulated evidence on COVID serial intervals from five studies [57] [58] [59] [60] [61] . A gamma density on the serial interval (SI) is assumed, and information on mean and standard deviation of the SI, or on quantiles of SI, is converted to gamma density parameters; for use of SI quantiles in this regard, see [62] . Large samples (of a million) from each the five densities are taken, and parameters of the pooled gamma density are estimated from the pooled sample of five million; the pooled gamma density has shape 1.38 and rate 0.36. is density is then converted to a discretised form (with 16 bins) to provide an informative prior on the SI to the model of White and Pagano [63] , which updates the prior serial interval density using new case data for the UK. e study in [64] recommends that the initial, approximately exponential, epidemic phase be used in estimation, and we use UK new case data up to time 24-04-2020, when UK new cases peaked (see 55, page 3). e updated mean serial interval is estimated as 3.5 days with standard deviation 3.1. e discretised serial interval is estimated as in [63] , with J in Section 3.5 set at J � 16. Table 1 compares parameter estimates from the observed (training) data under the three error assumptions for the three cross-validation analyses at M � 80, M � 100, and M � 120 (i.e., nine scenarios). It should be noted that predictions of cases and deaths should be based not on the posterior mean or median parameter values in Table 1 , but on sampled posterior predictive replicate data at each iteration. ese are based on sampling new data from the Poisson means (13) and (14) , and from the Richards model parameter profile at each particular iteration. e predicted values of cases and deaths are very close to actual values: for example, for the PLS model at M � 120, the average absolute deviation (over t � 2, . . . , 120) between actual new cases and predicted new cases is under 1. Table 2 compares the WAIC fits to the observed data under the nine options, while criteria regarding 20-day-ahead forecasts are shown in Table 3 . Table 1 shows broad consistency between the three distributional options in terms of estimated final epidemic size K c and eventual death total K d . e posterior density of these parameters may be skewed, with posterior mean exceeding median. For the PLN and PLS options, the estimate of K c increases as M does. is reflects the protracted nature of the UK downturn in cases after the peak in the first wave. Posterior mean estimates of the turning points τ c and τ d vary slightly and tend to be higher for M � 100 and M � 120, but for both outcomes and all M values are between 72 and 86. Turning point estimates for new cases τ c are also mostly higher under the PLS option. Table 2 shows that the PLS option has better fit, with lower WAIC values. Hence, its estimates of epidemic size and turning points are preferred and provide a better description of the slow decline in cases from their peak. Table 3 shows generally better cross-validation performance for the PLS option. As discussed above, satisfactory prediction would have 0.05 < ω c < 0.95, and 0.05 < ω d < 0.95, with ω c or ω d over 0.95 indicating overprediction, and with ω c or ω d under 0.05 indicating underprediction. Both PG and PLN options show underprediction for higher values of M, and in policy terms, may be misleading in suggesting a faster decline in cases and deaths than actually occurred. By contrast for M � 100 and M � 120, the 95% interval for predicted average cases under the PLS model comfortably includes the actual average new cases.

Estimates of policy relevant parameters also come from the longer-term scenario when the log-Student Richards model is applied to UK COVID new case data up to day T �190 (8th August, 2020). e shift point in the two-phase model is estimated as 179.1 with 95% interval from 179.0 to 179.3. Figure 4 plots out the posterior mean estimated effective reproductive ratios (5-day moving averages), distinguishing those estimates significantly above 1. e estimates hover around 1 throughout June and the first half of July but, in late July/early August, tend to exceed 1. e impression from this is that success in fully effective containment is then in doubt. e estimates of the reproduction ratio, and their path through June and July, are similar to those for the UK available at https://epiforecasts.io/covid/posts/national/ united-kingdom/#national-summary and based on the methods in [65] .

To illustrate the potential for longer-term forecasts of severity indicators, we also use observed data up to day 150 to make forecasts 40 days ahead through to day 190. is analysis uses the PLS model option for Poisson mixing. Figure 5 shows the resulting forecast of the case fatality ratio, with 95% intervals. e interval for t � 190 includes the observed value of 0.151.

e existing epidemic modelling literature has recognized the need for overdispersed distributions to deal with erratic incidence counts [29, 66, 67] . us, the study in [67] shows that use of a negative binomial distribution is more appropriate than the Poisson for describing emerging infections with overdispersed case distributions due to superspreading events. However, so far as the authors are aware, there has been no evaluation of the negative binomial as compared to other methods of representing overdispersion in epidemic counts. Hence, one contribution of this paper rests on a comparison of the negative binomial (Poisson-gamma mixture) against alternative Poisson mixture models [30] . For example, the Poisson log-normal distribution is a longer tailed alternative to the negative binomial distribution and may better fit overdispersed count data [68, 69] . e analysis here suggests such alternatives may usefully be considered as alternatives to the negative binomial and may both improve fit to actual observations and provide more accurate forecast performance.

Tackling overdispersion is one issue present in modelling count data associated with epidemics. Another is the outcome focus: the choice is between incidence alone (the usual approach), or on taking account of both incidence and related outcomes. We have developed here a bivariate approach to jointly modelling epidemic outcomes, using priors that link parameters between outcomes. For example, here, the eventual death total K d is linked (via the Bayesian prior specification) to the eventual epidemic size K c .

is approach has been used here to study the interrelationships between incidence and mortality but can be readily extended to more outcomes, such as incidence, mortality, and hospitalisations. Forecasts from such joint outcome models enable forecasts of epidemic severity (e.g., case fatality rates, case-hospitalisation ratios, deaths-hospitalisations ratios, and effective reproduction ratios) that assist in pandemic severity assessment. Severity assessment, as well as the development of summary indicators of epidemic severity, is being recognized as an important aspect of epidemic monitoring and modelling [70, 71] . e present paper also considers the gain in applying phenomenological models to later stages of incomplete epidemics, especially (in the case of COVID-19) after lockdown relaxations, where there may be a protracted period of nonzero incidence or further more pronounced waves of infection. e above analysis has therefore applied the Poisson-log-Student to late epidemic data for the UK (up to early August 2020), when there was some evidence of an incidence upturn, but not yet a full-blown second wave. e focus of this analysis was on the effective reproduction rate and the case fatality ratio, both indicators of epidemic severity. As mentioned in [72] , "the COVID-19 pandemic has shown that the effective reproduction rate of the virus R t is a crucial determinant not only of public health, but also of public policy." ere are a number of ways of estimating this quantity, including novel approaches such as using Google mobility data [73] . Here, an analysis using estimates of R t based on a two-phase Richards model suggests an upturn in transmission in the UK by late July/early August 2020. Hence, the added value of the study is provided by the following: (a) proposing a multivariate framework readily applied for assessing severity, (b) providing an approach for assessing alternative methods for representing overdispersion, and (c) providing a methodology for assessing long-term trends and making longer-term forecasts (e.g., of case fatality and the effective reproduction rate) in a multiwave situation. 

Many applications of phenomenological models have been to complete epidemics. However, evaluation of such models based simply on their fit to observed data may give only a partial picture. Also relevant to epidemic model assessment, particularly for policy application, is the accuracy of medium-term forecasts for incomplete epidemics. Arguably, evaluation in this case is better done using a cross-validation approach, where only some of the observed data are used to estimate parameters, and a hold-out sample can then assess the accuracy of forecasts. e analysis here of UK epidemic data in the first half of 2020 has shown that fit to training data and the cross-validation fit are consistent in their choice of preferred model option.

e contribution of the present paper is to illustrate a bivariate approach to two epidemic outcomes, and how prior information (under a Bayesian approach) can be applied to interlink the parameters governing each outcome. e benefits of a bivariate (and potentially multivariate) approach include the ability to forecast severity measures such as the case fatality ratio, and the borrowing of strength over outcomes in making forecasts [26] . A further contribution is that a focus on incidence and new deaths rather than cumulative outcomes has brought into sharper focus the question of adequately representing Poisson overdispersion. e latter is caused by often erratic fluctuations in the observed series, apparent in UK COVID data on new cases and deaths. e analysis has provided new evidence on the relative fit and forecasting performance of different ways of representing Poisson overdispersion in epidemic count data, and potential gains through using heavier tail alternatives to the negative binomial. We have also considered Bayesian analysis of longer-term epidemic trends, where multiple waves may exist, and illustrated monitoring of case fatality and epidemic reproduction. e implications of the research (e.g., in planning healthcare provision) are that effective medium-term forecasts of COVID-19 incidence and mortality can be provided by the proposed methodology. Extension of the bivariate approach (e.g., to include incidence, mortality, and hospitalisations) provides scope for forecasting other indicators relevant to severity assessment. e methodology also presents a way to monitor longer-term infection numbers leading to early detection of incipient upturns in infection numbers, via continuous monitoring and forecasting of the effective reproduction ratio. Assessing whether the latter is confined to values below 1 is important for strategic epidemic containment. e present study has some possible limitations. Comparative analysis of alternative ways of representing overdispersion has been limited here to UK data, and to the first wave COVID-19 epidemic, and will require validation with other epidemic time series, both for COVID-19 and other infectious diseases. A caveat, although not a limitation per se, is that, with regard to Bayesian estimation, relatively informative priors may be needed to guarantee stable estimation and ensure convergence. For example, diffuse gamma priors on the eventual epidemic size parameter caused convergence problems in the analysis reported here.

Regarding future research, as just pointed out, the methodology for comparing overdispersion approaches should be assessed with other epidemic datasets. Other Poisson mixtures may be considered such as the Poisson log skew-normal [74] or Poisson mixtures with other densities [75, 76] . Other types of analysis regarding forecasting potential can also be envisaged, such as forecast combinations, for instance, forecasts based on combining the logistic, Richards, and Gompertz curves [77, 78] . It would also be useful, especially in planning hospital care capacity, to apply a bivariate approach to COVID-19 incidence and hospitalisations, or a trivariate approach to incidence, mortality, and hospitalisations. A trivariate approach will provide model estimates and forecasts of indicators central to severity assessment.

e data used to support the findings of the study were obtained from the European Centre for Disease Prevention and Control (https://www.ecdc.europa.eu/en/publicationsdata/download-todays-data-geographic-distribution-covid-19-cases-worldwide).

e author declares no conflicts of interest.

An interactive web-based dashboard to track COVID-19 in real time

Estimating epidemic exponential growth rate and basic reproduction number

Hospital capacity and operations in the coronavirus disease 2019 (COVID-19) pandemic-planning for the nth patient

Time series modelling to forecast the confirmed and recovered cases of COVID-19

Centre for Disease Control and Prevention (CDC), Pandemic Severity Assessment Framework (PSAF), Centre for Disease Control and Prevention (CDC)

Novel framework for assessing epidemiologic effects of influenza epidemics and pandemics

Comparative analysis of phenomenological growth models applied to epidemic outbreaks

A unified approach to the Richards-model family for use in growth analyses: why we need only two model forms

Why is it difficult to accurately predict the COVID-19 epidemic?

Simple framework for real-time forecast in a data-limited situation: the Zika virus (ZIKV) outbreaks in Brazil from 2015 to 2016 as an example

A comparative analysis between a SIRD compartmental model and the Richards growth model

Applying the unified models of ecology to forecast epidemics, with application to COVID-19

Mathematical models to characterize early epidemic growth: a review

Forecasting COVID-19 impact on hospital bed-days, ICU-days, ventilator days and deaths by US state in the next 4 months

Real-time forecast of multiphase outbreak

Foot and mouth disease revisited: re-analysis using Bayesian spatial susceptible-infectious-removed models

Endemic-epidemic models with discrete-time serial interval distributions for infectious disease prediction

Impact of nonpharmaceutical interventions (NPIS) to reduce COVID-19 mortality and healthcare demand

COVID-19 virus outbreak forecasting of registered and recovered cases after sixty day lockdown in Italy: a data driven model approach

Forecasting of COVID-19 cases and deaths using ARIMA models

Forecasting the novel coronavirus COVID-19

Estimating the number of SARS-COV-2 infections and the impact of social distancing in the United States

Short-term forecasts of the COVID-19 epidemic in Guangdong and Zhejiang

Estimating and forecasting COVID-19 evolution in the UK and its regions

Modeling the incidence and death rates of COVID-19 pandemic in different regions of the world

Univariate versus multivariate time series forecasting: an application to international tourism demand

Model selection and evaluation based on emerging infectious disease data sets including A/ H1N1 and Ebola

Cholera modeling

Avoidable errors in the modelling of outbreaks of emerging pathogens, with special reference to Ebola

Mixed Poisson distributions

Estimation of the final size of the coronavirus epidemic by the logistic model

Logistic growth modelling of COVID-19 proliferation in China and its international implications

Richards model: a simple procedure for real-time prediction of outbreak severity

Richards model revisited: validation by and application to infection dynamics

When will the coronavirus pandemic peak

Prediction and analysis of coronavirus disease

Using phenomenological models to characterize transmissibility and forecast patterns and final burden of Zika epidemics

A generalizedgrowth model to characterize the early ascending phase of infectious disease outbreaks

e math of epidemic outbreaks and spread (part 3) least squares fitting of Gompertz growth models

Lognormal and gamma mixed negative binomial regression

Robust empirical Bayes analyses of event rates

A comparative study of four estimators for analyzing the random event rate of the Poisson process

Predicting future purchases with the Poisson log-normal model

Testing species abundance models: a new bootstrap approach applied to Indo-Pacific coral reefs

Why estimating population-based case fatality rates during epidemics may be misleading

Getting less deadly?

A Bayesian formulation of exploratory data analysis and goodness-of-fit testing

A new framework and software to estimate time-varying reproduction numbers during epidemics

e effective reproduction number as a prelude to statistical estimation of timedependent epidemic trends

Robust mixture modelling using the t distribution

A Bayesian analysis of two-piece distributions based on the scale mixtures of normal family

Regression analysis of hierarchical Poisson-like event rate data: superpopulation model effect on predictions

Widely applicable Bayesian information criterion

WinBUGS-a Bayesian modelling framework: concepts, structure, and extensibility

General methods for monitoring convergence of iterative simulations

Epidemiology and transmission of COVID-19 in 391 cases and 1286 of their close contacts in Shenzhen, China: a retrospective cohort study

e early phase of the COVID-19 outbreak in Lombardy

Transmission dynamics, serial interval and epidemiology of COVID-19 diseases in Hong Kong under different control measures

Serial interval of novel coronavirus (COVID-19) infections

Serial interval of COVID-19 among publicly reported confirmed cases

Computing gamma distribution parameters by gamma.parms.from.quantiles(

A likelihood-based method for real-time estimation of the serial interval and reproductive number of an epidemic

e R0 package: a toolbox to estimate reproduction numbers for epidemic outbreaks

Estimating the time-varying reproduction number of SARS-CoV-2 using national and subnational case counts

Using a negative binomial regression model for early warning at the start of a hand foot mouth disease epidemic in Dalian

e IDEA model: a single equation approach to the Ebola forecasting challenge

Maximum likelihood fitting of the Poisson lognormal distribution

Methods for fitting the Poisson-lognormal distribution to microbial testing data

Assessing the severity of COVID-19

Situation analysis and an insight into assessment of pandemic COVID-19

e stochastic reproduction rate of a virus

COVID-19 data analysis, part 6: insights from forecasting COVID-19 transmission rates in Mexico with Google mobility data

Analysis of over-dispersed count data with extra zeros using the Poisson log-skewnormal distribution

Heteroscedastic nonlinear regression models using asymmetric and heavy tailed two-piece distributions

Poisson-beta exponential distribution: properties and applications

Forecast combination for euro area inflation: a cure in times of crisis?

Real-time forecasting of epidemic trajectories using computational dynamic ensembles