key: cord-0463429-04s1hcic
authors: Levene, Mark
title: A skew logistic distribution with application to modelling COVID-19 epidemic waves
date: 2022-01-31
journal: nan
DOI: nan
sha: 2ba4bd85daa39b8b32589ca882a0b65bb757332b
doc_id: 463429
cord_uid: 04s1hcic

A novel yet simple extension of the symmetric logistic distribution is proposed by introducing a skewness parameter. It is shown how the three parameters of the ensuing skew logistic distribution may be estimated using maximum likelihood. The skew logistic distribution is then extended to the skew bi-logistic distribution to allow the modelling of multiple waves in epidemic time series data. The proposed skew-logistic model is validated on COVID-19 data from the UK, and is evaluated for goodness-of-fit against the logistic and normal distributions using the recently formulated empirical survival Jensen-Shannon divergence (${cal E}SJS$) and the Kolmogorov-Smirnov two-sample test statistic ($KS2$). We employ 95% bootstrap confidence intervals to assess the improvement in goodness-of-fit of the skew logistic distribution over the other distributions. The obtained confidence intervals for the ${cal E}SJS$ are narrower than those for the $KS2$ on using this data set, implying that the ${cal E}SJS$ is more powerful than the $KS2$.

In exponential growth the population grows at a rate proportional to its current size. This is unrealistic, since in reality growth will not exceed some maximum, called its carrying capacity. The logistic equation [Bac11,  Chapter 6] deals with this problem by ensuring that the growth rate of the population decreases once the population reaches its carrying capacity [Pan14] . Statistical modelling of the logistic equation's growth and decay is accomplished with the logistic distribution [JKB95] [ Kri15, Chapter 22] , noting that the tails of the logistic distribution are heavier than those of the ubiquitous normal distribution. The normal and logistic distributions are both symmetric, however, real data often exhibits skewness [Das10] , which has given rise to extensions of the normal distribution to accommodate for skewness, as in the skew normal [AC14] and epsilon skew normal [MH00] distributions. Subsequently, skew logistic distributions were also devised, as in [Nad09, SB16] .

Epidemics, such as COVID-19, are traditionally modelled by compartmental models such as the SIR (Susceptible-Infected-Removed) model and its extension the SEIR (Susceptible-Exposed-Infected-Removed) model, which estimate the trajectory of an epidemic [Li18] . These models formulate the solution to the maximum likelihood estimation of the parameters of the skew logistic distribution. In Section 4, we make use of an extension of the skew logistic distribution to the bi-skew logistic distribution to model a time series of COVID-19 data items having more than a single wave. In Section 5 we provide analysis of daily COVID-19 deaths in the UK from 30/01/20 to 30/07/21, assuming the skew logistic distribution as an underlying model of the data. The evaluation of goodness-of-fit of the skew logistic distribution to the data makes use of the recently formulated ESJS, and compares the results to those when employing the KS2 instead. We observe that the same technique, which we applied to the analysis of COVID-19 deaths, can be used to model new cases and hospitalisations. Finally, in Section 6, we present our concluding remarks. It is worth noting that in the more general setting of information modelling, being able to detect epidemic waves may help supply chains in planning increased resistance to such adverse events [SJ22] . We note that all computations were carried out using the Matlab software package.

Here we introduce a novel skew logistic distribution, which extends, in straightforward manner, the standard two parameter logistic distribution [JKB95] [Kri15, Chapter 22 ] by adding to it a skew parameter. The rationale for introducing the distribution is that, apart from its simple formulation, we believe its maximum likelihood solution, presented below is also simpler than those derived for other skew logistic distributions, such as the ones investigated in [Nad09, SB16] . This point provides further justification for our skew logistic distribution when introducing the bi-skew logistic distribution in Section 4. Now, let µ be a location parameter, s be a scale parameter and λ be a skew parameter, where s > 0 and 0 < λ < 2. Then, the probability density function of the skew logistic distribution at a value x of the random variable X, denoted as f(x;λ,µ,s), is given by

noting that for clarity we write x−µ above as a shorthand for (x−µ), and κ λ is a normalisation constant, which depends on λ. When λ = 1, the skew logistic distribution reduces to the standard logistic distribution as in [JKB95] and [Kri15] (Chapter 22), which is symmetric. On the other hand, when 0 < λ < 1, the skew logistic distribution is positively skewed, and when 1<λ<2, it is negatively skewed. So, when λ=1, κ λ =1, and, for example, when λ=0.5 or 1.5, κ λ =2/π. For simplicity, from now on, unless necessary, we will omit to mention the constant κ λ as it will not effect any of the results.

The skewness of a random variable X [Das10, Kri15] , is defined as

and thus, assuming for simplicity of exposition (due the linearity of expectations [Das10] ) that µ=0 and s=1, the skewness of the skew logistic distribution, denoted by γ(λ), is given by

First, we will show that letting λ 1 =λ, with 0<λ 1 <1, we have γ(λ 1 )>0, that is f(x;λ 1 ,0,1) is positively skewed. We can split the integral in (2) into two integrals for the negative part from −∞ to 0 and the positive part from 0 to ∞, noting that when x=0, the expression to the right of the integral is equal to 0. Then, on setting y =−x for the negative part, and y =x for the positive part, the result follows, as by algebraic manipulation it can be shown that

implying that γ(λ 1 )>0 as required.

Second, in a similar fashion to above, on letting λ 2 =λ 1 +1=λ, with 1<λ 2 <2, it follows that γ(λ 2 ) < 0, that is f(x;λ 2 ,0,1) is negatively skewed. In particular, by algebraic manipulation we have that exp(−λ 2 y)

implying that γ(λ 2 )<0 as required.

The cumulative distribution function of the skew logistic distribution at a value x of the random variable X is obtained by integrating f(x;λ,µ,s), to obtain F (x;µ,s,λ), which is given by

where 2 F 1 (a,b;c;z) is the Gauss hypergemoetric function [AS72, Chapter 15]; we assume a,b and c are positive real numbers, and that z is a real number extended outside the unit disk by analytic continuation [POP17] .

The hypergeometric function has the following integral representation [AS72, Chapter 15],

where c>b. Now, assuming without loss of generality that u=0 and s=1, we have that

where x is a real number.

Therefore, from (7) it can be verified that: (i) 2 F 1 (1,2 − λ;3 − λ;−exp(x)) is monotonically decreasing with x, (ii) as x tends to plus infinity, 2 F 1 (1,2−λ;3−λ;−exp(x)) tends to 0, and (iii) as x tends to minus infinity, 2 F 1 (1,2−λ;3−λ;−exp(x)) tends to 1, since

We now formulate the maximum likelihood estimation [WA18] of the parameters µ,s and λ of the skew logistic distribution. Let {x 1 ,x 2 ,...,x n } be a random sample of n values from the density function of the skew logistic distribution in (1). Then, the log likelihood function of its three parameters is given by

In order to solve the log likelihood function, we first partially differentiate lnL(λ,µ,s) as follows:

It is therefore implied that the maximum likelihood estimators are the solutions to the following three equations:

which can be solved numerically. We observe that the equation for µ in (10) does not contribute to solving the maximum likelihood, since the location parameter µ is equal to the mean only when λ=1. We thus look at an alternative equation for µ, which involves the mode of the skew logistic distribution.

To derive the mode of the skew logistic distribution we solve the equation,

to obtain

Thus, motivated by (12) we replace the equation for µ in (10) with

where m is the mode of the random sample.

We start by defining the bi-skew logistic distribution, which will enable us to model more than one wave of infections at a time. We then discuss how we partition the data into single waves, in a way that we can apply the maximum likelihood from the previous section to the data in a consistent manner.

We present the bi-skew logistic distribution, which is described by the sum,

of two skew logistic distributions. It is given in full as

which characterises two distinct phases of logistic growth (c.f. [Mey94, SMF04] ). We note that (14) can be readily extended to the general case of the sum of multiple skew logistic distributions, however for simplicity we only present the formula for the bi-skew logistic case. Thus while the (single) skew logistic distribution can only model one wave of infected cases (or deaths, or hospitalisations) the bi-skew logistic distribution can model two waves of infections, and in the general cases any number of waves.

In the presence of two waves the maximum likelihood solution to (14) would give us access to the necessary model parameters, and solving the general case in presence of multiple waves, when the sum in (14) may have two or more skew logistic distributions, is evidently even more challenging. Thus we simplify the solution for the multiple wave case, and concentrate on an approximation assuming a sequential time series when one wave strictly follows the next. More specifically, we assume that each wave is modelled by a single skewed logistic distribution describing the growth phase until a peak is reached, followed by a decline phase; see [CH82] who considers epidemic waves in the context of the standard logistic distribution. Thus a wave is represented by a temporal pattern of growth and decline, and the time series as whole describes several waves as they evolve.

To provide further clarification of the model, we mention that the skew-bi logistic distribution is not a mixture model per se, in which case there is a mixture weight for each distribution in the sum, as in, say, a Gaussian mixture [Bis06, Chapter 9 ]. In the bi-skew logistic distribution case we do not have mixture weights, rather we have two phases or in our context waves, which are sequential in nature, possibly with some overlap, as can be seen in Figure 1 (c.f. [Mey94, SMF04] ); however, strictly speaking, the bi-skew logistic distribution can be viewed as a mixture model where the mixture weights are each 0.5 and a scaling factor of 2 is applied. Thus, as an approximation, we add a preprocessing step where we segment the time series into distinct waves, resulting in a considerable reduction to the complexity of the maximum likelihood estimation. We do, however, remark that the maximum likelihood estimation for the bi-skew logistic distribution is much simpler than that of a corresponding mixture model, due to the absence of mixture weights. In particular, although we could, in principle make use of the EM (expectation-maximisation) algorithm [RW84] [Bis06, Chapter 9] to approximate the maximum likelihood estimates of the parameters, this would not be strictly necessary in the bi-skew logistic case, cf. [Mac21] . The only caveat, which holds independently of whether the EM algorithm is deployed or not, is the additional number of parameters present in the equations being solved. We leave this investigation as future work, and focus on our approximation, which does not require the solution to the maximum likelihood of (14); the details of the preprocessing heuristic we apply are given in the following section.

Here we provide a full analysis of COVID-19 deaths in the UK from 30/01/20 to 30/07/21, employing the ESJS goodness-of-fit statistic and comparing it to the KS2 statistic. The daily UK COVID-19 data we used was obtained from [GOV21] .

As a proof of concept of the modelling capability of the skew logistic distribution, we now provide a detailed analysis of the time series of COVID-19 deaths in the UK from 30/01/20 to 30/07/21.

To separate the waves we first smoothed the raw data using a moving average with a centred sliding window of 7 days. We then applied a simple heuristic, where we identified all the minima in the time series and defined a wave as a consecutive portion of the time, of at least 72 days, with the endpoints of each wave being local minima apart from the first wave which starts from day 0. The resulting four waves in the time series are shown in Figure 1 ; see last column of Table 1 for the endpoints of the four waves. It would be worthwhile, as future work, to investigate other heuristics, which may for example allow overlap between the waves to obtain more accurate start and end points and to distribute the number of cases between the waves when there is overlap between them.

In Table 1 we show the parameters resulting from maximum likelihood fits of the skew logistic distribution to the four waves. Figure 2 shows histograms of the four COVID-19 waves, each overlaid with the curve of the maximum likelihood fit of the skew logistic distribution to the data. Pearson's moment and median skewness coefficients [DS11] for the four waves are recorded in Table 2 . It can be seen that the correlation between these and 1−λ is close to 1, as we would expect. We now turn to the evaluation of goodness-of-fit using the ESJS (empirical survival Jensen-Shannon divergence) [LK21, Lev21] , which generalises the Jensen-Shannon divergence [Lin91] to survival functions, and the well-known KS2 (Kolmogorov-Smirnov two-sample test statistic) Table 2 : Pearson's moment and median skewness coefficients for the four waves, and the correlation between 1−λ and these coefficients.

[GC21, Section 6.3]. We will also employ 95% bootstrap confidence intervals [ET93] to measure the improvement in the ESJS and KS2, goodness-of-fit measures, of the skew-logistic over the logistic and normal distributions, respectively. For completeness we formally define the ESJS and KS2.

To set the scene we assume a time series [CX19] , x={x 1 ,x 2 ,...,x n }, where x t , for t=1,2,...,n is a value indexed by time, t, in our case modelling the number of daily COVID-19 deaths. We are, in particular, interested in the marginal distribution of x, which we suppose comes from an underlying parametric continuous distribution D. 

where I is the indicator function. In the following we will let P (z)= S(x)[z] stand for the empirical survival function S(x)[z], where the time series x is assumed to be understood from context; we will generally be interested in the empirical survival function P , which we suppose arises from the survival function P of the parametric continuous distribution D, mentioned above.

The empirical survival Jensen-Shannon divergence (ESJS) between two empirical survival functions, Q 1 and Q 2 arising from the survival functions Q 1 and Q 2 , is given by

where

We note that the ESJS is bounded and can thus be normalised, so it is natural to assume its values are between 0 and 1; in particular, when Q 1 = Q 2 its value is zero. Moreover, its square root is a metric [NV15] , cf. [LK21] .

The Kolmogorov-Smirnov two-sample test statistic between Q 1 and Q 2 as above, is given by

where max is the maximum function, and |v| is the absolute value of a number v. We note that KS2 is bounded between 0 and 1, and is also a metric. For a parametric continuous distribution D, we let φ = φ(D, P ) be the parameters that are obtained from fitting D to the empirical survival function, P , using maximum likelihood estimation. In addition, we let P φ = S φ (x) be the survival function of x, for D with parameters φ. Thus, the empirical survival Jensen-Shannon divergence and the Kolmogorov-Smirnov two-sample test statistic, between P and P φ , are given by ESJS( P ,P φ ) and KS2( P ,P φ ), respectively, where P and P φ are omitted below as they will be understood from context. These values provide us with two measures of goodness-of-fit for how well D with parameters φ is fitted to x [Lev21] .

We are now ready to present the results of the evaluation. In Table 3 we show the ESJS values for the four waves and the said improvements, while in Table 4 we show the corresponding KS2 values and improvements. In all cases the skew logistic is a preferred model over both the logistic and normal distributions, justifying the addition of a skewness parameter as can be see in Figure 2 . Moreover, in all but one case is the logistic distribution preferred over the normal distribution; this is for wave 3, where the KS2 statistic of the normal distribution is smaller than that of the logistic distribution. We observe that, for the second wave, the ESJS and KS2 values for the skew logistic and logistic distribution are the closest, since as can be seen from Table 1 Table 3 : ESJS values for the skew logistic (SL), logistic (Logit) and normal (Norm) distributions, and the improvement percentage of the skew logistic over the logistic (SL-Logit) and normal (SL-Norm) distributions, respectively.

In Tables 5 and 6 we present the bootstrap 95% confidence intervals of the ESJS and KS2 improvements, respectively, using the percentile method, while in Tables 7 and 8 we provide the 95% confidence intervals of the ESJS and KS2 improvements, respectively, using the bias-corrected and accelerated (BCa) method [ET93] , which adjusts the confidence intervals for bias and skewness in the empirical bootstrap distribution. In all cases the mean of the bootstrap samples is above zero with a very tight standard deviation. As noted above the second wave is more or less symmetric, so we expect that the standard logistic distribution will provide a fit to the data which as good as the skew logistic fit. It is thus not surprising that in this case the improvement percentages are, generally, not significant. In addition, the improvements for the third wave are also, generally, not Table 4 : KS2 values for the skew logistic (SL), logistic (Logit) and normal (Norm) distributions, and the improvement percentage of the skew logistic over the logistic (SL-Logit) and normal (SL-Norm) distributions, respectively.

significant, which may be due to the starting point of the third wave, given our heuristic, being close to its peak; see Figure 1 . We observe that, for this data set, it is not clear whether deploying the BCa method yields a significant advantage over simply deploying of the percentile method.

In Table 9 we show the mean and standard deviation statistics of the confidence interval widths, of the metrics we used to compare the distributions, implying that, in general, the ESJS goodnessof-fit measure is more powerful than the KS2 goodness-of-fit measure. This is based on the known result that statistical tests using measures resulting in smaller confidence intervals are normally considered to be more powerful, implying that a smaller sample size may be deployed [Liu13] . Table 5 : Results from the percentile method for the confidence interval of the difference of the ESJS between the logistic (Logit) and skew logistic (SL), and between the normal (Norm) and skew logistic (SL) distributions, respectively; Diff, LB, UB, CI, Mean and STD stand for difference, lower bound, upper bound, confidence interval, mean of samples and standard deviation of samples, respectively.

As mentioned in the introduction, we obtained comparable results to the above when modelling epidemic waves with the epsilon skew normal distribution [MH00] as opposed to using the skew logistic distribution; see also [KN15] for a comparison of a skew logistic and skew normal distributions in the context of insurance loss data, showing that the skew logistic performed better than the skew normal distribution for fitting the data sets tested. Further to the note in the introduction, that the skew logistic distribution is a more natural one to deploy in this case due to its heavier tails, we observe that in an epidemic scenario the number of cases counted can only be non-negative, while the epsilon skew normal also supports negative values. Table 7 : Results from the BCa method for the confidence interval of the difference of the ESJS between the logistic (Logit) and skew logistic (SL), and between the normal (Norm) and skew logistic (SL) distributions, respectively; Diff, LB, UB, CI, Mean and STD stand for difference, lower bound, upper bound, confidence interval, mean of samples and standard deviation of samples, respectively.

We have proposed the skew-logistic and bi-logistic distributions as models for single and multiple epidemic waves, respectively. The model is a simple extension of the symmetric logistic distribution, which can readily be deployed in the presence of skewed data that exhibits growth and decay. We provided validation for the proposed model using the ESJS as a goodness-of-fit statistic, showing that it is a good fit to COVID-19 data in UK and more powerful than the alternative KS2 statistic. As future work, we could use the model to compare the progression of multiple waves across different countries, extending the work of [DCDW20] . Table 9 : Mean and standard deviation (STD) statistics for the confidence interval (CI) widths using the percentile (P) and BCa methods.

Institute Of Mathematical Statistics Monographs

Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables

A Short History of Mathematical Population Dynamics

Pattern Recognition and Machine Learning. Information Science and Statistics

The Kullback-Leibler divergence as an estimator of the statistical properties of CMB maps

Methods for the measurement of epidemic velocity from time-series data

Power Analysis: An introduction for the life sciences

The Analysis of Time Series: An Introduction with R. Text in Statistical Science

Fundamentals of Probability: A First Course. Springer Texts in Statistics

The scale and dynamics of COVID-19 epidemics across europe

On an Sir epidemic model for the COVID-19 pandemic and the logistic equation. Discrete Dynamics in Nature and Society

Effects of non-pharmaceutical interventions on COVID-19 cases, deaths, and demand for hospital services in the UK: A modelling study

Measuring skewness A forgotten statistic

An Introduction to the

A bi-logistic growth model for conference registration with an early bird deadline

Nonparametric Statistical Inference. Marcel Dekker

COVID-19) in the UK, Download data

Tracking the mutant: Forecasting and nowcasting COVID-19 in the UK in 2021

Forecasting for COVID-19 has failed

Wiley Series in Probability and Mathematical Statistics

A comparison between skew-logistic and skew-normal distributions

Handbook of Statistical Distributions with Applications

A hypothesis test for the goodness-of-fit of the marginal distribution of a time series with application to stablecoin data

An Introduction to Mathematical Modeling of Infectious Diseases. Mathematics of Planet Earth

Divergence measures based on the Shannon entropy

Comparing sample size requirements for significance tests and confidence intervals

Empirical survival Jensen-Shannon divergence as a goodness-of-fit measure for maximum likelihood estimation and curve fitting

Is EM really necessary here? Examples where it seems simpler not to use EM

Bi-logistic growth

The epsilon-skew-normal distribution for analyzing near-normal data

The skew logistic distribution

Non-parametric Jensen-Shannon divergence

Growth Curve Modeling: Theory and Applications

Logistic equation and COVID-19

Numerical methods for the computation of the confluent and Gauss hypergeometric functions

Estimation of COVID-19 dynamics "on a back-of-envelope": Does the simplest SIR model provide quantitative parameters and predictions?

Mixture densities, maximum likelihood and the Em algorithm

A new skew logistic distribution: Properties and applications

The synthesis model as a planning tool for effective supply chains resistant to adverse events

Bi-phasic growth patterns in rice

An evaluation of R 2 as an inadequate measure for nonlinear models in pharmacological and biochemical research: a Monte Carlo approach

Maximum Likelihood for Social Science: Strategies for Analysis

A second wave? What do people mean by COVID waves? -a working definition of epidemic waves