Weibull-Type Incubation Period and Time of Exposure Using γ-Divergence

TL;DR

This paper introduces a robust statistical method to estimate exposure time and incubation period during disease outbreaks, even when data includes secondary infections.

Contribution

A novel robust estimation framework using a three-parameter Weibull distribution and γ-divergence for exposure time and incubation period estimation.

Findings

The proposed method outperforms conventional approaches in bias and mean squared error in simulations.

The method shows practical advantages when applied to real-world COVID-19 surveillance data.

It is robust to data contamination from secondary or tertiary infections.

Abstract

Accurately determining the exposure time to an infectious pathogen, together with the corresponding incubation period, is vital for identifying infection sources and implementing targeted public health interventions. However, real-world outbreak data often include outliers—namely, tertiary or subsequent infection cases not directly linked to the initial source—that complicate the estimation of exposure time. To address this challenge, we introduce a robust estimation framework based on a three-parameter Weibull distribution in which the location parameter naturally corresponds to the unknown exposure time. Our method employs a γ-divergence criterion—a robust generalization of the standard cross-entropy criterion—optimized via a tailored majorization–minimization (MM) algorithm designed to guarantee a monotonic decrease in the objective function despite the non-convexity typically present in robust formulations. Extensive Monte Carlo simulations demonstrate that our approach outperforms conventional estimation methods in terms of bias and mean squared error as well as in estimating the incubation period. Moreover, applications to real-world surveillance data on COVID-19 illustrate the practical advantages of the proposed method. These findings highlight the method’s robustness and efficiency in scenarios where data contamination from secondary or tertiary infections is common, showing its potential value for early outbreak detection and rapid epidemiological response.

Figures

Full text

1. Introduction

Identifying the precise time of exposure to a newly emerging infectious pathogen using symptom onset data is a fundamental step for locating infection sources and implementing effective public health measures such as contact tracing, particularly in diseases transmissible via human-to-human contact [1,2]. Once the exposure time is well estimated, one can also determine the incubation period, which is defined as the time from the exposure event to symptom onset [3]. Historically, a variety of parametric distributions including exponential, lognormal, and Weibull distributions have been employed to model this period [4,5]. Nishiura (2007) provides a comprehensive account of their application in infectious disease epidemiology [6].

However, real-world outbreak surveillance systems rarely supply “clean” data in which all individuals share the same single-exposure event when estimating the incubation period. In practice, many datasets are inherently mixed: they contain not only secondary infections (Case 1), who were exposed at the original infection source, but also tertiary or later infections (Case 2), whose exposure stems from secondary transmissions. Data from Case 2 do not conform to the single-exposure assumption and thus act as outliers in the estimation process. Although excluding these outliers a priori would be ideal, the detailed investigations to confidently remove them are often expensive and time-consuming tasks, especially in urgent pandemic surveillance contexts. As a motivating example, Figure 1, which will be revisited in Section 4 of this paper, shows an epidemic curve of COVID-19 in Tianjin, China (21 January–12 February 2020), where four tertiary cases deviate considerably from the main cluster yet are not straightforward to exclude.

Our previous work tackled a similar situation using a three-parameter lognormal framework, proposing a robust estimation approach based on gamma-divergence, which is a robust divergence measure generalizing standard cross-entropy [8] that mitigated the impact of outliers on incubation period estimates [9]. While that approach proved valuable for distributions with long right tails, other scenarios may call for alternative parametric assumptions. In particular, the three-parameter Weibull distribution may be more suitable if the data lack extremely heavy tails or if epidemiological insight suggests that the hazard rate of onset changes monotonically over time rather than following the shape implied by the lognormal form. Nonetheless, just as with the lognormal model, conventional maximum likelihood methods applied to the Weibull distribution remain susceptible to contamination from tertiary infections. In addition, a general $[eqn]$ -divergence-based estimation framework, which can in principle be applied to various parametric models, has been proposed by Okuno (2023) [10]. Although this approach can be extended to outbreak datasets, it does not specifically target the optimization challenges posed by a three-parameter Weibull model under severe data contamination (e.g., tertiary infections). In this study, we focus on a dedicated method tailored to the three-parameter Weibull distribution and its estimation procedure. We develop a specialized majorization–minimization (MM) algorithm that guarantees a monotonic decrease in our $[eqn]$ -based objective function and also derive its covariance matrix. This specialization allows us to handle mixed or outlier-contaminated outbreak data more stably in practice while still benefiting from the robust properties of $[eqn]$ -divergence.

To address this issue, the current study adapts a $[eqn]$ -divergence-based robust methodology to the three-parameter Weibull setting, interpreting its location parameter as the unknown exposure time. We also develop a tailored MM algorithm to optimize a $[eqn]$ -cross-entropy criterion, enabling a stable estimation of the shape, scale, and location parameters—even under mixed or contaminated data conditions. Simulation experiments demonstrate that our approach substantially reduces bias and variance compared to traditional estimators, while an application to real-world surveillance data for COVID-19 further highlights its practicality.

This article is organized as follows: In Section 2, we review $[eqn]$ -divergence and introduce the associated objective function tailored to the three-parameter Weibull distribution of our interest. We then present an optimization method based on the MM algorithm. Next, a sandwich-type estimator for the covariance matrix is proposed using the theory of M-estimation. Simulation and real-world data analyses employing epidemiological surveillance data for COVID-19 and hepatitis A are described in Section 3 and Section 4, respectively. The article ends with a discussion in Section 5.

2. Method

2.1. Three-Parameter Weibull Distribution for Estimating the Exposure Time to Infectious Source and Incubation Period

Let $[eqn]$ be the disease onset timing of the ith individual ( $[eqn]$ ) and assumed to have the probability density function (PDF) given by

$$[eqn]$$

The cumulative distribution function (cdf) and hazard function are, respectively,

$$[eqn]$$

and

$$[eqn]$$

The three-parameter Weibull distribution reduces to the two-parameter (i.e., conventional) Weibull distribution when $[eqn]$ . Figure 2 illustrates the PDFs in Equation (1) when $[eqn]$ .

Here, we assume the single and simultaneous exposure to the infectious source: i.e., every individual was supposed to be exposed to one infectious source at the same time. Under this assumption of single-point exposure, the notable advantage of the three-parameter Weibull distribution is the fact that the support of Y (the disease onset timing) ranges $[eqn]$ ; therefore, $[eqn]$ can be interpreted as “the timing of exposure to the infectious source”. Once the parameters are estimated, the average and q-percentile of the incubation period can be calculated, respectively, as

$$[eqn]$$

where q is the 100q% percentile of the three-parameter Weibull distribution.

2.2. Brief Introduction of γ-Divergence

The $[eqn]$ -divergence was defined for two PDFs by Fujisawa and Eguchi (2008) [11]. Let $[eqn]$ and $[eqn]$ be the PDFs of the data-generating and parametric model distributions of x, respectively. The $[eqn]$ -divergence is defined by

$$[eqn]$$

where $[eqn]$ is the parametric PDF of interest. Note that $[eqn]$ , which is the Kullback–Leibler (KL) divergence. This divergence satisfies the following two properties: (i) $[eqn]$ and (ii) $[eqn]$ , where c is a positive constant. Lower values of $[eqn]$ approach traditional likelihood-based methods such as KL divergence, which are efficient but less robust to outliers, while higher values prioritize robustness at the cost of efficiency.

Based on the above $[eqn]$ -divergence, the empirical version of $[eqn]$ -cross entropy between $[eqn]$ and $[eqn]$ is defined as

$$[eqn]$$

where $[eqn]$ is the empirical PDF. The robust estimator ( $[eqn]$ -estimator) is then defined as

$$[eqn]$$

Now, we consider the case where the data-generating distribution is contaminated with outliers (i.e., Case 2) and given by

$$[eqn]$$

which is a mixture of the target distribution $[eqn]$ and certain contamination distribution $[eqn]$ , and $[eqn]$ denotes the proportion of outliers. The most important assumption here is

$$[eqn]$$

which assumes the practical situation where the outliers mostly lie on the tail of the target distribution. Kanamori and Fujisawa (2015) show the robust properties from a viewpoint of latent bias [12].

2.3. γ-Entropy and MM Algorithm for Optimization

Using (1), the $[eqn]$ -cross entropy function can be written as

$$[eqn]$$

Notably, the second term in Equation (8) can be written in the simple form here, while it is not possible in many cases. To obtain the minimizer, we propose the iterative algorithm of the majorization–minimization algorithm (MM algorithm) as follows.

Let us prepare the majorization function h for the cross-entropy satisfying

$$[eqn]$$

where $[eqn]$ is the parameter value of the t-th iteration step for $[eqn]$ . The MM algorithm applies the iterative procedure by

$$[eqn]$$

It is possible to show that the objective function $[eqn]$ monotonically decreases at each step because

$$[eqn]$$

Here, we propose the majorization function for Equation (8) using Jensen’s inequality as follows:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ is a term that does not depend on the parameter $[eqn]$ . The first term on the target function $[eqn]$ is a mixture of densities, which is not easy to be optimized in general, while the first term in $[eqn]$ is a weighted log-likelihood and is easy to be optimized by using the derivatives in the Appendix A.

Using the t-th iteration values of $[eqn]$ , the $[eqn]$ -th iteration values can be obtained as follows:

$$[eqn]$$

Additionally, at the $[eqn]$ -th iteration, the values of $[eqn]$ and $[eqn]$ can be obtained by finding the root of the following equations:

$$[eqn]$$ $$[eqn]$$

Remark 1 (Identification of outliers). By plugging the estimated $[eqn]$ into Equation (1) (or the associated likelihood function), we can visually inspect the presence or absence of outliers. If a sample is an outlier, it will be plotted at the tail of the distribution. By using the estimated ε, we can identify $[eqn]$ % cases as the smallest estimated values of $[eqn]$ as outliers.

Remark 2 (Selection criterion for tuning parameter $[eqn]$ ). From the result of Sugasawa and Yonekura (2021) [13], we use the following selection criterion:

$$[eqn]$$

where

$$[eqn]$$

Then, we estimate the optimal γ as $[eqn]$ . In practice, as in the simulation section, several γ values are prepared a priori, and the one with the smallest $[eqn]$ is selected. Other simple decision rules might be based on previous studies using gamma divergence [14] and density power divergence [15,16], in which the value of 0.5 was selected as a reasonable middle ground for practical applications, offering sufficient robustness without a significant loss of efficiency. In addition, other practical approaches for tuning γ, such as cross-validation or utilizing external validation datasets, are also viable [14]. Note that $[eqn]$ is known as the Hyvarinen score, which is a model selection criterion defined via the score function (from a Bayesian perspective). Further, note also that while γ-divergence is known to exhibit relative affine invariance [11], the selection criterion $[eqn]$ itself does not necessarily share this property; thus, rescaling data (e.g., doubling the observed time) may alter the optimal value of γ.

2.4. Initial Value of MM Algorithm

The proposed MM algorithm ensures the monotonic decreasing property of the objective function. However, when Equation (8) has several local minima, it converges to a local minimum rather than the global minimum. Hence, the selection of the initial value is essential. A simple approach is to start from the estimated values using the maximum likelihood method or method of moments, which is applied in the simulation section. Another and more complex approach is to run the MM algorithm with various initial values and select the best run with the smallest value of $[eqn]$ . The procedure for creating the initial values is as follows. First, a subsample is created by randomly selecting q samples from N observations. Next, the median of each subsample is used as the initial value for $[eqn]$ , the median absolute deviation as the initial value for $[eqn]$ , and $[eqn]$ as the initial value for $[eqn]$ , and then we calculate the minimum value of Equation (8). The value of $[eqn]$ should be determined beforehand based on an expert opinion or other criteria. Repeat the above process M times and select the initial values that yield the smallest value of Equation (8).

2.5. Asymptotic Properties of θ^γ

We consider the estimation of the covariance matrix of $[eqn]$ . Let us assume some regularity conditions, which are common in the M-estimator (more precisely, in the theory of normalized estimating equation [17]). The asymptotic normality of $[eqn]$ is given as $[eqn]$ , where $[eqn]$ , $[eqn]$ and $[eqn]$ . The asymptotic covariance matrix, $[eqn]$ , can then be estimated using the sandwich-type estimator:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ can be empirically estimable. The detailed derivations are in the Appendix B.

3. Monte Carlo Simulation Experiments

3.1. Simulation Setup

To assess the performance of our approach, we conducted Monte Carlo simulation experiments, varying three key parameters: the proportion of outliers, $[eqn]$ ; the value of $[eqn]$ in the distribution that generates the outliers (thus it is denoted as $[eqn]$ ); and the number of non-outlier samples, $[eqn]$ . In particular, the main body of the data—denoted by $[eqn]$ in Equation (6)—is assumed to follow a three-parameter Weibull distribution with parameters $[eqn]$ , $[eqn]$ , and $[eqn]$ . Conversely, the outliers are drawn from the same three-parameter Weibull distribution except that $[eqn]$ takes values from the set $[eqn]$ . Note that in principle, the distribution of outliers could be any form as long as it appears in the tail of the main distribution; thus, the Weibull assumption for the outliers is adopted here primarily for simplicity rather than necessity. Moreover, since the same pathogen is implicated in both secondary (Case 1) and tertiary or subsequent (Case 2) infections, applying the same distribution is justified. Overall, considering all combinations of $[eqn]$ , and $[eqn]$ yields a total of 27 scenarios (see Table 1). For each scenario, 1000 Monte Carlo simulations were performed.

For each scenario, the procedure to generate an individual dataset for the kth scenario, $[eqn]$ , is as follows. First, we randomly generate N samples from a three-parameter Weibull distribution with parameters $[eqn]$ . Next, we randomly generate $[eqn]$ samples from a three-parameter Weibull distribution with parameters $[eqn]$ corresponding to the kth scenario.

We evaluated the performance in terms of the bias and mean squared error (MSE) of estimated mean and 95% percentile values of the true distribution, which were calculated from Equation (4) and $[eqn]$ . The comparison methods included estimates based on (1) ml: maximum likelihood method [18,19], (2) mm: method of moments [20], and (3) mps: the method of maximum product spacing [21,22]. These conventional methods are easily implemented in the R program [23] using the fitWeibull() function in the ForestFit package [24].

3.2. Simulation Results

The Monte Carlo simulation results show that the proposed method consistently outperforms the conventional approaches across most scenarios, providing less biased and more efficient estimates for both mean and 95% percentile of the true distribution (defined in Equation (6)) and $[eqn]$ . As shown in Table 2 and Table 3, our method achieves the smallest bias and MSE on average. Specifically, the bias in the mean under our approach achieved 65% reduction compared with that of the conventional approaches; the bias for our approach ranges from $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) compared to ranges of $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for ml, $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for moment, and $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for mps. The difference in performance between our approach and conventional approaches is even more pronounced in the estimation of the 95% percentile, achieving a 68% reduction compared with that of the conventional approaches; the bias for our approach ranges from $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) compared to ranges of $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for ml, $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for moment, and $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for mps. Regarding the bias of $[eqn]$ , our approach also outperforms the conventional approaches, achieving nearly a 90% reduction (excluding the results of mps, which show extremely large biases): the bias for our approach ranges from $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) compared to ranges of $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for ml, $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for moment, and $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for mps.

In terms of MSE, our method also provides markedly smaller values in the estimation of the mean (a 70% reduction overall), ranging from $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) compared to ranges of $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for ml, $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for moment, and $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for mps. Again, the performance gap is even enlarged for the 95% percentile, achieving a 76% reduction in MSE compared with the conventional approaches; MSE of our approach range from $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) compared to ranges of $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for ml, $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for moment, and $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for mps. Regarding the MSE of $[eqn]$ , our approach also outperforms the conventional approaches, achieving an 88% reduction: the MSE for our approach ranges from $[eqn]$ to $[eqn]$ (mean $[eqn]$ ), compared to ranges of $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for ml, $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for moment, and $[eqn]$ to $[eqn]$ (mean $[eqn]$ ) for mps.

Overall, conventional methods tend to suffer in performance under small sample sizes, when the proportion of outliers is high, or when outliers are not concentrated in the tail of the target distribution (i.e., small $[eqn]$ ). In contrast, our method remains robust and yields stable estimates even under these challenging conditions. In addition, although all methods show improved performance (i.e., reduced bias and MSE) as N increases, our method continues to provide similar or superior performance compared to the conventional methods even at $[eqn]$ and smaller outlier proportions $[eqn]$ or $[eqn]$ (Scenarios 3 and 6).

4. Application for Real-World Data: Epidemiological Surveys for COVID-19

This section applies both the proposed method and the comparison methods to contact tracing surveillance for COVID-19, aiming to identify infection sources and estimate the incubation period. Contact tracing surveillance refers to investigations primarily conducted by local public health authorities to prevent disease spread in communities where infections have occurred. The data and corresponding R code (https://www.r-project.org/, accessed on 17 February 2025) are available on the corresponding author’s GitHub page (https://github.com/, accessed on 17 February 2025).

We focus on a COVID-19 outbreak in Tianjin Province, mainland China, from 21 January to 12 February 2020, as detailed by Wang and Teunis (2020) [7]. The dataset, consisting of 112 confirmed cases, highlights a mixture of secondary, tertiary, and subsequent infection cases with the transmission network inferred from symptom onset dates. Figure 1 displays the epidemic curve: the first exposure is designated as day 0, and the reported onset dates span day 1 to day 24. In the figure, 31 secondary infection cases are shown in gray, whereas 39 tertiary or subsequent infection cases are shown in black.

To estimate the exact exposure time, the tertiary and subsequent cases would ideally be excluded, but doing so requires extensive investigation, which is typically time consuming and expensive in the context of emerging infectious diseases. Consequently, using a “mixed” dataset without removing these cases is common in practice as is treating the tertiary and subsequent cases as outliers in the analysis.

Table 4 briefly summarizes a comparison of the estimated time of exposure to the infectious sources, $[eqn]$ , across our approach and the conventional methods. With the selected $[eqn]$ of 1, our method produces the estimates of $[eqn]$ with the corresponding 95% CI for $[eqn]$ , mean and 95% percentile of the distribution of incubation period: $[eqn]$ , and the mean and 95% percentile are 3.96 and 9.77. We note that the estimated $[eqn]$ is quite close to the actual time of exposure (day 0). In contrast, the conventional approaches produce MLE values of $[eqn]$ ranging from −0.49 to 2.20. Clearly, our method succeeds in returning the preferable estimated value of $[eqn]$ closer to the realistic exposure time. In terms of the distribution of the incubation period, our method provides a mean and 95% percentile of 3.96 and 9.77, respectively, whereas the conventional methods provide estimates ranging from 3.89 to 4.28 for the mean and from 6.96 to 10.63 for the 95% percentile.

Our robust estimation method has practical implications for understanding the biological evolution of COVID-19, particularly regarding changes in incubation period distributions associated with different SARS-CoV-2 variants. Accurately estimating exposure times and incubation periods in contaminated datasets enables epidemiologists to better capture subtle shifts in viral characteristics, such as transmissibility and generation intervals, across successive infection generations. Such insights can be crucial when assessing how viral mutations or emerging variants alter epidemiological parameters, influencing the trajectory of outbreaks and informing timely public health responses.

5. Discussion

We have introduced a novel robust approach for estimating both the exposure time to infectious sources and the incubation period based on the $[eqn]$ -divergence approach for the Weibull distribution. This approach maintains robustness even under substantial contamination, which often arises when unexpected secondary or tertiary cases are captured in rapid epidemiological surveillance. A frequent challenge in robust estimation lies in developing an efficient algorithm, especially given the non-convex and non-differentiable nature of many robust objective functions. In this study, we devised a practical estimation method that guarantees a monotonic decrease in the objective function by leveraging the MM algorithm.

Although our analysis assumed a contaminated density of the form $[eqn]$ , we note that this setup can be generalized to a more intricate mixture, $[eqn]$ , using essentially the same framework. Numerical simulations and applications to real-world data consistently indicate that our method surpasses conventional estimators in terms of bias, MSE, and both the mean and 95% percentile of the true distribution.

A key remaining question in the realm of infectious disease surveillance involves the practical selection of the incubation period distribution. Currently, we assume that the incubation period follows the Weibull distribution. However, the accuracy of the estimation strongly depends on the validity of this assumed distribution. If the true distribution deviates from the assumed one, the estimates of exposure time may be biased or imprecise, resulting in an inaccurate estimation of the incubation period. Our robust approach can be extended to other types of distributions. Although the idea is straightforward, further efforts are required to derive a new MM algorithm and covariance matrices. This will be the subject of our future research.