Many efforts have been made by the French Regional Health Agencies in Mayotte to collect the data. The present study is based on a reported dataset of confirmed cases of COVID-19 from 13 March 2020 to 11 January 2022. The reported temporal daily cases and total confirmed cases are illustrated in Figure 1. We also consider a subset of the database containing the daily observed cases of confirmed hospitalization, intensive care, and death cases from 13 March 2020 to 26 February 2021 in Figure 1.

Some major studies have been made on the determination of generation times through serial interval distributions and incubation time from real data, refer to for instance [13, 14]. In this section, we will focus on the data collected by He et al. [13] where the authors report temporal patterns of viral shedding in 94 patients with laboratory-confirmed COVID-19 and modeled COVID-19 infectiousness profiles from a separate sample of 77 infector-infectee transmission pairs. The variation between individuals (infector/infectee) and the transmission process is summarized, respectively by the distribution of incubation period and serial intervals. Figure 2 (picture above) contains the number of infected cases (in red), the lower limit of infectors (in black), and the upper limit of infectors (in green).

The variation between individuals and chains of transmission is summarized, respectively by the distribution of incubation periods and the distribution of serial intervals (refer to definition below). For simplicity, we consider non-negative serial interval values which represent 93% of the data. Figure 2 shows the histogram (picture below) of the serial interval from the underline database.

The maximum likelihood estimation (refer to [17] is a statistical method used to estimate the parameters of the probability law of a given observation x = (x₁, .., x_n) with probability densities f(x_i, θ) by determining the values of the parameters maximizing the log-likelihood function given by:

where θ is the (vector) parameter to be estimated. The likelihood function then represents the joint density of the individual observations x_i for any given level of the distribution parameters. The maximum likelihood estimate is the distribution parameter values, which maximize the likelihood function since these same values of the estimators also maximize the log-likelihood function.

We present a mixture estimation which is a flexible tool (see [18, 19]) to model a known smooth probability density function as a weighted sum of parametric density functions f_j(x, θ_j) in a possibly multivariate independent observations x = (x₁, .., x_n) drawn from k several densities:

where k represents the number of classes or groups of the mixing, p_j>0 is the relative proportion of observations from each class j such that ∑j=1npj=1. The mixture density will allow finding the parameter Θ = (p_j, θ_j) of the model where θ_j is the (vector) parameter from class j. The log-likelihood of the mixing density is defined by:

The maximum likelihood estimates are:

subject to

For k>1, the sum of the term appearing inside a logarithm makes this optimization quite difficult. Solution of this non-linear optimization problem has long been a difficult task for applied researchers. There are, of course, many general iterative procedures that are suitable for finding an approximate solution to the likelihood equations.

Understanding the development of an epidemic is important as well as the basic reproduction number ℝ₀ at the start of an epidemic and the time varying reproduction numbers. It is historically defined as the average number of new cases of infection generated by an individual during a period of infectivity (refer to [20]). There are several methods for computing this parameter, e.g., the next-generation matrix method. Let us mention here a non-parametric approach following Wallinga and Lipsitch [21]. Let η be the probability that an individual will remain infectious over a unit of time after being infected (i.e., age of infection) and β the transmission rate at age of infection a. Then τ(a) = η(a)β(a) is the transmissibility of an infectious individual at the age of infection a, assuming that the entire population is susceptible. Hence,

We can normalize τ(a) to be a probability density function:

so that

Let Γ(t) be the number of new infections during the time interval ]t; t+dt[ and s(t) the proportion of susceptibles in the population. Note that the new infections at time t are the sum of all infections caused by infectious individuals infected at a unit of time (i.e., at time t = a) if they remain infectious at time t (with an infectious age a). In other words

For t≠0, we have a non parametric formula of the time varying reproduction number

For practical calculation, this formula offers a discrete version as follows

implemented in the package R0, Obadia et al. [22] and Boelle et al. [23] and package EpiEstim, [24, 25]. In this study, we will compare our results on the reproduction numbers with the R₀ package. For more on the calculation of the R₀, we refer to Batista [26]. Other improved estimation methods also exist in the literature, refer to Demongeot et al. [27] and Waku et al. [28], and references therein.

A regression method (refer to [29]) is a predictive statistical approach for modeling the relationship between a dependent variable say Y with a given set of observed variables say X = (X₁, X₂, …, X_n) under the form Y = f(X, η)+ω, where η is an unknown (vector) parameter and ω representing an additive error that may stand for white noise. The prediction models provide the insights that result in tangible outcomes.

The predictive statistical approach for disease research work found in the literature deals with ARIMA models [30, 31], polynomial regression and hierarchical polynomial regression models [32–34]. Recall that a process Y_t is said to be ARIMA(p, d, q) if

where Y_t−1, Y_t−2, ..., Y_t−p are the lags (past values); ϕ₁, ϕ₂, ..., ϕ_p are lag coefficients that are estimated by the model; ω_t−1, ω_t−2, ..., ω_t−q are error terms of the model for the respective lags; θ₁, θ₂, ..., θ_q their coefficients, and

with μ the mean of the process. When analyzing a model, it is important to measure its performance in order to draw the best conclusion and interpretation of the data. We will use the coefficient of determination or R² as a metric given by

where for n values marked y₁, ..., y_n, each associated with a fitted value ŷ₁, ..., ŷ_n, the sum of squares of residuals SSres=∑in(yi-ŷi), and the total sum of squares SStot=∑i=1n(yi-ȳ).

We propose a simple model that generalizes the SEIR model commonly used for virus disease modeling. The total population is partitioned into the following compartments: Susceptible S, Exposed E, Infected but not yet infectious, in a latent homogeneous period λ⁻¹; symptomatic Infected I_s with infectious period γ1-1 and asymptomatic Infected I_a in an infectious period γ2-1; Recovered R cases, and disease Dead cases D. We assume that every Exposed subject that enters I_s cohort later develops symptoms with probability 1−q and those who do not enter the I_a cohort with probability q. This assumption comes from real information in Mayotte and can be seen also in some studies, refer to Mizumoto et al. [35, 36]. We assume that a protected fraction ω of recovered individuals move to S (loss of immunity and re-infection). Given initial conditions S(0)>0, E(0)>0, I_s(0)>0, I_a(0)>0, R(0)>0, and D(0)≥0, the differential equations that govern the trajectories of our model are formulated as follow

where N = S+E+I_s+I_a+R+D is the fixed total population. With the initial conditions, the above system is well posed, with all variables remaining non-negative and N positive constant.

The epidemic equilibrium X₀ of the system is obtained by setting all the derivatives to zero with I_s = I_a = 0, which yields to X₀ = (N, 0, 0, 0, 0). To obtain the basic reproduction number, we use the next generation matrix (Dieckmann) by writing the Jacobian matrix of the flows of individuals between the different compartments at equilibrium:

The Jacobian matrix FV⁻¹ has two zero and one positive eigenvalues so that

Considering the negligible death rate, we have

Thus

where ratio r:=γ1γ2. We assume in this article that, due to the isolation the symptomatic individual reduces their infectiousness so we shall consider values of r less than one.

The basic method of moments (Figure 3A) and maximum likelihood estimation (Figure 3B) are summarized in Figure 3 for a set of classical distributions.

As we can see in the method of moments, the «Gamma »distribution seems to fit better than others. Note that the method of moments does not give a good fit since it does not use any optimization techniques such as the Maximum Likelihood Estimation (MLE) method. The fit with the exponential distribution is not good and in the sequel, we will focus on «Gamma », «Lognormal », and «Weibull »distributions. The Kolmogorov-Smirnov test result for them is given also in Table 1 for the method of moments.

The maximum likelihood estimate is done using the R software function nlm, refer to Obadia et al. [22]. In Figure 3B, we represent several fitted densities estimated by this method.

We can see that, compared to the method of moments, the maximum likelihood improves in some way the estimation, refer to Table 2 for the Kolmogorov test. Akaike's Information Criterion (AIC) and the Bayesian Information Criterion (BIC) are written as follows:

where k is the number of parameters to estimate, L is the maximum of the likelihood function of the model, and N is the number of observations in the sample. The log-likelihood values, widely applicable information AIC and BIC criteria, and estimate parameters are summarize in Table 2. We note that the log-normal distribution is the best to fit the data in this context.

In the case of a mixture, a two-component mixture model is a reasonable model based on the bi-modality. Using the R software function optim, it is possible to minimize the log-likelihood functions of the different laws taken as arguments in addition to the initialization of the parameters to be optimized for the serial interval distributions. In Figures 3, 4, we summarize the estimation results of the mixture models.

The Log-likelihood, AIC, and BIC values are given in Table 3, and we conclude that the best-fit mixture is given by the mentioned mixture of «Gamma »distributions. This best-fit mixture is better than the MLE.

In this section, we propose to modify the generation time function provided in the R0 package (refer to [22, 23]) to take into account a mixture modeling framework. Indeed, the choices of series interval estimation models in the R0 package are fixed distributions («Gamma », «Lognormal », and «Weibull ») and do not allow mixture distributions of the serial interval (SI). Furthermore, note that there is no reason for a mixture of Gamma distributions with different non-integer parameters to be Gamma distributed. Thus, using the epidemic incidence curve in Mayotte, we derive a generation time distribution in order to estimate the time varying reproduction number R₀(t) using the non parametric formula in Equation (1) and the best-fit mixture of the SI. The picture in Figure 5 represents the evolution of the time varying reproduction number R₀(t); we make the curve smooth using estimated values.

The picture in Figure 6 shows not only the evolution of the time-varying reproduction number R(t) with package R0 but also a comparison with our method of estimation in a given period. We can see that at the start of the epidemic the basic reproductive number is around 3. This is consistent with the literature [1, 6, 37] and show that, at the start of the epidemic, an infected individual could contaminate an average of 3 individuals. The public policy measures taken by the authorities and health agency of Mayotte in response to the COVID-19 with a total lockdown on March 17 made it possible to observe a value of this indicator below 1 around the middle of April 2020, showing a temporary decline in the epidemic curve during this period. After then, unfortunately, the dynamic changed over time since the end of the lockdown at the end of May 2020. Subsequently, the Mayotte authorities implemented several measures including non-pharmaceutical intervention measures such as the establishment of curfews, the wearing of mandatory masks in public spaces, and now vaccination. As we can see two other maximum values (explained by the appearance of new variants being more contagious over time) were observed at the end of 2020 and 2021 which justified restriction decisions by the local authorities. Figure 7 shows the relative variation of the reproduction numbers under other mixtures of the SI distributions.

We give some estimated results of the transmission rates for some given periods such as March 2020, December 2020 and 2021. The following pictures in Figure 8 show the range values of the time-varying transmission rates β(t) using the formula (2) according to the time-varying reproduction numbers R₀(t) and infectious rate γ₁. We make a sensitivity analysis for some values of the asymptotic ratio q. The information we have about the data allow to choose q = 0.2 or q = 0.3. For an infectious ratio r~0.7, one observed transmission rate going around up to 1.6 in March 2020 for a range from 20 % to 30 % of asymptomatic ratio. In the same conditions, the range of transmission rate goes up to 0.9 in December 2020 and up to 2.6 in December 2021. This show strong contact rates in December 2021, which can be explained by the traditional weddings in December in Mayotte. In December 2020, the contact rates were not strong and this can be explained by perhaps less contagiousness of the first variant at this time or by the fact that there were fewer traditional weddings in this same period of 2020 because of the fears of the population.

The decomposition of the time series of daily hospitalized, intensive care unit, and death cases, allowed us to choose the best models summarized in Table 4. The prediction of the above models is given in Figure 9.

We analyze in Figure 10, the performance of a polynomial regression model on the cumulative cases with a common R² = 0.998.. The model equations read as follows: