Introduction

Imported cases from China contributed to the first COVID-19 wave in Malaysia from January 25, 2020 to February 26, 2020 [1]. This first wave had a total of 22 cases out of which 20 were directly linked to foreign travel while the remaining two cases were local transmissions [2, 3]. The second wave of COVID-19 in Malaysia was largely attributed to a four-day mass gathering held in Sri Petaling from February 27, 2020, which contributed to an exponential rise of COVID-19 cases in the country. This gathering involved over 16, 000 participants including a large number of foreign visitors from countries that later registered COVID-19 cases; see, for example, [4] which refers to the official Facebook link of the Ministry of Health, Malaysia [5]. On March 18, 2020, the Malaysian government introduced a nationwide lockdown which was the Phase 1 Movement Control Order (MCO) throughout the country to stem the spread of COVID-19. Phase 1 MCO was enforced for a 2-week period starting from March 18, 2020 to March 31, 2020. During this phase, various non-pharmaceutical interventions (NPIs) were strictly enforced by means of movement restrictions, wearing of face masks, social distancing and hand hygiene practices to reduce disease transmission. The Phase I MCO was then evaluated after 2 weeks based on case trends and model forecasts by the Ministry of Health (MOH) Malaysia. As local transmission persisted, the MCO was extended a total of four times until May 12, 2020. Phases 2, 3 and 4 of the MCO covered periods starting from April 1, 2020 until May 12, 2020 (Phase 2—April 1, 2020 to April 14, 2020, Phase 3—April 15th to April 30, 2020, and Phase 4—May 1, 2020 to May 12, 2020). Subsequently, from May 4, 2020, the MCO was eased into the Conditional MCO (CMCO) until June 9, 2020. However, during CMCO, there were still several identified hot-spots of COVID-19 which were placed under Enhanced MCO (EMCO) with the aforementioned strict movement control restrictions.

The implementation of MCO proved to be effective—the reported COVID-19 cases reached its peak around the first week of April and subsequently started to reduce. However, concerns remained whether a rebound in transmission would occur when the MCO was lifted and if compliance to NPIs were not followed strictly at that time. In order to gain a quantitative understanding of the effect of MCO, we developed a class of mathematical models to capture the dynamics of COVID-19 spread before and after the MCO implementation. A variant of the Susceptible-Exposed-Infected-Recovered (SEIR) model is proposed and developed in this paper which incorporates heterogeneity in the transmission dynamics, additional compartments for asymptomatic transmission and a change-point, chosen adaptively based on data, to reflect the shift in spread dynamics after the MCO implementation. The models developed in this paper are able provide a quantitative assessment of the extent of COVID-19 spread during the pre-MCO (large gathering) and MCO periods by means of a measure of infectivity developed from them. This measure is similar to the basic reproduction number, commonly denoted by ℛ₀, but can be calculated for more complex epidemic models such as the ones proposed here.

Deterministic compartmental models, such as the Susceptible-Infected-Recovered (SIR) or the Susceptible-Exposed-Infectious-Recovered (SEIR) models, provide a good theoretical framework to study infectious disease spread, and have been widely used and reported in the literature. However, more complex versions of these models, and their stochastic counterparts require data-rich inputs to model all aspects of the disease dynamics. Data-rich inputs, if lacking, can be compensated using reliable and informative prior elicitation. Considering the acute nature and scale of the pandemic as well as the urgent need for a multisectorial response, comprehensive data availability of the pandemic was limited in Malaysia. For example, the open source website outbreak.my initially reported a transmission network for all cases; however, it was unable to cope with the scale of the pandemic when it intensified. Factoring in this data limitation, we choose to develop models that are deterministic, rather than stochastic, while ensuring that they are able to capture salient transmission dynamics satisfactorily. As mentioned earlier, we enhance the deterministic models by incorporating compartments for asymptomatic transmission and a change-point to reflect the shift in disease dynamics. We also take into account heterogeneity in the disease spread such as varying contacts among susceptibles within the closed population. The starting point of our proposed models are the class of epidemic models with power transmission dynamics which are shown to incorporate heterogeneity (see [6, 7]).

Several studies in the literature [8–12] have analyzed the effects of NPIs in reducing the number of COVID-19 cases. In [8], the effects of different types of NPIs on COVID-19 cases are modeled using a negative-binomial distribution whose underlying parameters incorporate country information, type of NPI implemented and change-point effects. The associated Bayesian analysis is carried out using Markov Chain Monte Carlo (MCMC) algorithms to arrive at posterior parameter estimates and credible intervals. No epidemic models are considered in this work. A generalization of the SEIR epidemic model is considered in [10] to understand the dynamics of transmission in New York, USA, under various NPI settings. However, the model is complex and requires data-rich inputs for the estimation of all unknown parameters. As a result, the authors derive baseline epidemiological parameters from published literature and not from actual observed cases in New York, and in the end conduct a simulation study based on the assumed parameter values. The study in [9] extends the work of [13] and computes a time-varying basic reproduction number as a way of gauging the effect of NPIs over time. Both these works assume that serial intervals (i.e., the time between onset of symptoms for the infector and the infectee) can be computed for each case, which is another situation requiring data-rich input.

Studies that use compartmental epidemic models as a way of gauging the time-varying effects of NPIs have emerged over the course of the pandemic [14–16]. Compartmental epidemic models naturally model disease spread via contact rates which directly quantify the extent of NPIs (since, as mentioned earlier, NPIs are designed to reduce person-to-person transmission). Thus, epidemic models provide a natural approach for considering time-varying effects of the MCO period. Further, in this paper, the estimation of SEIR parameters is carried out based on local considerations and local data; they are not obtained from published literature based on studies conducted elsewhere where their local dynamics can be vastly different.

We seek to address one important aspect of Malaysia’s multifaceted response to the COVID-19 pandemic, that is, to inform the health officials at MOH and aid them in their decision-making. Thus, our model was developed under local considerations using local data. Our model and related analyses are able to provide a quantitative assessment of (1) the impact of the Sri Petaling gathering, and (2) the extent of decreasing transmission during the MCO period by incorporating a time-varying contact rate parameter, which is estimated using locally available data. In essence, the proposed models here are being used as a lens to interpret the observed data in terms of when, and to what extent, a reduction in COVID-19 transmission occurred as result of the implementation of MCO.

Materials and methods

Data collection

Daily situation reports on COVID-19 cases in Malaysia are published by the National Crisis Preparedness and Response Centre (CPRC) of MOH, as well as other official websites (such as outbreak.my). The data on daily COVID-19 cases have been published since 21 January 2020 and are publicly available. The reports consist of confirmed daily and cumulative cases, recovered cases and deaths, as well as cases requiring ICU care and ventilator support. Cases by states are also available for the 13 states and 3 federal territories. In this study, we studied characteristics of the second COVID-19 wave in Malaysia starting from March 1, 2020 (corresponding to the final day of the Sri Petaling gathering). Data used for the current study are confirmed daily cases for Malaysia, and for two states: Selangor and Sarawak. These states were chosen to illustrate the aggressive transmission propagated by the Sri Petaling gathering. Selangor is the state where Sri Petaling is located and from where a majority of the participants originated, whereas Sarawak represents a state which was essentially not affected by this gathering. The time period of study is between March 1, 2020 (end of Sri Petaling gathering) and April 28, 2020, covering the period immediately after the Sri Petaling gathering and the first three phases of the MCO. Our study duration is further divided into two periods. The first period ranges from March 1, 2020 until March 18, 2020, which covers the subsequent 17 days after the gathering. The second period is taken from March 18, 2020 until April 28, 2020, covering the three successive Phases 1, 2 and 3 of the MCO. Fig 1 gives the trajectories of reported daily COVID-19 cases between March 1, 2020 and April 28, 2020 for Malaysia, and the states of Selangor and Sarawak.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Reported daily cases.

Reported daily cases for (a) Malaysia, (b) Selangor and (c) Sarawak. The time period considered is from March 1, 2020 to April 28, 2020.

https://doi.org/10.1371/journal.pone.0252136.g001

All COVID-19 cases reported by MOH were confirmed by real-time reverse transcriptase-polymerase chain reaction (real-time RT-PCR) tests. A positive case was reported when the person in question was found to be positive for SARS-CoV-2 via a real-time RT-PCR test. Upon confirmation, the individual was isolated at COVID-19 designated hospitals and healthcare facilities [17]. Active cases are defined as infected persons who were currently undergoing treatment, and hence, isolated and removed; the individual is assumed to be unable to infect other susceptibles in the population, and hence “removed” from further modelling steps. This study did not consider transmission from positive isolated COVID-19 patients to health personnel as there was no evidence of this type of transmission occurring in the COVID-19 designated hospitals in Malaysia.

The modified SEIR model

The aforementioned SEIR model does not account for heterogeneity, asymptomatic transmission and change points. To this end, we propose models with power transmission dynamics that incorporate heterogeneity in the disease parameters; see, for example, [6, 7]. Further, the infectious compartment in (3) of the SEIR model is now split into two, I_o and I_c, for symptomatic (or, observed) and asymptomatic (or, cryptic) individuals, who, respectively, exhibit and do not (or, mildly) exhibit symptoms but are nevertheless infectious. Correspondingly, the R compartment is also split into two, as mentioned earlier, to accommodate quarantined and un-quarantined cases. It is assumed that the proportion of individuals transiting from E to I_o is p. The remaining exposed individuals (a proportion of 1 − p) transition into the I_c compartment and remain undetected throughout their disease experience.

From now on, we consider only renormalized state values which represent proportions, and not actual numbers, of the population. The modified SEIR model with asymptomatic and observed infectiousness is given by the following system of ODEs: (6) (7) (8) (9) (10) (11) where the rate of new infections is now given by (12) as opposed to βS(t)I(t) in (5) for the SEIR model. Fig 2 gives the flow chart of the modified SEIR model and its associated parameters.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Flow chart of the modified SEIR model.

Compartments of the modified SEIR model and their associated parameters.

https://doi.org/10.1371/journal.pone.0252136.g002

A key difference of the model formulation in (6)–(12) is the power transmission dynamics used to model heterogeneity in the population; see [6, 7]. The epidemic models developed in [6, 7] incorporate parametric heterogeneity, that is, these models describe a heterogeneous population where individuals in the population have different trait values from one person to another. In this framework, disease characteristics related to the population, such as the likelihood of a successful transmission for a susceptible or the rate of infectivity of an infectious individual, can be taken to vary heterogeneously. This heterogeneity is characterized via probability distributions, and is known as distributed susceptibility and infectivity in [6, 7]. When a susceptible and infectious trait, respectively, has a Gamma(k₁, ν₁) and Gamma(k₂, ν₂) distribution (where Gamma(k, ν) refers to the Gamma distribution with shape parameter k > 0 and scale parameter ν > 0), it is shown in [6, 7] that the resulting epidemic model has power law transmission characteristics with the rate of incidence (of new infections) given by (13) where the powers on S and I are q₁ = 1 + 1/k₁ and q₂ = 1 + 1/k₂ depending on the shape parameter of the corresponding Gamma distribution assumed for distributed susceptibility and infectivity, respectively. Since k₁, k₂ > 0, it follows that q₁, q₂ > 1. The values of q₁ = 1 and q₂ = 1 corresponding to (5) can be thought as the limiting case when k₁ → ∞ and k₂ → ∞ representing homogeneity. From this discussion, it follows that the powers v, w_o and w_c on S(t), I_o(t) and I_c(t), respectively, are such that v ≥ 1, w_o ≥ 1 and w_c ≥ 1. The lower bounds on v, w_o and w_c recover the original SEIR model dynamics with no heterogeneity. The remaining parameters have the following interpretation: (1) α represents a small yet significant force of infection that starts the local infection process but is eventually overwhelmed by it. In our context, α represents the initial force of infection arising from, say, foreign infectious individuals attending the large gathering at Sri Petaling starting February 27th. (2) The parameters γ_o and γ_c have the same interpretation as γ in (3), that is, they are rates of infectiousness but for the observed and asymptomatic compartments, respectively. (3) The parameter δ is the same as before, namely, it is the rate of incubation associated with the exposed compartment. (4) The functions β_o(t) and β_c(t) are time varying transmission rates for the observed and asymptomatic compartments, respectively, in the modified SEIR model with β_c(t) = μβ₀(t) and μ ∈ [0, 1]. In other words, we assume that the transmission rate for asymptomatic individuals is smaller than that of symptomatic individuals; this is a plausible assumption to make as asymptomatic individuals generally possess a lower viral load which leads to lower chances of a successful transmission. On the other hand, a longer asymptomatic infectious period may compensate for the lower transmission rates for an asymptomatic individual, and this possibility is captured by the model via the parameter γ_c.

A change-point is incorporated into the modified SEIR model to capture the shift in disease dynamics before and after the start of the MCO. Several works in the COVID literature have incorporated change-points to study the effectiveness of interventions; see, for example, [24–26]. To incorporate our change-point, an unknown threshold date, T*, is chosen so that observed daily cases fall either to the left or right of T*. We denote the observation window to the left of T* by consisting of dates from March 1st up to and including T*. The window to the right of T* is denoted by which consists of dates after T* up to and including April 28th, 2020. The change-point date T* is inferred from data; it is not taken as March 18th, 2020 which is the date of the start of MCO. Choosing T* in this data-driven way is justified as the impact of MCO on observed cases may not be realized immediately. The modified SEIR model with change-point T* is governed by the ODE system (6)–(12) for all time points . For the parameters of the modified SEIR model in , we denote them using the same symbol as before but with a subscript L. For , the same ODE system (6)–(12) is considered but now with a different set of parameter values in compared to which are denoted using a subscript R.

The transmission rates in are assumed to be constant, that is, and β_c,L = μ_L β_o,L. The transmission rate for the I_o compartment in is taken as (14) to model possible changes in disease transmission over time. The functional form of β_o,R(t) in (14) has an initial value of at t = T* after which it decreases (provided a₁ < 0) to the asymptotic value of a₂ γ_o,R as t → ∞. Thus, a₂ γ_o,R represents the residual disease transmission that may be present even during the MCO period, for example, due to close contact between family members in the same household. The general functional form of β_o,R(t) subsumes the constant disease transmission rate model as a special case by taking a₁ = a₂ = 0 in (14). The constant rate submodel has the advantage of not explicitly assuming any functional form for the change in disease transmission over time. On the other hand, it can only ascertain if there is an overall change (increase or decrease) in transmission after the change point T*. Similar to the relationship β_c,L = μ_L β_o,L in , we assume β_c,R(t) = μ_R β_o,R(t) for for a possibly different μ value.

In the subsequent text, several redundancies are removed based on the interpretation of the parameters involved. These simplications are reasonable and lead to a reduction in the number of unknown parameters. This reduction, in turn, improves the model fitting and inference procedures. It is reasonable to assume that γ_c,L = γ_c,R ≡ γ_c, δ_L = δ_R ≡ δ, μ_L = μ_R ≡ μ and p_L = p_R ≡ p; in other words, intervention does not change the values of these parameter in and . These parameters can be assumed constant since they relate to intrinsic characteristics of the disease and of the underlying population. The parameter γ_c relates to asymptomatic infections which are not detected in both and , δ relates to the incubation period of the disease, μ refers to the transmissibility ratio of asymptomatic to symptomatic individuals (which depends on the viral load of infectees only), and p refers to the proportion of symptomatic versus asymptomatic infectees. Hence, these parameters are characteristics of the overall population and the disease, and not the intervention. However, we do not assume that γ_o,L = γ_o,R for the symptomatic infectious periods, but rather 1/γ_o,R ≤ 1/γ_o,L since the process of detecting and isolating symptomatic individuals during MCO was generally quicker and more efficient in compared to , resulting in a shorter symptomatic infectious period in compared to .

In the Results section, the constant transmission rate submodel is first fitted to observed data to determine if the MCO is succesful in reducing the overall disease transmission. Only when this is established, the full change-point model, with the explicit decay form of (14) and a₁ < 0, is fitted to the observed data. Further, in the model formulation of (6)–(11), only the R_o compartment is modelled directly using a likelihood function based on daily observed cases. The other compartments of the modified SEIR model remain latent and do not have any direct observation processes for modelling based on likelihoods; see the Bayesian Inference section for further details.

Bayesian inference

Model fitting and inference is carried out in a Bayesian framework.

The likelihood.

The total numbers of daily cases, , which constitute noisy observations from the R_o-compartment in the modified SEIR model, are assumed to be distributed as (17) where NB(x;ν, τ) denotes the negative binomial probability function with mean ν and variance ν + τν². The parameter τ measures overdispersion with respect to the Poisson probability function which is recovered from NB(x;ν, τ) when τ = 0; see, for example, [30, 31]. The observed data on daily cases in Malaysia exhibited significant overdispersion, and as a result, the Poisson likelihood did not fit the observed cases well. More discussion on this aspect is presented later.

Assignment of priors.

Several parameters of the modified SEIR model in (6)–(12) take different values in and , and reflect changes in the disease trajectory due to intervention. Other parameters, including those that relate to the intrinsic characteristics of the disease, remain constant throughout and . Defining for H ∈ {L, R}, we denote the collection of all parameters of the full change-point model in by where all notation used for the parameters is as before. In the subsequent discussion, separate priors are elicited either on individual components of Θ, or on specific subsets of Θ where the parameters within these subsets satisfy certain restrictions. A full prior on Θ then follows from assuming independence between the separate prior assignments. For each ξ₁ ∈ {p, μ}, the prior distribution on ξ₁ is assumed to be uniform, with hyperparameters and , representing the lower and upper bounds, respectively. Since each ξ₁ represents a proportion, we have and ; however, the choices of hyperparameters used in the experiments for all parameters, including ξ₁ ∈ {p, μ}, will be discussed later. For now, we provide only the explicit forms of the priors assumed on the components of Θ.

The parameters that relate to the infectious periods, γ_o,L, γ_o,R and γ_c, satisfy the relations 1/γ_o,R ≤ 1/γ_o,L (as pointed out earlier) and 1/γ_o,L < 1/γ_c. The latter condition follows from the fact that the asymptomatic infectious period, being undetected, will generally be longer than its symptomatic counterpart in , since symptomatic individuals will be isolated once they test positive. The joint prior distribution on {γ_o,L, γ_o,R, γ_c} can be specified according to the following generation scheme: We first generate , then set 1/γ_o,L = 1/γ_o,R + ξ₂, and finally take 1/γ_c = 1/γ_o,L + ξ₃ where with and , independently for j = 2 and 3. We refer to the infectious periods reported in literature, namely [18–20, 32], to choose reasonable values for , , and . Based on the reported values, infectious periods within 2–8 days are found to provide a reasonable range for our experiments. The prior on δ, where 1/δ is the incubation period, is taken as 1/δ ∼ U(A_δ, B_δ). To select reasonable values for the lower and upper bounds, we refer to the reported literature (for example, see [33]) and consider an incubation period between 2–13 days to be a satisfactory range for our study.

The proposed Bayesian methodology calls for Θ-samples to be generated from its full prior specification as described above. For this purpose, explicit mathematical experssions of the joint prior density function on parameter subsets of Θ, such as {γ_o,L, γ_o,R, γ_c}, are not required. One only needs to be able to simulate the reciprocals based on the generation scheme described above, and then invert the reciprocals to obtain prior samples of γ_o,L, γ_o,R and γ_c from their implicit joint density specification.

Next, the prior on τ is taken as τ ∼ U(A_τ, B_τ) where 0 ≤ A_τ < B_τ, and B_τ is a pre-specified positive constant representing the maximum extent of overdispersion apriori. The prior on α_L is assumed to be , where and represents the maximum value of a small but constant force of infection arising, for example, from infectious foreign attendees of the Sri Petaling gathering. The parameter α_R is set to 0 as there is no such force of infection in . For the parameters w_o,L, w_o,R and v_L, independent uniform priors are considered with lower and upper bounds given by and for w_o,L, and for w_o,R, and and for v_L, respectively. Since these parameters have a natural lower bound of 1 that represents homogeneity, the hyperparameters of the uniform priors are chosen in the range [1, ∞). The priors on w_c,L, w_c,R and v_R are taken using the following relations: w_c,L = w_o,L + ξ₄, w_c,R = w_o,R + ξ₅, and v_R = v_L + ξ₆ where with and , independently of each other for j = 4, 5 and 6.

The prior on the change-point T* is taken to be uniform on dates from March 18th, 2020, to March 31st, 2020, both inclusive. Since we define T* as the number of days after March 1st 2020, T* ∼ U(17, 30) a priori.

The constant contact rate β_o,L during the pre-MCO period cannot be determined by direct observation. Hence, putting a prior directly on β_o,L is difficult. However, since , the quantity has the interpretation of a reproduction number, a more fundamental quantity for infectious diseases compared to β_o,L. Thus, one can obtain a range of values for a_0,L from previous reports on similar flu-like epidemics (see, for example, [34]) where we found the range 0–3.5 appropriate for covering the true value of a_0,L in our case. We take the prior on a_0,L as where and based on benchmark values in [34]. The function β_o,R(t) has the form in (14) which depends on coefficients a_0,R, a₁ and a₂. For the prior on a_0,R, we consider two cases: When using the constant rate submodel to determine if an overall reduction in transmission has occurred or not, the prior on a_0,R is chosen as with and . This ensures that the prior support includes both positive and negative values to represent possible increase and decrease in the transmission rate. On the other hand, if the full model of (14) is considered, a_0,R is chosen deterministically (hence, no prior) to ensure continuity of the infection process before and after T*. Specifically, a_0,R is chosen so that the incidence rate of new infections, h(S(t), I_o(t), I_c(t)), in and coincide at t = T*. In the case of the full change-point model, a₁ is given a uniform prior with support on . We take and if a reduction in disease transmission is first established by the constant rate submodel. For the prior on a₂, note that . Hence, a reasonable prior to put on a₂ is .

Next, we describe the prior assignment on the initial number of symptomatic infectious and exposed individuals, denoted by i₀ and e₀, respectively. Given δ and γ_o,L, whose priors were elicited earlier, the prior on i₀ is taken as , where the midpoint and half-wdith depend on δ and γ_o,L. Reasonable choices of the midpoint and half-width are made based on the differential equation for the R_o compartment. Noting that for from (9), we get . To obtain an estimate of , a second-order polynomial is fitted using least squares to the trajectory of cumulative cases in a window of m ≥ 3 days starting from day 0. is then estimated by taking the first derivative of the fitted polynomial and evaluating it at time t = 0. Subsequently, is set as . The half-width is fixed to be reasonably large to represent the maximum extent of prior uncertainty around . Similar to i₀, the prior on the initial number of exposed individuals, e₀, is taken as for the midpoint and half-width . The midpoint is determined similarly as before using the value of the second derivative of the fitted polynomial evaluated at time t = 0 and the generated value of δ. The half-width is once again fixed to be large to represent the maximum extent of prior uncertainty around .

The model fitting experiments were conducted either by choosing default values for the lower and upper bounds, or by choosing values guided by literature when such information is available. For all other hyperparameters where no information is available, reasonable values were chosen at the start of the experiments, and the hyperparameters were adjusted by repeated trial and error runs to achieve the best fit models to the observed data akin to the empirical Bayes approach [35], at least informally. Table 1 gives the values of the hyperparameters used to obtain the plots and inference in the Results section.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Hyperparameter values for the reported results.

https://doi.org/10.1371/journal.pone.0252136.t001

Computational algorithm.

Based on the negative binomial likelihood and prior elicitation in The Likelihood and Assignment of Priors, respectively, the posterior of Θ can be derived using Bayes theorem as (18) (19) where is the collection of observed daily cases from March 1st until April 28, 2020, L(Θ) is the entire likelihood for D and π(Θ) is the prior on Θ as described in Assignment of Priors. Bayesian inference is carried out using Monte Carlo importance sampling. A total of M samples, Θ_i for i = 1, 2, ⋯, M, are generated from the prior specification π(Θ). The likelihood, L(Θ_i), is computed for each Θ_i and normalized to obtain weights (20) for i = 1, 2, ⋯, M. The ensemble of weights and samples constitute a weighted sample from the posterior distribution π(Θ|D) and any summary posterior measure can be obtained from the ensemble via Monte Carlo. For any function g of Θ, the posterior expectation of g(Θ), is approximated by the corresponding ensemble average the latter being an unbiased estimator of E(g(Θ)) with variance proportional to 1/M which tends to 0 as M → ∞. Thus, becomes a more accurate estimate of E(g(Θ)) as M → ∞. Different choices of g are used to obtain the posterior means, variances and credible intervals of parameters in Θ.

Resampling from the ensemble converts it to a new ensemble of weights and samples where the new weights are now uniform, i.e., . The usual sample mean and sample variance of can now be used to provide an (unbiased) estimate of E(g(Θ)) and its corresponding variance. The above importance sampling methodology is easy to code and implement but one has to be cautious so that the distribution of weights is not too extreme. If the prior on Θ is chosen such that it spans a considerably larger region compared to where the likelihood concentrates its mass, the distribution of weights will be become extreme. As a result, almost all of the weights will be extremely small (close to zero) and only a few will have large positive values. Resampling these weights will select Θ_i-samples associated to these large weights multiple times and this will degrade the estimation of E(g(Θ)). Thus, it is important to judiciously choose the prior ranges, where possible, so that it adequately covers regions of high likelihood values. This is achieved in the Results section by running different selections of the hyperparameters and judiciously adjusting the range of the priors (where possible and respecting the information on the disease and intervention dynamics) to obtain the best fitting models. Based on this importance sampling method as well, an approximation to the Maximum-a-Posteriori (MAP) estimator of Θ can be found since (21) say, for large M. To see why the above approximation is valid, we partition Θ into Θ = (Ψ, T*) where Ψ consists of all continuous parameters in Θ and T* is the change point which takes on finitely many values apriori. Correspondingly, partition as . Further assume that the posterior is continuous at Ψ = Ψ_MAP and lies in the interior of supp[π(Θ|D)] where

The occurence of the event for an arbitrarily small ϵ > 0 implies that there exists a neighbourhood of which is contained in supp[π(Θ|D)] with , and none of the M samples lie in . Hence, as M → ∞. In other words, as M becomes large, the sampling based maximum will be close to the true MAP estimate with probability tending to 1. The approximation in (21) could fail if any one of the assumptions made above is not valid: The posterior could be discontinuous at Ψ = Ψ_MAP or may not lie in the interior of supp[π(Θ|D)]. However, discontinuity is not a concern here as is a continuous (in fact, differentiable) function of Ψ. To ensure that lies in the interior of supp[π(Θ|D)], the experiments were run based on judicious choices of the prior hyperparameters which has been mentioned earlier as well. The resampled values are used for plotting and providing a quantification of uncertainty around (see, for example, Figs 10 and 11 in the Results section). Numerical measures of uncertainty can be obtained, for example, by computing the usual sample variance based on to describe the variablity around the sample mean which in this case is an unbiased estimate of E(g(Θ)).

The Bayesian computational algorithm described in this section is developed using R and RStudio^®. The likelihood L(Θ) for a specific Θ is evaluated numerically from the solution of the ODE system (6)–(12) obtained using the deSolve package in R.

Results

The study period is from March 1, 2020 until April 28, 2020 with T₀ = 0 and T₁ = 59. The impacts of the Sri Petaling gathering and MCO implementation are analyzed here using the proposed model described in The modified SEIR model section. Reported daily cases at the national level as well as for Selangor and Sarawak are used for model fitting and parameter estimation. Note that all states in Malaysia implemented MCO Phases 1–3 using the same guidelines and protocols. Thus, one can gauge the impact of the Sri Petaling gathering on COVID-19 spread in Malaysia based on a comparison between states and the national experience. Here, Selangor and Sarawak are chosen as two such representative states with high and low population densities, respectively.

First, we investigate if the MCO implementation had an overall effect of reducing COVID-19 transmission rates. For this purpose, the constant rate submodel of (14) is used and the prior on a_0,R in is chosen to be uniform with support on both positive and negative values. The Bayesian inference methodology described in the Bayesian Inference section is carried out with M = 50, 000 to obtain and samples from the posterior of Θ in (19). The curves of F(t), defined as and ΔR_o(t)≡R_o(t) − R_o(t − 1) for each day t are obtained based on and are displayed in Fig 3. This submodel captures broad features (increasing and decreasing trends) of the reported cases trajectories in all three panels for Malaysia, Selangor and Sarawak. To quantify the overall change in transmission before and after MCO implementation, Δ_t (see (16)) is obtained for t in and . The plots of Δ_t versus t are shown in Fig 4. A consistent feature of the plots in all three panels is that they first increase for time points t ≤ T* followed by a significant drop for t > T*. Hence, we conclude that an exponential rise in cases occurred right after the completion of the Sri Petaling gathering on March 1st 2020, and the implementation of the MCO successfully stemmed the exponential rise at the national, Selangor and Sarawak levels.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Constant rate submodel—Daily cases.

Reported daily cases (red line), and overlay plots of ΔR_o(t) (blue line) and F(t) (green line) based on for the constant rate submodel for (a) Malaysia, (b) Selangor and (c) Sarawak.

https://doi.org/10.1371/journal.pone.0252136.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Constant rate submodel—Δ_t.

Plots of Δ_t versus t for the constant rate submodel for (a) Malaysia, (b) Selangor and (c) Sarawak.

https://doi.org/10.1371/journal.pone.0252136.g004

With reduced disease transmission established in , we next proceed to utilize the functional form (14) of β_o(t) as a quantitative model for transmission decay in . The Bayesian inference methodology in the Bayesian Inference section is applied to the full model with M = 50, 000 to obtain and samples from the posterior of Θ in (19). The curves of F(t) and ΔR_o(t)≡R_o(t) − R_o(t − 1) for each day t are obtained based on and are displayed in Fig 5. Daily cumulative cases and the curve of R_o(t) are displayed in Fig 6. We note from these figures that the proposed model captures broad features of the observed data and is an improvement over the constant rate submodel. Uncertainty estimates are obtained for all unknown parameters in Θ based on the ensemble . Variability estimates can be obtained for all parameters and their functions. As an illustration, we demonstrate the extent of variability inherent in the posterior visually for the F(t) curve (see (17)) for t ∈ [T₀, T₁]. This is displayed in Fig 10 which shows that most of the reported case numbers are well within the limits of variability of the posterior of F(t). Hence, the proposed model together with the negative binomial likelihood are able to explain the variability in the reported case numbers. However, there are a few exceptions, the most notable being the reported case number on Day 14 for Malaysia in Fig 10(a). We present an explanation for this outlying case later in the section Discussion on Reporting Delays, Case Redistribution and Overdispersed Likelihoods.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Modified SEIR model: Daily cases.

Reported daily cases (red line), and overlay plots of ΔR_o(t) (blue line) and F(t) (green line) based on for (a) Malaysia, (b) Selangor and (c) Sarawak.

https://doi.org/10.1371/journal.pone.0252136.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Modified SEIR model: Cumulative cases.

Reported cumulative cases (red line), and overlay plots of R_o(t) (green line) based on for (a) Malaysia, (b) Selangor and (c) Sarawak.

https://doi.org/10.1371/journal.pone.0252136.g006

Further results from the Bayesian analyses are summarized in Tables 2–4. These tables give the MAP estimates of parameters and their corresponding 95% credible intervals for Malaysia, Selangor and Sarawak, respectively. The correlation plots of selected pairs of parameters are provided in Figs 7 and 8. From these figures and others which we did not provide due to space limitations, we find that the posterior samples of parameter pairs are not strongly correlated with each other, except for the pairs for which correlation was induced apriori via the joint prior elicitation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. MAP estimates and associated credible intervals for Malaysia.

https://doi.org/10.1371/journal.pone.0252136.t002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. MAP estimates and associated credible intervals for Selangor.

https://doi.org/10.1371/journal.pone.0252136.t003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. MAP estimates and associated credible intervals for Sarawak.

https://doi.org/10.1371/journal.pone.0252136.t004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Correlation plots 1.

Correlation plots of selected pairs of parameters of the modified SEIR model for Malaysia (top row), Selangor (middle row) and Sarawak (bottom row).

https://doi.org/10.1371/journal.pone.0252136.g007

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Correlation plots 2.

Correlation plots of a different set of pairings of parameters of the modified SEIR model for (a) Malaysia, (b) Selangor and (c) Sarawak.

https://doi.org/10.1371/journal.pone.0252136.g008

We provide a summary of the salient findings of our inference procedure. The symptomatic and asymptomatic infectious periods as well as the incubation periods are found to be around 6–8 days for Malaysia, Selangor and Sarawak. These findings are similar to values reported in the literature for other countries; see, for example, [18–20, 32]. Change-points T* are estimated not too far away from the date of MCO implementation, March 18th 2020. For Malaysia, T* = 17 is the MAP estimate which corresponds to March 18th, 2020, and the associated 95% credible interval is (17, 23). For Selangor, the MAP estimate of T* is T* = 24 (March 25th 2020) with (20, 28) being the 95% credible interval. For Sarawak, the transition date is less precise. The MAP estimate is T* = 28 (March 29th 2020) but the 95% credible interval (19, 33) is much larger indicating higher uncertainty in T*. This can be attributed to the fact that the trajectory of reported case numbers for Sarawak shows a slower and more gradual increase, then decrease, compared to Malaysia and Selangor (see Fig 5).

The plots of Δ_t versus t for Malaysia, Selangor and Sarawak are provided in Fig 9. We note that all three panels in Fig 9 indicate a decay in the transmission rates after T*. The measure Δ_t is calculated to be approximately 3.39, 2.50 and 2.08 for Malaysia, Selangor and Sarawak, respectively, at the start of , that is, when t = T₀. Starting from t = T₀, Δ_t showed an increase in , reaching values of4.55, 3.42 and 2.50 at t = T*, respectively, for Malaysia, Selangor and Sarawak. Hence, before the MCO was introduced, Selangor achieved a rate of increase higher than that of Sarawak. After T*, Δ_t for Malaysia, Selangor and Sarawak registered a decay demonstrating the effectiveness of the MCO. Δ_t declined sharply to a value around0.33, 0.0003 and 0.31, respectively, for Malaysia, Selangor and Sarawak at t = T* + 10, and after that, it declined more gradually to its corresponding asymptote. Based on Δ_t, it is seen that the initial transmission rates tend to be higher for areas with a higher population density (comparing Selangor and Sarawak). On the other hand, based on the MAP estimates of a₁ in of −1.756 and −0.325 for Selangor and Sarawak, respectively (see Tables 3 and 4), higher population density areas also experience a faster decline in the transmission rates under an effective implementation of the MCO. Although the MAP estimate of a₁ for Malaysia (a₁ = −0.371 from Table 2) is not as negative as it should be, we will show in the next section that a redistribution of cases further improves this estimate of a₁ and brings it closer to that of Selangor.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Modified SEIR model: Δ_t.

Plots of Δ_t versus t: (a) Malaysia, (b) Selangor and (c) Sarawak.

https://doi.org/10.1371/journal.pone.0252136.g009

Discussion on reporting delays, case redistribution and overdispersed likelihoods

Fig 10 indicates the presence of outliers that fall outside the limits of variability of the posterior. The most notable outlier is the total number of new cases reported on Day 14 for Malaysia. Generally speaking, such outliers highlight a mismatch between the proposed model and the observed data, and point towards model inadequacy. However, we wish to emphasize that this is not the case here. One key consideration is the effect of delay, that is, whether or not the reported case numbers coincide with the day of testing. It is highly likely that a lag occurred in the reporting of cases since the COVID-19 experience was new to Malaysia. Based on the report [36], it is reasonable to assume that delays in testing and reporting were expected during the initial days of the COVID-19 outbreak in Malaysia. The peak on Day 14 seem to suggest a significant backlog of reporting of cases coupled by the fact that the day was a Monday, the first day of the week.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Variation of modelled daily cases.

Illustration of the variability of F(t) from the posterior: Reported cumulative cases (red line), and overlay plots of F(t) for (a) Malaysia, (b) Selangor and (c) Sarawak.

https://doi.org/10.1371/journal.pone.0252136.g010

The effects of reporting delays on observed case trajectories and parameter inference are illustrated here based on a simulation study. A delay-in-reporting model based on the multinomial distribution is assumed: Let X ∼ Mult(D_t;p₁, p₂, ⋯, p_K), where X = (X₁, X₂, ⋯, X_K) with K = 5 and X_k is the number of cases (out of the total reported cases on day t, D_t) that is to be redistributed to day t − k + 1 for k = 1, 2, ⋯, K. The probabilities p_k, k = 1, 2, ⋯, K are chosen according to a truncated geometric distribution taking values in k − 1 for k = 1, 2, ⋯, K with success probability 0.4. The cases redistribution model is applied to new cases reported from Day 10 until Day 15. The redistributed reported case trajectory, the best fit curves and associated variabilities are shown in Fig 11. Comparing Figs 10(a) and 11(b), one can immediately notice that the reported case numbers in Fig 11(b) are better explained by the variabilities of the underlying model and the negative binomial likelihood. Parameter estimates and credible intervals for the redistributed case numbers are given in Table 5. The new infectious periods are still within the 6–8 day range and are consistent with previously reported literature. The redistribution of case numbers have also reduced the uncertainty around the MAP value of T* = 20: The credible interval for T* in Table 5 is narrower compared to that in Table 2. The MAP estimate of a₁ is now −1.693, which is about halfway between Selangor and Sarawak.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Effect of redistribution.

Panel (a) shows the redistributed daily cases and the corresponding best fit curves of F(t) (blue line) and ΔR_o(t) (green line) based on . Panel (b) shows the variability of the fit based on the ensemble set .

https://doi.org/10.1371/journal.pone.0252136.g011

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Summary results for Malaysia (with redistributed cases).

https://doi.org/10.1371/journal.pone.0252136.t005

A final point to be discussed is our preference for the negative binomial likelihood compared to the more traditional Poisson likelihood for modelling COVID-19 case numbers. Our initial investigation used the Poisson likelihood for reported case numbers but we found that the underlying model together with the Poisson likelihood was not able to capture inherent variabilities in the observed data. Hence, we opted for the overdispersed negative binomial likelihood which was able to satisfactorily represent the observed data via its overdispersion parameter τ. This is evidenced by the variability bands presented in Figs 10 and 11(b) which successfully enclose most of the reported case numbers. This coverage is further improved in Fig 11(b) by a redistribution of delayed cases. We also provide the loglikelihood values corresponding to the Poisson and negative binomial observation models in Table 6 for Malaysia (with original case numbers), Malaysia (with redistributed case numbers), Selangor and Sarawak. Note that the negative binomial loglikelihood values are consistently larger than the Poisson counterparts indicating a better model fit to observed data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. Loglikelihood values of the NB and poisson likelihoods.

https://doi.org/10.1371/journal.pone.0252136.t006

Conclusion

Quantitative models and assessment of the impacts of the Sri Petaling gathering and implementation of MCO on COVID-19 spread in Malaysia are developed in this paper. The MCO implementation is found to be highly effective in containing (an exponential rise of) the COVID-19 outbreak in Malaysia. The analysis here quantitatively demonstrates how quickly transmission rates fall under effective NPI implementation within a short time period. Higher disease transmission is found in Selangor (a state with higher population density) compared to Sarawak. We also found that under MCO, the decline in transmission was faster in Selangor compared to Sarawak. The rise and fall of disease transmission in Selangor mirrored the national level whereas Sarawak showed a more gradual increase and decrease in COVID-19 transmission. The change points were mostly found to be close to the date of MCO implementation (18th March 2020) although Sarawak exhibited a larger uncertainty around that date due to its gradual and slower increasing and decreasing trends of reported case numbers. Our study developed a new model to represent COVID-19 spread in Malaysia that accounts for heterogeneity and asymptomatic transmissions. We found that reported case numbers in Malaysia exhibited large variabilities which can possibly be attributed to a delay in reporting, particularly during the early stages of the pandemic as the experience with handling COVID-19 was new to the country. Nevertheless, the model developed here together with the overdispersed negative binomial likelihood are able to capture salient features of COVID-19 spread in Malaysia and provide reliable quantitative assessments even under the challenges of limited and delayed data.

Supporting information

Acknowledgments

The authors would like to thank the Director General of Health Malaysia for his permission to publish this article.