Right-censoring, the most commonly occurring type of censoring, happens when the event of interest is not observed, but is known to happen after a certain time-point, for instance after the end of the study period. In most longitudinal infectious disease studies individual data are routinely collected at regular time points, like visits to a clinic. At these visits it is observed whether an individual has already experienced the event or not. Then, the exact time of an event is not observed, but rather known to be occurred between two visits. This type of censoring is known as interval-censoring. If the start time is interval-censored as well, for instance if the origin event is an initial infection, the data are double-interval-censored. A well-known example of double-interval-censored data in literature is the study of AIDS induction time, i.e., the time from HIV infection to onset of AIDS. The time of HIV infection and the time of onset of AIDS are both interval-censored [6–8]. Also in the illustrating malaria data we have to deal with double-interval-censoring: the time of malaria infection and the time of initiation of gametocyte production are both interval-censored.

In survival analysis we often have to deal with a competing event, i.e., an event that prevents the event of interest from being observed [9]. This happens for many infectious diseases where the origin event (infection) can be cleared before the event of interest is observed. For instance, in the malaria data in which we are interested in studying the time from infection to gametocyte production, the infection might be cleared before gametocyte initiation. Therefore, clearance of the infection can be treated as a competing event. When competing events are present in the data, they should be modeled as such using competing risks models [9, 10]. For infections that cannot resolve, like HIV, death could be seen as a competing event.

To the best of our knowledge there is no available literature on regression models for modeling competing events with double-interval-censored data, and no available software. Adamic [11] looked at non-parametric modeling of double-censored data with competing risks but did not consider interval-censored data. Sankaran and Sreedevi [12] looked at a model for double-censored data (also not considering interval-censoring) with competing risks using an extension of the Gray [13] approach. In this paper we aim to develop and compare models for double-interval-censored competing risks data of infectious diseases. We restrict ourselves to situation in which one is interested in the time from infection to an event of interest with clearance of the infection as a competing risk. We illustrate the methodology through a simulation study and through a real-world example in malaria transmission epidemiology. Through the simulation study we also show how well these models compare to standard competing risks models when the exact times are observed (or right-censored) instead of interval-censored. Through the malaria example, we show how these models can be applied to infectious diseases data where the start time, which is not known exactly, is defined by an infection that can be cleared. We share available code for reproducing such analyses.

In Section 2, we describe the methodology for analyzing double-interval-censored data with competing events. We propose different statistical models for the analysis. Thereafter, in Section 3, we examine and compare the performance of these models under various scenarios with a simulation study. In Section 4, we look at the illustrative example using malaria transmission data. We end the paper with a discussion and conclusion in Section 5.

Let the survival time T be a non-negative continuous random variable representing the time from an origin event until the time of an event of interest. The survival function S(t) = P(T>t) is defined as the probability the event happens after time t and equals the compliment of the cumulative incidence at time t. It is either estimated non-parametrically by the Kaplan-Meier estimator or through its relationship with the hazard function λ, i.e., S(t) = exp[−Λ(t)] where Λ(t)=∫0tλ(s)ds is the cumulative hazard function. The hazard function is the instantaneous rate at which an individual experiences an event of interest over time conditional that he/she has remained event-free until that time, i.e.,

In regression modeling, where we aim to study the relationship between a covariate vector x and the survival time T, the emphasis is on the hazard function. The most commonly used regression model for time-to-event data is the Cox proportional hazards (CPH) model [14]. In this semi-parametric model the hazard function is defined as λ(t)=λ0(t)exp(β′x), where β is the vector of unknown regression parameters and λ₀ the unknown baseline hazard function. The model assumes that covariate effects β are constant over time [i.e., proportional hazards (PH)]. Flexible models that address violations in the proportional hazards assumption with applications to infectious diseases have been discussed in greater detail [15, 16]. Many other extensions of the CPH model have been proposed to accommodate complexities in the data generating mechanisms and/or the data collection limitations such as left-truncation, left-censoring, interval-censoring, competing risks, recurrent events, etc. [17–21], but in none of them the combination of double-interval-censoring and competing risks, which regularly is the case in infectious disease data, is studied.

The baseline hazard function can be specified through a parametric distribution for the survival time T. The exponential and the Weibull distributions are common (and easily implemented) specifications. An exponential distribution assumes that the baseline hazard rate is constant over time, while a Weibull distribution implies a monotonically increasing/decreasing baseline hazard [22]. The parametric forms of these distributions, while useful, may be too simple for the true baseline hazards which may lead to biases in the estimated covariate effects. A spline based model for the baseline hazards gives much more flexibility, so much so that these models are often considered non-parametric [23]. Flexible models using penalized splines are gaining momentum in advanced survival analysis models [15, 16, 24].

For the exponential model the baseline hazards is specified as a constant, i.e.,

where λ>0. Its cumulative baseline hazard equals Λ₀(t) = λt. The Weibull specification is

where λ, ρ>0. Note, when ρ = 1 it reduces to an exponential baseline hazard. The cumulative baseline hazard is Λ0(t)=(λt)ρ. The B-spline baseline hazard function is given by a linear combination of m B-spline basis functions and a set of corresponding coefficients, i.e.,

where θ_j>0. The cumulative baseline hazard can be specified using higher order B-splines basis functions as an approximation to the integrated B-splines basis functions [25] such that

where B_j are integrated B-splines basis functions. According to Wood [26] and other authors [15, 16], m = 10 is usually sufficient for good estimation.

Competing events are events that prevent the event of interest from being observed. They are easily modeled in a framework that redefines the survival function as the probability of not experiencing either the event of interest or the competing event before time t, and defining cumulative incidences for each event as the measure of interest (the probability of experiencing each event over time). Ignoring a competing event (i.e., by considering the time that competing event occurs as right-censored) results in inflated cumulative incidences for the event of interest [9, 10] as well as a missing piece of information in the transmission epidemiology. We use the multi-state approach for modeling competing risks [27]. Here we model the effects of covariates on the cause-specific hazard rates.

As shown in Figure 1, we define two hazard rates: one for the time from infection to the event of interest, λ1(t|x1)=λ01(t)exp(β1′x1), and one for the time from infection to its clearance (the competing risk), λ2(t|x2)=λ02(t)exp(β2′x2). Here λ₀₁ and λ₀₂ are the baseline hazard functions for both events. For each event k, k = 1, 2, x_k is a vector of covariates and β_k is a vector of regression coefficients to be estimated. The cumulative hazard rates for event k is

The marginal survival probability (probability of not experiencing either event before time t) is given by

and the cumulative incidence probabilities up to time t for event k is given by

As mentioned before, a survival time is interval-censored if the exact time is not observed, but known to lie within an interval. In many applications researchers ignore interval-censoring in the data by using the midpoint or the right-end point of the interval as the exact observed event time. This usually leads to biased estimates of the regression parameters and of the standard error [28–30]. In the case of double-interval-censoring these issues may be more pronounced.

In the likelihood the interval-censoring is dealt with by integrating over the interval in which the event occurred. If the time of the origin event is observed and the time from this origin event to the event of interest is known to fall in an interval (L, R) its leverage to the likelihood function is ∫LRf(s|x)ds=S(L|x)-S(R|x) where f(t|x) is the conditional density of T. Then the full likelihood for n individuals is given by the product

where i indicates the individual.

Parameters are estimated by maximizing the log-likelihood function. This is straightforward when the baseline hazard function is specified. Flexible models for dealing with interval-censoring using splines for baseline hazards specification are becoming increasingly popular [24, 31, 32].

In this paper, we consider double-interval-censored data. To illustrate different scenarios of the time-to-event data, the longitudinal follow-up data in calendar time of three individuals is presented in the upper panel of Figure 2. Participant 1 has an interval-censored infection time (infection occurred between L₀ and R₀) and is right-censored for the event of interest at time L₁—we do not observe an event time for this participant but we know that it had not happened before the time of the last follow-up visit at L₁. Participant 2 has an interval-censored infection time (infection happened between L₀ and R₀) and time of the event of interest (happened between L₁ and R₁), and participant 3 has an interval-censored infection time and event of interest time and the observation intervals overlap completely (both events occurred between L₀ = L₁ and R₀ = R₁). In general, L₀ is the date of the last visit before infection, R₀ is the date of the visit when the infection was confirmed to have happened before, L₁ is either the date of the last visit prior to the event of interest (or clearance of the infection) or the date of last follow-up when no event or clearance was observed (right-censored), and R₁ is the date of the visit at which it is being confirmed that the event have occurred before. For right censored data R₁ is not observed. The data in the upper panel of Figure 2 are in calendar time. For the analysis of the data we transform these calendar data to follow-up data by defining a starting point at which the follow-up time starts. The starting point t₀: = 0 is defined as the last observed moment infection had not yet taken place (so L₀ in the upper panel). By this transformation new intervals for the events are defined as can be seen in the middle panel of Figure 2; (t₀, t₁) for the infection time and (t₂, t₃) for the event of interest. For the third scenario where the two intervals overlap, we specifically note that t₂ = t₀: = 0 and t₃ = t₁. For right censored data t₃ = ∞. In the middle panel of Figure 2, we denote the unobserved times for infection and the event of interest by W and U. By definition W∈(t₀, t₁), U∈(t₂, t₃) and the unobserved survival time is T = U−W. The admissible survival times (T) for the three participants given the intervals (t₀, t₁) and (t₂, t₃) are shown in the lower panel of Figure 2. In some studies, a simple transformation of double-interval-censored data into single interval-censored data have been used, i.e., the time from (L₀, R₀) to (L₁, R₁) is converted to (L₁−R₀, R₁−L₀). De Gruttola and Lagakos [7] argues that this is not appropriate and is only valid when the density function for the time to the origin event is uniform in chronological time. In our likelihood formulation described in Section 2.5, we assume that the infection time is uniformly distributed within the interval (L₀, R₀) and make no assumptions about the density function for the time to the origin event prior to L₀.

A few methods that take double-interval-censoring into account in their analysis have been proposed: these can be described as non-parametric approaches [7], proportional hazards approaches [33, 34], and multiple imputation approaches [35] among others. Kim et al. [33] proposed a discrete time proportional hazards model in an extension of the Turnbull [36] estimator. In their work, Sun et al. [34] estimated the effects of covariates in a proportional hazards framework using estimating equations. However, they focused on interval-censored start times with an exactly observed or right-censored event-times. In the presence of a few interval-censored event times, they assumed the time exactly equal to the right-boundary of the interval and performed some sensitivity analysis. Sun [6] gives an overview of some of these approaches and discuss issues related to them with a focus on left truncation and also proportional hazards violations. In this paper, we look at modeling the hazard rates in a regression modeling approach for continuous time-to-event data. The baseline hazards are specified by well-known parametric distributions or more flexible shapes using B-splines. Estimation involves maximizing the full-likelihood or penalized likelihood when using splines. This is the topic of the next subsection.

For individual i we observe the data Ω_i = (t_0i = 0, t_1i, t_2i, t_3i, δ_1i, δ_2i, α_i, x₁, x₂), with t_0i, t_1i, t_2i, t_3i as defined in the previous subsection. The functions δ_1i and δ_2i indicate whether the event of interest or the competing event (clearance of the infection) has occurred: (δ_1i, δ_2i) = (1, 0) for the event of interest, (δ_1i, δ_2i) = (0, 1) for the competing event, (δ_1i, δ_2i) = (0, 0) none occurred (the data are right-censored) and (δ_1i, δ_2i) = (1, 1) is not possible by definition. For right-censored data t_3i: = ∞. We also define α_i = I[(t_0i, t_1i) = (t_2i, t_3i)] to indicate whether or not the event of interest was observed in the same interval as the infection. We assume that the observation intervals are sufficiently narrow to exclude the case in which infection and its clearance happen in the same interval, so that all infections are observed. Not making this assumption would complicate the methodology dramatically, as then latent infections could happen in any interval. In the malaria example the assumption is likely satisfied. The vectors x₁ and x₂ are observed covariates that might be associated with the time to the event of interest and to the competing event.

We assume that the times of the visits in calendar time (observation times of the individuals) are statistically independent of the actual unobserved event times U and W (and of the competing event). By this assumption, it is reasonable to assume that time to infection W within the interval (t_0i = 0, t_1i) (given t_1i) follows a uniform distribution at this interval with the density f(w)=1t1i.

In the following, we derive the likelihood function for the observed data of individual i (conditional all observation (visit times)). First, we do this for the four possible scenario's indicated by δ_1i, δ_2i, and α_i. The full likelihood for individual i can be found by combining these four terms.

Right censored, no overlap (δ_1i = 0, δ_2i = 0, α_i = 0)

Competing event, no overlap (δ_1i = 0, δ_2i = 1, α_i = 0)

Event of interest, no overlap (δ_1i = 1, δ_2i = 0, α_i = 0)

Event of interest, with overlap (δ_1i = 1, δ_2i = 0, α_i = 1)

The resulting log-likelihood is thus

where ϕ is a vector of the baseline parameters for either the exponential, Weibull, or B-spline formulations that are to be estimated, and the vector of regression coefficients corresponding to covariates sets x₁ and x₂. By independence of the observations of the individuals in the data, the full log likelihood function for all data equals the sum of the individual log likelihood functions: l(ϕ|Ω)=∑il(ϕ|Ωi) for Ω = ∪_iΩ_i.

For the exponential model the integrals in Equation (1) can be computed explicitly. For the Weibull and B-spline regression models this is not the case and approximations using quadrature methods could be used. More specific, for the Weibull model we used more precise quadrature methods, i.e., Gauss-Konrad quadrature [37] and for the B-splines model we used a simpler Simpson's approximation. The latter allowed for easier computation of the gradient and the information matrix which is needed for the algorithm to choose the penalty factor in the penalized likelihood shown in Equation (2).

For the exponential and Weibull hazards models, estimates for the parameters can be found by standard maximum likelihood estimation procedures. To make use of parallel computing we made use of the R package optimParallel [38] to optimize the log-likelihood using the limited-memory Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm [39]. This is an approximation of the BFGS that uses less computer memory and allows bound constraints. The parallel version, is much faster if the log-likelihood takes more than 0.1 s to compute.

For the B-splines baseline hazards formulation we add a penalty to the log likelihood for the spline parameters with a smoothing parameter that controls the amount of smoothing [32]. The penalized log-likelihood is given by

where S is a block diagonal matrix with blocks γ_kS_k for events k = 1, 2 and zeros elsewhere and γ = (γ₁, γ₂) a vector of penalty factors. Matrix S_k is the integrated second derivative of the B-spline basis function. There are different ways to select the smoothing/penalty factors. Due to the complexity of the likelihood and the presence of multiple smoothing parameters we used the approach adapted by Machado et al. [32]. They modeled multi-sate interval-censored data and used a two-step algorithm, following work from Marra et al. [40] for estimation as a whole. First, for a given set of smoothing factors they estimated the spline parameters and covariate effects. In the second step, they used the parameter and covariate effects estimates to estimate the smoothing parameter that minimizes the unbiased risk estimator [UBRE, [41]] which can be seen as an approximation to AIC [42]. When the difference between the parameter estimates in two consecutive iterations are less than a certain threshold (say 1 × 10⁻⁴) convergence is reached. For mathematical details of the UBRE formulation, see Machado et al. [32]. We adapted this approach for faster computation time such that in step 1 the parameters are estimated given a set of smoothing parameters using the parallel BFGS optimization described earlier. Furthermore, we estimated the information matrix required for the UBRE optimization via the empirical Fisher's information matrix.

To obtain the standard errors of the covariate effects an approximation of the Fisher's information matrix can be obtained through numerical methods. Point-wise 95% confidence intervals of the hazard and cumulative incidence estimates can be obtained by a bootstrap procedure where the analysis can be re-run on many datasets selected with replacement from the original data. The 2.5th and 97.5th percentiles can then be used as confidence interval limits.

In infectious disease applications baseline hazards are often expected to be unimodal, i.e., the hazards steadily increase to a certain point and decreases thereafter. An exponential or Weibull hazard may be to simple. In these cases a more flexible hazard should be chosen. In this simulation study, motivated by the malaria example described in Section 4, where much fewer incident infections result in the event of interest compared with the clearance event, we assume a “flatter” baseline hazard curve for the event of interest and one with a more defined peak for the clearance event. If a baseline hazard function is expected to be approximately flat, an exponential or a Weibull model may perform just as well.

With this simulation study we aim to evaluate

For every scenario we simulated N = 500 datasets and fitted three models based on each dataset: (1) an exponential PH model; (2) a Weibull PH model; and (3) a B-spline PH model with m₁ = m₂ = 10 basis functions for the event of interest and the competing/clearance event. Convergence in the algorithm is reached when the tolerance (i.e., the difference in parameter estimates between two consecutive iterations) is <1 × 10⁻⁴ or after a maximum of 30 iterations. Moreover, the data with the exact time-to-event (non-interval-censored) data are analyzed with a CPH model. The estimates found from the latter analyses are seen as the best possible results and used for comparison. We will simply refer to these models as the exponential, Weibull, spline, and CPH model.

For all scenarios the performance of each model across simulations are summarized in Table 1. The table gives the bias (β^-β), the empirical standard error (EMPSE, which is the standard deviation of the estimated β's across the different simulations), the model averaged standard errors (MODSE, the average of the SE's across simulations), and the % coverage (cov.) of the true β in the 95% confidence intervals for β. For the estimates of the bias, EMPSE, MODSE, and coverage Monte Carlo standard errors are included in Supplementary Table 1. To investigate sensitivity to the number of basis functions m_k, an additional summary was made for the performance measures (with Monte Carlo standard errors in brackets) of a B-spline model with m₁ = m₂ = 20 basis functions for scenario 1 (see Supplementary Table 2).

Table 2 shows the estimated mean square error for all estimates from the double-interval-censored model relative to the estimates from the CPH model based on the exact time-to-event data across all simulation scenarios. That is, cMSE=N-1∑j=1N(β^j,cph-β^j,model)2 for β^j,cph and β^j,model the estimates in the jth iteration.

Two scatter plots are used to summarize the association between the double-interval-censored models and the CPH model in terms of the estimated coefficients (Figure 4) and the estimated standard errors of the coefficients (Figure 5) across simulations for scenario 1. The blue line in each plot is the least squared line and the red line is the y = x line for one-to-one association between the CPH model estimates and the double-interval-censored model estimates. Figures for scenarios 2–5 are included in the Supplementary Figures 1–8. Figure 6 shows the estimated baseline hazards (Figures 6A,C) and cumulative incidence estimates (Figures 6B,D) from the double-interval-censored spline model across simulations for the event of interest (Figures 6A,B) and for the competing event/clearance (Figures 6C,D) in gray relative to the true baseline hazards and cumulative incidences (solid black lines). Figures for scenarios 2–5 are included in the Supplementary Figures 9–12. An additional plot was made for the hazards and cumulative incidence estimates for the B-spline model with m_k = 20 basis functions for each event k = 1, 2 for scenario 1 (see Supplementary Figure 13).

Malaria transmission epidemiology is used as an illustrative application for the models presented in Section 2. Malaria parasite life cycle involves two hosts: humans and mosquitos. Transmission from infected humans to mosquitoes begins when a subset of asexual parasites converts to female and male gametocytes. Through a blood meal, male and female gametocytes are ingested by anopheles mosquitos. Inside the mosquito parasites fuse, progress through several developmental stages until they multiply and develop into sporozoites. Once sporogenic development is complete, a mosquito becomes infectious and can infect more human hosts through a blood meal. Gametocyte density is a known driver of infectivity in P. falciparum mosquitos [45–47]. Their impact on mosquitoes infectivity was shown by feeding Anopheles coluzzi mosquitoes with natural P. falciparum gametocytes diluted samples [48].

This analysis used malaria incident data obtained from a cohort study conducted in Nagongera subcounty, Tororo district in Eastern Uganda between October 2017 and October 2019. In summary, 531 participants were enrolled from 80 randomly selected households and followed up every 28 days at a designated study clinic [49]. At each follow up, parasite densities were determined microscopically and by ultrasensitive varATS qPCR [50]. All patients that tested positive for microscopy were treated with a 3-day course of artemether-lumefantrine (Coartem, Novartis, Switzerland) while asymptomatic individuals were not treated. Incident malaria infections were estimated based on qPCR negativity. Infections were considered cleared by qPCR if a participant was parasite-free on three consecutive visits. For all qPCR positive samples, female (CCp4 mRNA transcripts) and male (PfMGET mRNA 229 transcripts) gametocyte densities were quantified by quantitative reverse transcriptase PCR [51]. At enrollment, whole blood was also collected to assess host genetic polymorphisms of the HBB gene [52], which is known to have a protective effect against the severe clinical consequences of infection [53, 54] and might also influence gametocytes [55]. Multiplicity of infection (MOI), the number of genetically distinct parasite strains co-infecting a single host, was assessed by apical membrane antigen-1(AMA-1) amplicon deep sequencing [56].

An incident infection was defined as a new malaria infection detected in an individual who was parasite-free by qPCR on three previous visits. Some follow-up times were significantly shorter if participants experienced clinical symptoms and received treatment. Since treatment interrupts the parasite transmission mechanism, we treated the moment of clinical symptoms prior to gametocyte production (if any) as a malaria clearance.

Here, we are interested in the evolution of the hazard rates for the initiation of gametocyte production over time since incident malaria infection, and factors that have a multiplicative effect on this hazard rate. Specifically, we assess the effects of host and parasite characteristics such as age at detected malaria infection (<5 years, 5–15 years, 16+ years), sex, parasite density (log10 transformed), sickle hemoglobin gene (HbS) (HbAA—wildtype HbS, HbSS—sickle cell anemia polymorphism, and HbAS—the human sickle trait polymorphism) and MOI (either single clone or multiclonal). Because of the routinely measured data collection criteria, the time of 104 incident malaria infections are interval-censored and thus the time to initiation of gametocyte production or malaria clearance are double-interval-censored. We also assume that an incident infection could not occur and sporadically clear within the 28 day window. Six individuals in the study had two incident infections each and some infections come from unique individuals in the same household. We do not adjust for this and assume independence between the time-to-event data between incident infections.

From an overview of the data it was clear that many incident infections cleared without first initiating gametocyte production. A competing risk model realizes that malaria clearance prior to gametocyte initiation impedes the observation of gametocyte initiation.

From the 104 incident malaria infections, 32 (31%) initiated gametocyte production, 68 (65%) had cleared their malaria without gametocyte production, and the remaining 4 (4%) were right-censored. Missing data were reported in the HbS gene information (3, 2.9%) and in the MOI data (37, 35.6%) due to sequencing failures. All three infected individuals with missing HbS data had cleared their infection without producing gametocytes. From the 37 infected individuals with sequencing failures for MOI, 34 had cleared their infection without gametocytes, 1 had initiated gametocyte production, and 2 were right-censored.

We first estimated the baseline hazards using the exponential, Weibull and spline models without any covariates and then performed univariable analysis with each of the five factors. In Figure 7, the estimated baseline hazards and cumulative incidences are shown. While we see that the shape of the hazards (Figures 7A–C) largely differ between models, the cumulative incidences (Figures 7D–F) are very much similar. From the latter estimates all three models suggest that by 100 days after incident infection approximately one-third of the infections initiate gametocyte production and the remaining two-thirds of infections clear, either naturally or through clinical intervention. Figure 7C suggests that the hazard of gametocyte initiation peaks around three weeks after incident malaria infection and declines rapidly thereafter. This is in line with other findings [57, 58], although they estimated closer to 2 weeks. A later increase in the hazards is also observed for longer event-free infections. In a descriptive analysis, Andolina et al. [59] showed that infection duration was strongly associated with gametocyte initiation. They showed that 82% (18/22) of the long infections (≥12 weeks/84 days) initiated gametocyte production earlier than 12 weeks from detection of incident malaria. From the remaining four infections, three had initiated gametocyte production at or after 12 weeks from detection of incident malaria infections.

In the simulation study in Section 3, we saw that increasing the number of basis functions improved the estimates of the peak hazards, without much impact on the estimated effects of covariates. We re-ran the malaria data analysis with 15 basis functions for each event and allowed a much larger maximum number of iterations (500) in the estimation algorithm before convergence is reached. The hazards and cumulative incidences are presented in Supplementary Figure 14. Here, we see that the peak hazard for gametocyte initiation is reached at 2 weeks and a slower but steady decline is observed afterwards until the later resurgence. The hazard estimates for malaria clearance without gametocytes remain largely unchanged.

The results of the effects of the covariates are summarized in Table 3. The estimated effects and corresponding p-values (for testing for association) across all models are fairly similar. A significant effect is seen for HbS and parasite density at detection of the initial infection on the initiation of gametocyte production. Specifically, human sickle-cell trait (HbAS) is associated with a nearly three times higher hazard of gametocyte initiation over time compared with wildtype (HbAA). A one unit increase in baseline log parasite density is associated with a nearly 20% increased hazard (p-value <0.001) of gametocyte initiation over time from the spline model. MOI showed a two times hazard of gametocyte initiation for those with more than one co-existing parasite genotype versus those with only one, and this was significant at a 10% level of significance (p-value <0.1) for the exponential and spline models.

From Table 3, we saw evidence of increased hazards for malaria clearance without gametocytes for higher baseline log parasite densities (HR = 1.14; 95% CI: 1.09–1.21, p-value < 0.001) and HbSS compared with HbAA at the 10% level of significance. Presentation of clinical symptoms was also considered a clearance event, and clinical symptoms are strongly associated with parasite density [60]. When modeling the effect of parasite density for the 88 initially asymptomatic incident infections only, we found no statistically significant effect of baseline parasite density (spline model HR = 1.02; 95% CI 0.93–1.12; p-value = 0.691, see Supplementary Table 3) on malaria clearance without gametocyte production.