Dynamic optimization of a two-stage fractional system in microbial batch process

1.   Introduction

Since the onset of the Covid-19 pandemic in early 2020 an enormous research effort has been underway on the modelling and prediction of various aspects of the pandemic. For accessible reviews concentrating on general features of this modelling, see, for example, Vespignani et al. (2020), Poletto et al. (2020) and Gnanvi et al. (2021), while for discussion of the growth models widely used for predicting Covid-19 infections and deaths, see Tovissodé et al. (2020) and Shen (2020). Central to this modelling, the analysis of the mass of time series data accumulating daily and weekly from the coronavirus pandemic has become ever more important as the pandemic has progressed through its numerous phases. Spiegelhalter and Masters (2021) provide an accessible introduction to such data issues, paying particular attention to the evidence emerging from the U.K.

It is becoming increasingly apparent that econometric techniques are particularly suited to analysing this data: see, for example, Li and Linton (2020), Manski and Molinari (2020) and the review by Dolton (2021). Much of the research using these techniques has focused on short-term forecasting of cases, hospital admissions and deaths, with notable examples being Doornik et al. (2020), Doornik et al. (2021) and Harvey et al. (2021). It has also been directed at generalising, to stochastic settings, compartmental epidemiological models, such as the well-known "susceptible (S), infected (I) and recovered or deceased (R)", or SIR, model, as in Korolev (2020) and Pesaran and Yang (2021).

The focus of the present paper is rather different, however, in that we investigate the changing dynamic relationship between infections, hospital admissions and deaths using daily data from England. Section 2 thus considers the relationship between hospital admissions and subsequent deaths and introduces two models that might be useful for this task: the recently proposed balanced growth model of Harvey (2020) and the more familiar autoregressive distributed lag/error correction model used widely to analyse economic time series (see, for example, Banerjee et al., 1993, for detailed development and Mills, 2019, chapters 12 and 14, for a more introductory treatment). Section 3 extends the analysis to examining the prior relationship between infections and hospital admissions, while Section 4 links the two sets of models together before discussing the advantages and disadvantages of the two modelling procedures.

2.   Modelling the relationship between hospital admissions and deaths in England

Figure 1 shows daily hospital admissions and deaths in England between 19th March 2020 and 31st October 2021.¹ Both admissions and deaths show pronounced multiple wave patterns with admissions obviously leading deaths, but the shifting nature of the relationship between the two series is clearly discernible. How might this evolving and dynamic relationship be modelled? Attention is focused in this paper on two approaches: balanced growth modelling and the use of autoregressive distributed lags.

¹The focus here is on data from England as U.K.—wide hospital admissions rely on different definitions across the home nations.

2.1.   Balanced growth modelling

Let daily deaths due to Covid-19 in England be denoted ${y}_{t}$, $t = \mathrm{1, 2}, \cdots, T$, with their cumulation being ${Y}_{t} = {\sum }_{j = 1}^{t}{y}_{j}$, so that the growth rate of daily deaths is ${g}_{y, t} = {y}_{t}/{Y}_{t-1}$. Similarly, denote daily hospital admissions due to Covid-19 by ${x}_{t}$, their cumulation by ${X}_{t} = {\sum }_{j = 1}^{t}{x}_{j}$, and their growth rate by ${g}_{x, t} = {x}_{t}/{X}_{t-1}$.

Following Harvey and Kattuman (2020), we initially assume that there is balanced growth between daily deaths and hospital admissions lagged $k$ days, which implies the regression model.

where ${\varepsilon }_{t}$ is an error term assumed to independently and identically distributed through time with zero mean and variance ${\sigma }_{\varepsilon }^{2}$, which is denoted ${\varepsilon }_{t}\sim$ IID$\left(0, {\sigma }_{\varepsilon }^{2}\right)$. The "equilibrium" relationship between the two growth rates is given by

Between daily deaths and admissions, ${y}_{t}$ and ${x}_{t}$, the equilibrium is then

Allowing for a lag structure in the leading admissions series in (1) gives

where $h < k$ and $\sum {\beta }_{j} = 1$, a restriction that may be imposed by rewriting (4) as

so that ${\beta }_{k} = 1-{\sum }_{j = h}^{k-1}{\beta }_{j}$, a restriction that ensures that there is indeed balanced growth. The corresponding equilibrium relationship is then

where ${\bar{g}}_{x, t-k}$ is the weighted geometric mean of ${g}_{x, t-h}, \cdots, {g}_{x, t-k}$. The levels equilibrium is thus

where ${\bar{x}}_{t-k}$ and ${\bar{X}}_{t-k-1}$ are the corresponding weighted geometric means of ${x}_{t-h}, \cdots, {x}_{t-k}$ and ${X}_{t-1-h}, \cdots, {X}_{t-1-k}$, respectively. Thus $\mathrm{\Delta }$ measures the long-run response of deaths to an increase in hospital admissions: if daily admissions increase by 100 then deaths will increase by 100$\mathrm{\Delta }$ after $k$ days.

When the two series are not on the same growth path, the model can be extended by replacing the intercept $\delta$ with a stochastic trend:

or

where ${\delta }_{t}$ is defined as

i.e., ${\delta }_{t}$ is a random walk with a drift that is itself potentially a random walk. If ${\sigma }_{\zeta }^{2} = 0$ then ${\gamma }_{t} = {\gamma }_{t-1}$ and the drift is constant. On the other hand, if ${\sigma }_{\eta }^{2} = 0$ then ${\delta }_{t} = 2{\delta }_{t-1}-{\delta }_{t-2}-{\zeta }_{t-1}$ and ${\delta }_{t}$ will tend to evolve very smoothly, being known as an integrated random walk, or IRW. The equilibrium relationship in (6) is

so that the dynamic relationship between the two growth rates is given by $\mathrm{e}\mathrm{x}\mathrm{p}\left({\delta }_{t}\right)$. In terms of daily deaths and admissions, ${y}_{t}$ and ${x}_{t}$, we have

Thus, an increase of 100 in hospital admissions will lead to an increase of ${100\mathrm{\Delta }}_{t}$ deaths in the following $k$ days and this long run response will shift through time. As Harvey (2020) shows, this model may be arrived at by assuming that deaths and admissions follow Gompertz processes separated by $k$ days, but such an assumption is not necessary.

Equation (9) may be fitted by casting it into state space form and employing the Kalman filter. Estimation is carried out by maximum likelihood using the predictive error decomposition and the estimates of ${\delta }_{t}$ used to compute the ${\mathrm{\Delta }}_{t}$ series shown in Figures 3 and 7 below are the "smoothed" estimates, obtained by running the Kalman filter first forwards from $t = k+1$ to $t = T$ and then backwards from $t = T$ to $t = k+1$. Mills (2019, chapter 17) provides an introductory discussion to state space modelling and Harvey (1989) is the classic exposition. A value for the lag $k$ must be obtained before estimation can be undertaken. This lag could be selected by considering the basic model (9) with $h = 0$, setting $k$ to an initial value ${k}_{max}$ and then running the regression (9) for $k = {k}_{max}, {k}_{max}-1, {k}_{max}-2, \cdots$, the sequence stopping when ${\widehat{\beta }}_{k}$ is significant at some pre-chosen level of significance. Alternatively, clinical considerations may suggest an appropriate value for $k$ and indeed $h$, and this simpler approach is also investigated. A further refinement may be to select only a subset of the lagged regressors $\mathrm{l}\mathrm{o}\mathrm{g}\left({g}_{x, t}\right), \cdots, \mathrm{l}\mathrm{o}\mathrm{g}\left({g}_{x, t-k}\right)$, thus determining $h$ and including only significant lags. This may be done, for example, by using a stepwise least squares algorithm or other sequential testing procedure. The balanced growth assumption may be checked by including $\mathrm{l}\mathrm{o}\mathrm{g}\left({g}_{x, t-k}\right)$ (or, indeed, any other lag) in (9) and testing for its significance.

2.2.   ARDL modelling

An autoregressive distributed lag (ARDL) model linking deaths and admissions may be specified in general as

where ${u}_{t}\sim$ IID$\left({0, \sigma }_{u}^{2}\right)$ is an error term. Typically, the lag lengths $m$ and $n$ will be unknown and must be determined from the data. An algebraically equivalent but often more convenient form of this ARDL$\left(m, n\right)$ model, particularly for model specification and inference, is the error correction model (ECM)²

²Banerjee et al (1993) sets out and analyses in detail the algebraic equivalencies existing between the ARDL and ECM formulations. It should be emphasised that the recasting of (14) as (15) is a purely algebraic transformation and is not predicated on any particular statistical properties of the data. With integrated data, cointegration leads from an ARDL to an ECM via Granger's representation theorem (see Engle and Granger, 1987). The series here are not integrated, as may be demonstrated from standard unit root tests, and so we are in a stationary world in which the ECM (15) is simply a more convenient representation for our purposes than the ARDL (14).

where $\nabla$ is the difference operator defined such that $\nabla {y}_{t} = {y}_{t}-{y}_{t-1}$ and where the coefficients of (14) and (15) are linked by the set of relationships

The error correction term is ${y}_{t-1}-d{x}_{t-1}$, which embodies the long-run, or equilibrium, relationship $y = dx$ between deaths and admissions, with $d$ being termed the long-run or total multiplier. An increase in hospital admissions of 100, say, will eventually increase deaths by $100d$.

The speed at which the increase in deaths approaches the total multiplier depends on $c$, the speed of adjustment parameter. The smaller is $c$, the faster the speed of adjustment and the quicker the total multiplier is arrived at. The actual time path of adjustment depends upon the lag coefficients ${\psi }_{i}$ in the distributed lag

where

This result is most easily obtained by utilising the lag operator $B$, defined such that ${B}^{j}{\mathcal{z}}_{t}\equiv {\mathcal{z}}_{t-j}$ (note that the difference operator introduced in (15) may then be written as $\nabla = 1-B$). This allows (14), on ignoring the error term as we are only interested in the systematic dynamics here, to be written as

where the lag polynomials $\phi \left(B\right)$ and $\theta \left(B\right)$ are defined as

Equation (19) can then be expressed as

where

Thus, the lag coefficients in (18) are obtained by equating coefficients of powers of $B$ in $\psi \left(B\right)\phi \left(B\right) = \theta \left(B\right)$. The total multiplier is then given by the sum of these lag coefficients, i.e., $d = {\sum }_{i = 1}^{\infty }{\psi }_{i}$ and the increase in deaths after $l$ days is given by the $l^{{\rm{th}}}$ interim multiplier ${d}_{l} = {\sum }_{i = 1}^{l}{\psi }_{i}$ $l = 1, 2, \cdots$ which will converge to $d$ as $l$ increases, i.e., ${d}_{l}\to d$ as $l\to \infty$.

2.3.   Fitting balanced growth and ARDL/ECM models to the relationship between hospital admissions and deaths in England

The logarithms of the growth rates of daily hospital admissions and subsequent deaths between 21^st March 2020 and 31^st October 2021 are shown in Figure 2. The series are evidently not on the same growth path so that a model of the form (9) rather than (5) is clearly required. Both clinical considerations and exploratory sequential testing along the lines suggested in section 2.1 above suggest setting $h = 0$, thus allowing for the death of some patients on the day of their admission to hospital, and $k = 7$, so that there is a one week delay between hospital admission and death.

Column (1) of Table 1 reports estimates of (9) with $\theta \mathrm{l}\mathrm{o}\mathrm{g}\left({g}_{x, t-7}\right)$ included as an additional term. If $\theta$ is non-zero then balanced growth does not hold: with this term included,

Several of the coefficients are estimated to be insignificantly different from zero, including $\theta$, so that balanced growth is confirmed. Column (2) of Table 1 reports the estimates of the model with ${{\beta }_{0}, \beta }_{1}, {\beta }_{3}, {\beta }_{5}$ and $\theta$ all restricted to be zero on using a 10% significance level; the remaining coefficients are all significantly positive (${\beta }_{7}$ is calculated to be $0.267\left(0.057\right)$). It is also found that ${\sigma }_{\eta }^{2} = 0$, so that ${\delta }_{t}$ will be a smoothly evolving IRW. Figure 3 shows the resulting estimate of ${\mathrm{\Delta }}_{t}$, the response of deaths to an increase in hospital admissions, calculated using the geometric mean

and it is indeed seen to evolve smoothly. 95% confidence interval upper and lower bounds for ${\mathrm{\Delta }}_{t}$ are also shown, and these indicate that ${\mathrm{\Delta }}_{t}$ is indeed estimated precisely.

${\mathrm{\Delta }}_{t}$ reaches a maximum at 0.38 on 27^th April 2020, just under three weeks after the peak in daily deaths. There are further local maxima of 0.28 on 5^th December 2020 and 0.31 on 26^th December 2020, both occurring during the second wave of deaths in late 2020/early 2021. ${\mathrm{\Delta }}_{t}$ reaches a minimum of 0.06 on 9^th June 2021, the previous minimum being 0.13 on 10^th August 2020, and at the end of October 2021 it stood at 0.15. This suggests that, while ${\mathrm{\Delta }}_{t}$ is positively related to the number of daily deaths, its magnitude has been declining relative to the number of deaths as clinical practice has improved over the course of the pandemic and the vaccination program has been rolled out. The long run response of deaths to an increase of 100 in hospital admissions has fallen from a maximum of 37 deaths (95% confidence interval $\left(34, 40\right)$) during the first wave of the pandemic in April 2020, to around 30 $\left(27, 33\right)$ during the second wave at the turn of the year and down to just 6 $\left(5, 7\right)$ by the summer of 2021, although this had increased to 15 $\left(13, 17\right)$ by the end of October.

ARDL/ECM models were developed to characterise the relationship between daily deaths and admissions for the complete sample period to 31^st October 2021 and for two sub-periods, the first ending on 31^st January 2021, the second beginning on 1^st February 2021, the "break point" being chosen to reflect the increasing vaccination uptake at this point. The models, estimated by nonlinear least squares, were initially selected using the AIC information criteria with insignificant coefficients (at the 10% level) then being sequentially removed to obtain a parsimonious ECM specification. Details of the chosen specifications are given in Table 2. HAC standard errors are shown in parentheses to accommodate any remaining autocorrelation and heteroskedasticity in the residuals. The most notable features of the models are the much lower estimate of the equilibrium parameter $d$ in the second sub-period than in the first, which is accompanied by a much larger estimate of $c$. Note that the estimates of these parameters are all highly significant and, in terms of goodness of fit, the model for the second sub-period has much the superior performance, with a higher ${R}^{2}$ and a lower equation standard error than the model for the first sub-period.

The time paths of the interim multipliers ${d}_{l}$ for the two sub-periods are shown in Figure 4. The paths have been smoothed to ensure that they are monotonically non-declining, since the number of deaths that follow a given increase in hospital admissions cannot fall with time! With $\widehat{d} = 0.328\left(0.035\right)$ in the first sub-period ending on 31^st January 2021, an increase of 100 in hospital admissions will therefore eventually lead to a further 33 deaths (95% confidence interval $\left(26, 40\right)$). From Figure 4 it is seen that approximately $95\%$ of these deaths $\left(31\right)$ occur within $14$ days of admission to hospital.

In contrast, for the second sub-period beginning on 1^st February 2021, $\widehat{d} = 0.164\left(0.010\right)$, so that now an increase of 100 in hospital admissions eventually leads to just 16 further deaths (95% confidence interval $\left(14, 18\right)$). Moreover, only 69% of these deaths $\left(11\right)$ occur within $14$ days of admission to hospital. It would thus appear that, as the vaccination programme was rolled out along with other improvements in clinical practice, so the hospital death rate was more than halved with patients being kept alive for longer.

3.   Modelling the relationship between Covid-19 cases and hospital admissions

As well as the linkage from hospital admissions to deaths, there is also the prior link from testing positive (i.e., becoming infected) for Covid-19 to admission into hospital. This has become particularly important to analyse since the roll-out of the vaccination program in the U.K: is there evidence that vaccination has "broken the link between cases and admissions", as has been stated several times by the government? To investigate whether this is indeed the case, we analyse English data on daily positive cases and hospital admissions from 1^st September 2020 to 31^st October 2021, as shown in Figure 5. Earlier data have been excluded both because of the limited extent of testing during the early months of the pandemic and the hiatus in cases and admissions during the summer months of 2020.

We first fit balanced growth models to the logarithms of hospital growth, ${g}_{x, t}$, and case growth, ${g}_{w, t} = {w}_{t}/{W}_{t-1}$, where ${w}_{t}$ and ${W}_{t}$ are daily and cumulative cases, respectively. These logarithms are shown in Figure 6. From both clinical considerations and exploratory sequential testing a lag of $k = 14$ was chosen, along with $h = 1$.

Estimates of the balanced growth Equation (9) are shown in Table 3 after the elimination of insignificant coefficients. Here ${\sigma }_{\eta } > 0$ so that, with both ${\sigma }_{\zeta }$ and ${\gamma }_{t}$ estimated to be very small, ${\nabla \delta }_{t}\approx {\eta }_{t}$, a driftless random walk. The consequent ${\mathrm{\Delta }}_{t}$ series, estimated under the assumption of balanced growth (which is questionable here since $\theta$ appears to be significantly negative: $\widehat{\theta } = -0.321\left(0.064\right)$), is shown in Figure 7, along with 95% confidence interval upper and lower bounds. These show that ${\mathrm{\Delta }}_{t}$ is estimated extremely precisely. The trend in ${\mathrm{\Delta }}_{t}$ values is generally upwards until the middle of February 2021, reaching a maximum of 0.12 (95% confidence interval $\left(11, 13\right)$) on 12^th December 2020 (i.e., an additional 120 hospital admissions result from an increase of 1000 positive test cases), after which it turns down, reaching a minimum of 0.02 $\left(0.019, 0.021\right)$ on 21^st July 2021 (20 additional admissions from an increase of 1000 cases), and was still around this value at the end of October. This suggests that there is indeed evidence that the rollout of the vaccination program has had an impact on the relationship between cases and subsequent hospital admissions.

ARDL/ECM models were fitted to the sub-periods 1^st September 2020–31^st January 2021 and 1^st February–31^st October 2021, with the resulting estimates and statistics shown in Table 4. The estimate of the long-run parameter $d$ is 0.087 (0.006) in the early period but just 0.027 (0.001) in the later period, i.e., up to the end of January 2021 an increase of 1000 positive test cases ultimately led to 87 (95% confidence interval $\left(75, 99\right)$) additional hospital admissions, whereas from February such an increase led to only an additional 27 (95% confidence interval $\left(25, 29\right))$ admissions. The estimate of the speed of adjustment parameter, $c$, on the other hand, has remained almost the same across the two periods, just increasing from 0.103 to 0.118. Again, the fit of the second sub-period model is vastly superior to that of the first.

The interim multipliers for admissions are plotted in Figure 8. In the period up to the end of January 2021, approximately $85\%$ of hospital admissions occur within $28$ days of a positive test. In contrast, for the second sub-period beginning on 1^st February 2021, only $79\%$ of these admissions occur within $28$ days of a positive test. It would thus appear to be clear that, as the vaccination programme was rolled out, along with other improvements in medical practice, so fewer people were admitted to hospital.

4.   Conclusions

Given the "causal structure" inherent in the relationship between cases, hospital admissions and deaths, the models in Sections 2 and 3 may be linked together. If the equilibrium responses from the balanced growth models are now denoted ${\mathrm{\Delta }}_{y, t}$ and ${\mathrm{\Delta }}_{x, t}$, to use an obvious extension of notation, then the response of deaths to an increase in cases is given by the product ${\mathrm{\Delta }}_{y, t}\times {\mathrm{\Delta }}_{x, t-7}$ to ensure an appropriate timing match-up, and this series is shown in Figure 9.³ This product reached a maximum of 33 deaths per thousand cases on 18^th December 2020 (95% confidence interval $\left(28, 38\right))$ and by the end of October 2021 had declined to just 3 $\left(2, 4\right)$.

³Using (7) we have ${y}_{t} = {\Delta }_{y, t}{x}_{t-7}$ and ${x}_{t} = {\Delta }_{x, t}{w}_{t-14}$, so that ${y}_{t} = {\Delta }_{y, t}{\Delta }_{x, t-7}{w}_{t-21}$.

Similarly, denoting ${d}_{x}$ and ${d}_{y}$ as the total multipliers from the ARDL/ECM models, with ${d}_{x, l}$ and ${d}_{y, l}$ being the accompanying interim multipliers, then the products ${d}_{x}\times {d}_{y}$ and ${d}_{x, l}\times {d}_{y, l}$ provide the total multiplier and set of interim multipliers for the response of deaths to an increase in cases, these being shown for the two sub-periods in Figure 10. The total multipliers are 28 deaths per thousand cases (95% confidence interval $\left(18, 38\right)$) for the first sub-period and 4 $\left(3, 5\right)$ for the second, both consistent with the balanced growth estimates.

Both modelling approaches provide a consistent story: as clinical practice has evolved and the vaccination program rolled out, so there has been a substantial decline in hospital admissions and subsequent deaths for a given level of infection. Are there grounds for preferring one approach over the other? Each have their advantages and disadvantages.

Balanced growth models have the advantage of providing a time varying equilibrium response which can therefore be tracked through time. They have the disadvantages of requiring balanced growth to hold, which might not be the case, and a fixed adjustment period dependent on the setting of the lag $k$. Their specification is also very precise with little flexibility in the setup and, moreover, estimation requires specialised software. ARDL/ECM models, on the other hand, have greater flexibility in their specification and may be estimated by routine regression software. The adjustment to the long-rum response is freely estimated rather than fixed but the multiplier is not time varying, with evolving relationships having to be investigated through fitting the models over sub-periods of the data, which have to be selected by the investigator or by using break-point procedures and tests, such as those developed by Bai and Perron (1998). Given the limited time period available and the clear waves in the data, we have chosen not to follow this latter route here. However, as more data becomes available, and particularly given the more recent wave of the pandemic associated with the omicron variant, formal testing for break-points would be a useful extension of the ECM modelling approach. Given these competing benefits and drawbacks, it would thus seem sensible to keep an open mind and to use both models to track the behaviour of key Covid variables through time.

While this paper has focused on the relationship between infections, hospital admissions and deaths in England for the period of the pandemic up to the end of October 2021, it is clear that the models may be used for similar data from other countries and time periods. They may also be used for data on infections, admissions and deaths disaggregated into age groups and regions, if such data are available, as it is for England, where such research is ongoing.

Conflict of interest

The author declares that there are no conflicts of interest in this paper.