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. These 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] . There is also a growing literature providing evidence on the parameters of the COVID-19 infection (for example, case fatality ratios, serial intervals, etc). Hence the potential occurs for epidemic models that are applicable to routinely collected data, that 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.
Commonly applied approaches to quantitative modelling of aggregate epidemic data differ in the data inputs they require, their assumptions, their estimability, and in their scope for practical application to making forecasts. Commonly applied methods include phenomenological growth models [2] , such as the Richards family of densities [3] , and variants of the susceptible-infected-recovered (SIR) compartment model [4] . Phenomenological models are parameterised in terms of epidemic trajectories and provide estimates of crucial epidemiological parameters [5] , while avoiding the complexity of more formal models of disease transmission which can be difficult to estimate and provide projections [4, 6] . As mentioned by Wu et al [7] , phenomenological models are "particularly suitable when significant uncertainty clouds the epidemiology of an infectious disease". By contrast, as noted in [8] , 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. They may also be complex to specify when an epidemic has more than one phase.
Here we focus on a practical forecasting approach using routinely available data on new cases and deaths (from ECDC) to estimate parameters of a bivariate version of the Reynolds model. Implementation is based on Bayesian inference principles, incorporating accumulated evidence on relevant parameters via informative priors. The 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 20-day ahead forecasts of cases and deaths, based on mid-epidemic data. Other policy relevant parameters such as the effective reproduction ratio are also estimated in analysis of longer term epidemic data.
A number of studies, with widely differing methodologies, have been made 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 [9] . This 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 [8] . This 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 have also been made, using ARIMA models [10, 11] , exponential smoothing [12] , or autoregression in cases; for example, the study by Johndron et al [13] postulates daily deaths as a lagged function of earlier new cases. Autoregression in cases, with potential for short term forecasting, is illustrated (for foot and mouth disease) by the first order autoregressive model of Lawson et al [14] , while (for multiple spatial units) the model of Bracher and Held [15] specifies a first order autoregression based on the mean incidence in adjacent areas.
Regarding phenomenological models, the study by Zhao et al [5] exemplifies their application to epidemic forecasting, though to the Zika epidemic in 2015 rather than the COVID-19 epidemic. Applications of such models to COVID-19 include [16] and [17] .
Basic phenomenological models for epidemic trajectories include the logistic, Gompertz, and Rosenzweig, which have been the basis for a range of generalizations [2, 3] . Application of the logistic model to COVID-19 is exemplified by the studies of Batista [18] and Shen [19] . 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. The Richards model [20] modifies the logistic incidence function to
The parameter a  0 modifies the incidence decline phase of the logistic, that is measures the extent of deviation from the standard logistic curve. The turning point  , when incidence peaks, is obtained when Ct equals K1  a −1/a [21] . The peak incidence is important for the healthcare planning, for example aligning the forecast peak with hospital bed capacity [22] .
Other commonly used models are the Gompertz model [23] with C ′ t  rCt logK/Ct, while the Rosenzweig [24] has
The incidence function represented by C ′ t can be used to define mean incidence in statistical likelihoods for new cases data. Thus time series of incidence counts can often be satisfactorily modelled as a Poisson, with means defined by C ′ t functions [25, 26] . Similarly the cumulative cases function Ct can be used to define mean epidemic size in models for cumulative case counts [24, 27] .
It may be preferable to estimate epidemic model parameters by analyzing new cases rather than cumulative cases, as the latter have serially correlated measurement error, leading to understatement of parameter uncertainty [28] . As stated in [28] , "independence of sequential measurement errors, ... is clearly violated when observations are accumulated through time". However, estimation using new cases and deaths 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, whereas cumulative cases and deaths are usually smooth functions.
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 [28] . Other overdispersed versions of the Poisson can be achieved by mixing the Poisson with a suitable density (e.g. a lognormal density) [29, 30] , and this may be beneficial in detecting unusual observations.
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 we have
. . , T. In practice many epidemic datasets are overdispersed relative to the Poisson, and many studies adopt a negative binomial model instead. We can specify an overdispersed model (including the negative binomial) by introducing multiplicative random effects [29, eqn 2] , such that the Poisson means for incidence are specified by
where the  t are positive random effects. For the Poisson-gamma model (equivalent to a Negative Binomial), the  t are gamma distributed with mean 1, namely
where 1/ is the overdispersion parameter mentioned by [28] . Note that the assumed parameterisation of the gamma density with random variable x is
Other options in (3) are to take u t  log t  as normally distributed [29] u t  Normal0,   2 , or Student t distributed [31] ,
where  is a degrees of freedom parameter. These two options define the Poisson-Lognormal (PLN) and Poisson-log-Student (PLS) options respectively [30] . The PLN and PLS representations may provide a more robust alternative to the Poisson-gamma [31, 32, 33, 34, 35] , 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. Thus denote d t and D t as new and cumulative deaths at time t. The 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 Ga1,  or G,  , with  small -was found to lead to convergence problems. As noted in [18] ,"... in the early stage, the logistic curve follows an exponential growth curve, so the estimation of K is practically impossible". This difficulty persists when an epidemic is past its peak but early in a downturn. However, Batista [18] 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 c e  as mean, and a suitable variance, such that the prior is still relatively diffuse. For example, suppose K c e is 250000, and the variance in the lognormal is set at 1. Then the 97.5 percentile for the lognormal prior is 1.73 million.
In the bivariate specification (for new cases and 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#the-current-case-fatality-rate-of-covid-19.
Alternatively one may link K c and K d using both the point estimator K c e and a case fatality ratio, namely K d e  K c e . Then a lognormal density prior can be taken for
logK c e   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 [36] , in which case the prior on  may be constrained to exceed the final  t based on observed data. With regard to COVID-19 not all epidemics fit this pattern, with the US (for example) showing a decline in CFRs at later epidemic stages. There is also evidence that the infection to mortality ratio (a more precise measure than the CFR) has fallen [37] .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.
Thus suppose the training sample is formed by observations up to time M  T, while Underprediction means under-forecasting of future cases and may lead to incorrect inferences regarding epidemic control, as it implies a lessening in incidence earlier than actually occurred.
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. Thus in the UK the choice on whether or not to relax the initial COVID lockdown restructions 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. The reproduction rate may also become especially relevant at later epidemic stages (post-lockdown), 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, with new cases remaining non-neglible, the concern is especially that there may be a resurgence in cases, and possibly also deaths, at some point. This is colloquially known as a "second wave", and in some European countries (e.g. France, Spain) there have been pronounced second waves in the COVID-19 epidemic. Such a resurgence indicates use of a multiphase model [39] , with a second phenomenological model applied to data after a switch-point between epidemic regimes.
Consider a model for new cases only. The two phenomenological models (before and after the switch-point) have distinct parameters. The switch-point is taken as a latent parameter, , such that for a two-phase Richards model
where IA  1 when condition A is true, and 0 otherwise. The parameter  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 switch points later for deaths than cases, due to delays in mortality upturns following incidence upturns. A three phase model would have two switch-points, with mean
An estimator of the effective reproduction rate R t at time t is based on predictions from the phenomenological model, possibly multiphase, and from an estimate of the serial interval density. The serial interval is the time between symptom onset in an infected subject and symptom onset in the infectee. The serial interval density can discretized in the form of weights  j , applied to serial interval lengths (in days) up to a maximum J . These can be used to estimate effective reproduction ratios R t within a phenomenological model to analyse new cases; see the papers [25] , [40] and [41] .
Thus R t can be estimated as
where c new,t are predicted new case data from the phenomenological model.
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. Figures 1a  and 1b 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, though with a more protracted decline as opposed to the initial upturn. Figures 2a and 2b show the relatively smooth evolution of cumulative cases and deaths. 
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 precision 1/  2 in the PLN and PLS options, and for the parameter  in the Poisson-gamma model, a gamma prior Gamma1, 0. 001 is assumed. For the PLS option, we take   4 as a preset option. The degrees of freedom can be difficult to estimate for relatively small datasets, and the option of the preset value   4 is a robust option [43, 44] . For the maximum cases (epidemic size) parameter K c we assume as lognormal prior, centred at logK c e , as discussed in section 3.2.
For the maximum deaths parameter K d , we assume K d  K c where the prior for  is a beta with mean defined by the EU-wide case fatality  EU with total prior count C set at 5, C  5 . Thus 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.
Fit to the observed training data is assessed using the Watanabe--Akaike information criterion (WAIC) [45] . Cross-validation fit is assessed by the probabilities of overprediction  c and  d , and by predictive coverage: whether the credible intervals for predicted cases and deaths (averages over each validation period) include the actual averages. Estimation uses the BUGS program [46] , with posterior estimates based on the second halves of two chain runs of 100,000 iterations, and convergence assessed using Brooks-Gelman-Rubin criteria [47] .
To provide a longer term perspective on epidemic containment, we apply the best performing model to the full set of observations as at 8th August 2020 (T=190) when the first epidemic peak had passed. At this point there were 309 thousand cumulative cases and 46500 cumulative deaths. New cases and deaths (as 7 day moving averages) had reached maxima of 4850 and 950 respectively when the UK epidemic peaked in April. As a result of lockdown measures imposed in late March, by July daily new cases averaged just over 500 daily. 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 ), though deaths continued to fall, averaging around 50 per day in early August.
To model the full time series, and since there is no upturn in deaths, we focus on new cases only. A two phase model is applied with shift parameter  assigned an exponential 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 2 we assume K 2  K 1 with  assigned an exponential prior with mean 1. Of policy interest here is epidemic containment, as summarised by the effective reproduction ratio: specifically is this ratio consistently below 1, and its 95% credible interval also entirely below 1.
To provide estimates of the effective reproduction ratio, we use accumulated evidence on COVID serial intervals from five studies [48, 49, 50, 51, 52] . 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 [53] . 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 density has shape 1.38 and rate 0.36. This density is then converted to a discretized form (with 16 bins) to provide an informative prior on the SI to the model of White and Pagano [54] , which updates the prior serial interval density using new case data for the UK. The paper [55] 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] . The updated mean serial interval is estimated as 3.5 days with standard deviation 3.1. The discretized serial interval is estimated as in [54] , with J  16 .
5.1 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. 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 replicate data at each iteration. These are based on sampling new data from the Poisson means (4) and (5), and from the Richards model parameter profile at each particular iteration.
The 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 .
Term Forecasts. 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 . The 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. This reflects the protracted nature of the UK downturn in cases. 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. 
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). The 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. The estimates hover around 1 throughout June and the first half of July, but in late July/early August tend to exceed 1. The impression from this is that success in fully effective containment is in doubt. The 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 [56] .
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 assess the accuracy of forecasts. The analysis here has shown that fit to training data and the cross validation fit are consistent in their choice of preferred model option.
Many applications of phenomenological models have also been to cumulative cases data. However, to quote [57] on cholera epidemics, "The [incidence] data from Zimbabwe show multiple peaks and other features characteristic of heterogeneously mixed populations [...]. The practice of fitting epidemic models to cumulative incidence curves rather than incidence curves can obscure these features while also violating statistical assumptions of independence between fitted data points". A focus on incidence rather than cumulative incidence does therefore bring into sharper focus the question of adequately representing often erratic fluctuations in the observed series, apparent in UK COVID data on new cases and deaths.
For larger counts, and overdispersed count data, the literature has mainly focussed on the negative binomial as the solution to such variability. So far as the authors are aware, the literature does not contain any assessment of alternative distributional choices to the Poisson-gamma mixture (negative binomial). Hence the contribution of this paper rests on a comparison of the negative binomial against alternative Poisson mixture models [30] . The analysis suggests these may usefully be considered as alternatives to the negative binomial, and may both improve fit to actual observations and provide more accurate forecast performance.
There is also value in applying phenomenological models to later stages of incomplete epidemics, especially (in the case of COVID-19) after lockdown relaxation where there may be second waves of infection. The above analysis has therefore applied the Poisson-log Student to late epidemic data for the UK,with a focus on the effective reproduction rate or effective reproduction ratio. As mentioned in [58] , "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". There are a number of ways of estimating this quantity, including novel approaches such as using Google mobility data [59] . Here an analysis using estimates of R t based on a two-phase Richards model suggest a upturn in transmission in the UK.
Competing Interests: The authors declare no competing interests. 