Estimation of trajectory of protective efficacy in infectious disease prevention trials using recurrent event times

In studies of infectious disease prevention, the level of protective efficacy of medicinal products such as vaccines and prophylactic drugs tends to vary over time. Many products require administration of multiple doses at scheduled times, as opposed to one‐off or continual intervention. Accurate information on the trajectory of the level of protective efficacy over time facilitates informed clinical recommendations and implementation strategies, for example, with respect to the timing of administration of the doses. Based on concepts from pharmacokinetic and pharmacodynamic modeling, we propose a non‐linear function for modeling the trajectory after each dose. The cumulative effect of multiple doses of the products is captured by an additive series of the function. The model has the advantages of parsimony and interpretability, while remaining flexible in capturing features of the trajectories. We incorporate this series into the Andersen‐Gill model for analysis of recurrent event time data and compare it with alternative parametric and non‐parametric functions. We use data on clinical malaria disease episodes from a trial of four doses of an anti‐malarial drug combination for chemoprevention to illustrate, and evaluate the performance of the methods using simulation. The proposed method out‐performed the alternatives in the analysis of real data in terms of Akaike and Bayesian Information Criterion. It also accurately captured the features of the protective efficacy trajectory such as the area under curve in simulations. The proposed method has strong potential to enhance the evaluation of disease prevention measures and improve their implementation strategies.


INTRODUCTION
In studies of infectious disease prevention, protective efficacy (PE) of medicinal products such as vaccines and prophylactic drugs is often defined as one minus hazard ratio (HR), where HR < 1 indicates success of the products in reducing disease incidence.The level of PE typical rises initially to a peak and subsequently decays due to various reasons such as waning immunological responses, changes in naturally acquired immunity, and elimination of drugs from the body.It is important to characterise this trajectory, for understanding the potential impact of the intervention and developing an appropriate strategy to deploy it.2][3][4][5] Disease incidence and person-time between the date of dose one and this start point may be included in intention-to-treat analysis, but the efficacy during this initial period, which can be weeks or months, is not usually reported separately.The Covid-19 pandemic has led to rethink about this practice.The global shortage of Covid-19 vaccines drew attention to the level and tempo of the initial rise of PE after dose one.It influenced the decision to delay the second dose of the primary series or not.If the initial rise of PE after dose one was sufficient, the second dose could be delayed and the limited vaccine supply could be prioritized for administration as the first dose to protect more people. 6,7This has highlighted the relevance of the early part of the PE trajectory.
The waning of PE after it reaches a peak level has long been a concern in the evaluation of preventive medicine and disease prevention strategies.][10][11] Accurate information on the rapidity of PE waning can guide the frequency and timing of booster doses or treatment re-administration.It can also guide the timing and method to deploy the first dose.For example, if PE is stable over a year, a vaccine may be suitable for delivery in immunization clinics according to some age schedule.In contrast, if PE wanes quickly, it is more suitable for delivery through a community campaign timed before the high disease transmission season.
Some studies have attempted to estimate "duration of protection" offered by preventive medical interventions.However, the search for a duration of protection may be inappropriate, as it assumes a sudden drop of PE from a protective to a non-protective level.It is more realistic and informative to estimate the trajectory of PE level over time.2][13] It can work well if the number of events is very large.Otherwise, it may encounter volatility in PE estimates if the intervals are narrow, or loss of information if the intervals are wide.Some researchers have used exponential decay function or its adaptation to capture time-varying effects. 14,15Such patterns are not suitable for studies of the initial rise and possibly a plateau of PE that precede the decay.Some researchers smoothed the densities of time from initiation of intervention to outcome event among those who have had the event time observed, and compared the smoothed densities between trial arms. 10This may generate biased results for various reasons, for example, different patterns of censoring between trial arms.
Infected individuals may acquire natural immunity, which may be short-lived or long-lasting, and may be partial or near complete.Prevention of infection or disease episodes thus may coincidentally lead to a lower rate of acquisition of natural immunity, which in the long term increases the disease risk.These effects may influence the trajectory of PE.This problem is known as rebound effects in infectious disease epidemiology. 16,17The concern for rebound effects delayed the World Health Organization's recommendation for and widespread deployment of seasonal malaria chemoprevention, which has been a major public health success since being recommended in 2012.But it could have been deployed 20 years earlier. 18,19ifetime immunity after one episode of an infectious disease is uncommon.Many infectious diseases can recur, and recurrences can be frequent within the duration of clinical trials.The World Health Organization Malaria Vaccine Advisory Committee highlighted the limitations of analysis of only the first or single event and called for methodological research on analysis of recurrent event data. 20The Andersen-Gill (AG) model is an extension of the Cox model for analysis of recurrent events.2][23][24] In the analysis of the first or single event, it is known that Cox model is biased towards the null value as time goes by due to unobserved heterogeneity. 25,26This bias distorts the trajectory.The analysis of recurrent events by the AG model mitigates against this distortion. 21,22he elementary form of the AG model and Cox model assumes time-constant intervention effects, or proportional hazard (PH) assumption.A previous study proposed a 4-parameter function that represents a monotonic, non-linear pattern of time-varying effect for use with time-to-event analysis model 27 : This function has the advantage that the four parameters, A, B, C and D, are interpretable, representing the level of short-term effect, rate of waning, shape of trajectory, and level of long-term effect, respectively.A negative D value suggests rebound effect.A disadvantage of this function is that it is monotonic.It assumes intervention effects to be at or near peak level when analysis time begins.This assumption makes the function more suitable for fast-acting products, such as monoclonal antibodies, or typical vaccine trial analysis that sets 14 or 28 days after completion of the primary series as the start point of analysis time.In contrast, a non-monotonic function that captures the trajectory from zero to peak to waning and rebound is more useful for the evaluation of slower-acting products and analysis starting from the time of dose one.
A study of the time-course of the impact of a dose of antibiotic on continuous outcome variables proposed a 4-parameter, non-monotonic function of intervention effect that was derived from a harmonic oscillator. 28This was also applied to time-to-event analysis. 29This function has the advantage of capturing the whole trajectory of intervention effects from dose one.A disadvantage is that since it originates from physics/engineering, the parameters have no biomedical interpretation.
The aim of this study is to propose and evaluate a non-monotonic, non-linear function of time-varying effects that is motivated by pharmacokinetic/pharmacodynamic (PK/PD) modelling for use with the AG model.To be clear, the proposed method involves PK/PD concepts, but it does not require PK/PD biological samples.It has the advantage of parsimony and interpretability yet being flexible in fitting data.In Section 2, we will describe the statistical model.In Section 3, we will apply the model to data from a malaria prevention trial and compare it with various alternatives.Section 4 will provide an evaluation of the method by simulation.In Section 5, we will provide some concluding remarks.

Incorporating time-varying effects and repeated dosing into AG model
Let N be the total number of subjects in a 2-arm trial, z i = 1 if the i-th subject is in the intervention group and z i = 0 otherwise, and 0 < t i1 < t i2 < … < t in i be the event times and n i be the number of events of the i−th subject that are observed during the study period.Incorporating time-varying effect into the AG model, the hazard of the outcome event for subject i at time t is: where  0 (t) ≥ 0 is an unspecified baseline hazard function that is common to all events, x i (t) is a vector of possibly time-varying covariates,  is a vector of regression coefficients, and G(t) is a function of the time-varying intervention effect that is the key interest.For a subject who have received m i doses of a medicinal product, let 0 ≤ d i1 < d i2 < … < d im i be the times at dosing.Deviation from scheduled dosing times is common because of non-compliance and logistic reasons.Such deviation can dilute the pattern of time-varying effect.To allow variation in time at dosing and prevent the dilution, G(t) is formulated as: where t is the time since the first dose, g is the time-varying effect of the j−th dose the i−th subject has received, and I(⋅) is an indicator function.Since g and G i (t) refer to intervention effect on ln(HR), both the individual and cumulative dose effects have a lower limit of HR = 0 and therefore the protective efficacy has an upper limit of PE = 1.The effects of multiple doses are additive on the log scale.Thus, their effects in terms of HR are multiplicative.This is standard formulation of Cox-type models, including the AG model.

Modeling time-varying effects after each dose
Pharmacodynamics (PD) is the study of how the body responds to medicinal products, with emphasis on dose-response relationship.Our proposed method begins with a sigmoid Emax PD model. 30Let PE at time t be: where Emax is the maximal response, C 50 is the concentration of a stimulus required to achieve half of Emax,  is a sigmodicity factor that represents the shape of the curve, and C(t) is the concentration of the stimulus over time.Since PE = 1 − HR and minimal HR is zero, it follows that Emax = 1 and Define g(t) as the ln[HR(t)], The relation between g(t) and time is via C(t), which is to be described by a pharmacokinetic (PK) model.Pharmacokinetics is the study of how the body affects the medicinal products administered, with emphasis on changes over time.As Upton and Mould noted, 31 "There is no 'correct' method for developing PD models.The nature of model is dependent on the data and the intended purpose of the model."Here, the purpose of developing the function g(t) is not to estimate pharmacokinetic parameters.Instead, the purpose is to develop a flexible and parsimonious tool for the estimation of the trajectory of time-varying intervention effects.We consider using a first-order absorption PK model, which is suitable for studies of products used in extra-vascular administration like intra-muscular injection of vaccines and oral drugs, for substitution into g(t).Some parameters in the PK model may be omitted as they are not of interest to us.In particular, the ratio of dosage to volume of distribution (D/V) is a scalar that does not affect the shape of C(t).So, it can be fixed at 1 without loss of information in the present context.Consequently, where  a is the absorption rate and  is the elimination rate in PK studies.Furthermore, the shape of C(t) mainly depends on the relative size of  a and .So, we make a simplification by constraining  = 1, leaving only one parameter  a (≠ 1) in the C(t) component of the model to be estimated.
To capture the long-term effect (rebound effect), similar to the term D in Equation ( 1), we can add a term  to form: Shaw et al. proposed to allow a long-term effect to grow asymptotically, with the growth rate fixed to equal another parameter of their model for simplicity. 28Along the same line, we replace  in g(t) by: such that the faster the short-term effect rises (larger  a ), the faster the long-term effect grows, or vice versa.
Finally, we define g 1 (t) as our proposed 4-parameter PK/PD motivated function for modelling PE over time: In this function, a smaller C 50 indicates a higher level of PE, a smaller  a indicates a slower rise of C(t), and therefore, g 1 (t), a smaller γ indicates a less sigmoid shape of the trajectory.A larger  with a negative sign indicates a stronger protective effect in the longer term, while a larger  with a positive sign indicates a stronger rebound effect.These features of the function are illustrated in Figure 1.
The function is then plugged into G(t) in Equation ( 3), denoted as G 1 (t), which in turn is plugged into the AG model in Equation (2).
The first derivative of g 1 (t) with respect to t is: If  is close to zero and is ignored, a closed form expression for the time to peak level of PE ( t p ) after a single dose is available: Then, peak level of PE is: Without assuming  = 0, there is no closed form expression for t p and PE ( t p ) .In the simulation studies and real data application that will be described in the next two Sections, we obtained the solutions numerically without assuming  = 0.
The time to PE dropping to a certain fraction of the peak PE is important for guiding the timing of booster or re-administration. 27To find the time to a fraction f of peak PE remaining We refer to t 0.5 as t p∕2 .The area under the PE(t) curve (AUC) after a single dose up till some defined time-point can be found by numerical integration.In the presence of a rebound effect ( > 0), the AUC is the sum of the areas above and below zero.As shown in Figure 1, PE (y-axis) is up to one.If the metrics for the time-scale (x-axis) is months and the maximum follow-up time is 8 months as in Figure 1, the upper bound of AUC is 1 × 8 = 8, which refers to complete elimination of the disease incidence over the 8-month duration.As another illustration based on this figure, the AUC for the trajectory indicated by the solid line was 0.97.This value indicates that the amount of disease incidence reduction is equivalent to that of removing the exposure time by 0.97 month, or 0.97/8 = 12.1% of the disease incidence in the 8-month duration, assuming absence of seasonality.Alternatively, it is possible to use the study duration as one unit of time for the metrics of the x-axis, such that a duration of 8 months equals 1.0 in this figure.Then maximum AUC would be 1.0.The advantage of using months or years for the x-axis is that it facilitates comparisons between estimates from studies with different durations.A negative AUC (> −∞) is possible, which indicates harm of the product.For brevity, in this article we mainly discuss AUC after one dose, calculated from g 1 (t).But it is straight-forward to extend to AUC of a multi-dose regime, calculated from G 1 (t).
For comparison purpose, we define the 4-parameter monotonic function in Equation (1) as g 2 (t) and the 4-parameter non-monotonic function derived from harmonic oscillator by Shaw et al. 28 as g 3 (t).Furthermore, we consider step function and cubic B-splines for comparison.We refer to the step and cubic B-spline functions as g 4 (t) and g 5 (t), respectively.Details of the step and spline functions will be discussed in Section 3.
Previous observations revealed that the last dose of the primary series of Haemophilus influenzae type b vaccine offered little additional protection. 32Accordingly, we consider an option to allow C 50 in g 1 (t) to be dose-specific.That is, changing C 50 to C 50,1 , C 50,2 , • • • , C 50,m for the m doses.Another aspect is that  may change sign over time as it is a balance between the intervention group being protected by the medicinal product and the control group acquiring natural immunity more than the intervention group.This balance may change over time.So, we also consider letting  in g 1 (t) to be dose-specific.
Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) can be used to compare models with different g h (t) or models with common or dose-specific parameters.

Estimation
Let  h (h = 1, 2, … , 5) be the unknown parameters in the AG model to be estimated from the data, where h represents the choice of the function g h (t).The log-partial likelihood is: ) to be positive,  is replaced by exp( * ), where  * = ln(), and so on.The estimators θh can be obtained by solving the equations U( h ) = l( h )∕ h = 0.Under regularity conditions similar to VII.2.1 and VII.2.2 in Andersen et al., 33 it can be shown that as sample size increases towards infinity, the estimator θh of  h is consistent and normally distributed.We used the quasi-Newton method (BFGS) available in R for the estimation. 34he square root of the diagonal elements in the inverse of the information matrix, I ( θh ) , are the standard errors calculated assuming independent observations.The AG model uses a robust sandwich estimator for clustered data to allow for multiple observations per person 35,36 : where We used the robust sandwich estimator for clustered data to estimate R codes for the model estimation and the simulation in Section 3 are available at https://github.com/cheungyb/PE-trajectory.

3
MALARIA PREVENTION TRIAL

Materials and methods
We re-analyzed data from a randomized controlled trial of malaria prevention in Ghana for illustration. 27Infants were recruited from immunization clinics at the age of about 2 months.They were enrolled and randomized to receive four doses of placebo or sulfadoxine-pyrimethamine (SP) for malaria prevention, at 1, 2, 7 and 10 months after enrolment.The dose-schedule was designed to align with the immunization schedule such that the trial did not necessitate additional visits to the clinics.The actual timing of the administration of placebo/SP was variable, with mean (SD) of the first to fourth doses at 1.0 (0.3), 2.1 (0.5), 7.6 (0.9) and 10.8 (1.0) months after enrolment, respectively.The infants were under surveillance for malaria disease episodes for up to 24 months.Clinical malaria was the primary endpoint, which was defined as a visit to a health care facility with (a) microscopy confirmation of malaria parasites in the blood and (b) either temperature ≥ 37.5 • C or a parental report of fever.To avoid double-counting of disease episodes, observations of clinical malaria within 7 days of a previous clinical malaria episode of the same person were ignored. 27There were 1442 clinical malaria episodes among 1044 infants in the SP trial arm and 1593 episodes among 1001 infants in the placebo arm.
In the data analysis, date of receiving the first dose of placebo/SP for each infant is the origin of analysis time (t = 0).In all models, age at enrolment and rainy season (July to November; dry season otherwise) were included as time-constant and time-varying covariates, respectively.
We began with fitting the AG model with SP trial arm versus placebo trial arm as a time-constant exposure variable, assuming proportional hazard (PH).Then we fitted AG models with different G h (t) to capture the time-varying effects of the four doses of SP.
To challenge the performance of the parametric functions, we set g 4 (t) as a step function with eight steps (weeks 1, 2, 3, 4, 5-6, 7-8, 9-12 and ≥13).Furthermore, we fitted a series of models with time-varying effect represented by cubic B-splines, 37 and chose the best fitting one among them for g 5 (t) according to AIC and BIC.A cubic B-spline with one (inner) knot has four coefficients to be estimated, which is the same as the parametric functions we consider.We fitted cubic B-spline models with one to four knots.Previous studies of SP have suggested that its efficacy tended to last no more than 4-6 weeks. 8,38Placing knots beyond 6 weeks is unnecessary.We fitted five models with one knot, at 7, 14, 21, 28 or 35 days since dosing; four models with two knots, at (5, 20), (7, 21), (7, 28), or (14, 28) days; two models with three knots, at (5, 15, 30) or (7, 21, 35) days; one model with four knots, at (5, 15, 25, 40).The fitting of multiple spline models increased the risk of over-fitting.It was used only for the purpose of challenging the performance of the proposed PK/PD motivated model.

Results
The PH model gave AIC 44613.0 and BIC 44631.1 (Table 1).The SP trial arm was estimated to have an HR of 0.87 (PE = 0.13) as compared to the placebo trial arm.The 4-parameter PK/PD motivated model, G 1 (t), gave much smaller AIC and BIC than the PH model (Table 1), indicating that the PH model was not tenable.Both AIC and BIC demonstrated the superiority of G 1 (t) over G 2 (t) to G 5 (t).AIC and BIC identified models with cubic-B spline functions with 3 and 2 knots, respectively, as the best among all the models using spline functions (details in Online Appendix S1).Both of them had larger AIC and BIC than the model with G 1 (t).
Allowing dose-specific C 50 in G 1 (t), that is, using three more parameters, led to larger AIC and BIC than the model with a common C 50 .Similarly, allowing dose-specific  did not improve the AIC and BIC either.We thus focused on the 4-parameter model.Table 2 shows the coefficient estimates and robust standard errors for the AG model with G 1 (t).Based on the estimates, we plotted the estimated time-varying PE in Figure 2. The figure was generated assuming the four doses were received at 0, 1, 6 and 9 months (from dose 1) according to schedule.According to the PK/PD function, t p , PE ( t p ) , t p∕2 and PE ( t p∕2 ) were 8.3 days, 0.93, 34.1 days and 0.465, respectively.While the product was highly protective in the short-term, the protection did not last long.The step function is robust, though not parsimonious, as it does not involve assumption of the shape of PE trajectory.We compared the trajectory features of the short-term effect estimated from the proposed method versus this robust alternative (weekly intervals: 1, 2, 3, 4, 5-6, 7-8, 9-12).The step function indicated the second week post-dosing as t p , where PE ( t p ) = 0.91.Furthermore, PE at the 5-6 weeks (29-42 days) interval was 0.472.All these features agreed well with the model-based estimated trajectory.
There was some sign of a rebound effect,  = 0.024, but it was far from statistically significant (SE 0.017).Graphically it was represented by the PE curve dropping below zero (Figure 2).The AUC of the 4-dose regime over the 24-month trial duration was 2.25 (or 2.25/24 = 9.4% reduction).This was the sum of areas between PE(t) and zero where PE(t) > 0 minus the sum of areas between PE(t) and zero where PE(t) < 0.
The time interval between dose 2 and 3 was wide, and the short-term effect apparently waned quickly.To evaluate the PK/PD-based estimate of long-term/rebound effect, we extracted the data in the time intervals from 3 months after dose 2 till dose 3 (interval 1) and from 3 months after dose 4 to end of study (interval 2), analysed them using a proportional hazard model, controlling for the same covariates, and compared their PE estimates for these intervals against those from the PK/PD function.The PE's (95% CI) for intervals 1 and 2 estimated from this subset of data using proportional hazard model were − 0.069 (−0.280, 0.107) and −0.106 (−0.243, 0.016), respectively.Since the decaying component of the PE(t) driven by C(t) in the PK/PD function had declined to about 0.1% by 3 months post-dosing, the model-based estimates for these two time intervals using all data were approximately 1 − exp(2×) = 1 − exp(2 × 0.024) = −0.049and 1 − exp(4 × 0.024) = −0.101,respectively.The direction of change in rebound effect over cumulative doses and their magnitude in the two intervals of time were quite accurately estimated by the model.
The analysis showed that the protection offered by the product reached peak level quickly but was short-lived.The reduction in malaria incidence was small not because the product was unprotective, but because the frequency of administration was too low.To achieve a protection level of PE ≥50% at any point in time while children are at high risk of malaria, the product should be administered every 4 to 5 weeks.Furthermore, the rebound effect, if any, was mild, and was clearly out-weighed by the benefit.Avoidance of the product due to concern of rebound effect was not justified.

Simulation methods
We considered three sets of simulation scenarios: a PE trajectory following the 4-parameter PK/PD function, with two set of parameters (C 50 ,  a , γ, ) = (0.4,3, 3, 0.1) or (0.4, 9, 3, 0.1).b PE trajectory following the 5-parameter PK/PD function that has no constraint on  = 1, with two set of parameters (, C 50 ,  a , γ, ) = (2, 0.4, 3, 3, 0.1) or (0.5, 0.4, 9, 3, 0.1).c PE trajectory following piecewise linear spline, with t p = 10 or 14 days: For each PE pattern, we simulated two sample sizes, 500 or 700 per trial arm (N = 1000 or 1400 in total, respectively), and 1 or 2 doses of intervention.The gap between doses independently followed a Uniform (1, 2) months distribution in the two-dose schedule.The study duration was 5 or 7 months for one-or two-dose scenarios.For scenarios in (a), we additionally simulated a three-dose schedule, with a study duration of 9 months.Number of events in the placebo group was approximately 500, 700 and 900 events for N = 1000 and 700, 1000 and 1250 events for N = 1400 for the one, two and three dose scenarios, respectively.Censoring was generated by having 80% of the persons completing the follow-up time and the remaining 20% were uniformly censored between 0.8 and 1.0 of the planned follow-up time.We used 500 replicates per simulation scenario.Assuming the actual coverage of the estimated 95% confidence interval is indeed approximately 95%, 500 replicates offered a margin of error approximately no more than ±2% (i.e., 1.96 × √ 0.05 × 0.95∕500) in the estimation of coverage probability.We used the thinning approach for generation of the recurrent event times (details in Online Appendix S2).
We fitted the AG model with the 4-parameter PK/PD function for analysis of the data generated in all the scenarios.In scenarios under (a), we evaluated the relative bias (%) of the estimators, ratio of mean SE to empirical SD of estimates (ASE/ESD), coverage probability (CP) of the 95% confidence interval, and root mean square error (RMSE).For all scenarios under (a), (b) and (c), we derived the peak PE level, t p , t p∕2 and AUC up to 3 months after a dose, and present the average estimates of these four features as compared to their respective true values determined by the true model parameter values.

Simulation results
Table 3 shows the simulation results for scenarios that the true PE was generated by the 4-parameter function g 1 (t).In all scenarios simulated, the estimates of log(C 50 ) were practically unbiased.The other three parameters usually had less than 10% bias, and the bias reduced as sample size or number of doses increased.For all parameters, RMSE consistently reduced as sample size or number of doses increased.In all scenarios, coverage probability for log(C 50 ) and  were close to the target 95% level.In scenarios with N = 1000 and only one dose, ln( a ) and ln(γ) clearly had CP < 95%.However, the ASE/ESD and CP improved as sample size or number of doses increased.Table 4 shows the mean of the estimates of the derived parameters about the four features of PE trajectory.For scenarios that were generated by the 4-parameter PK/PD functions, the estimates closely represented the four features regardless of sample size and number of doses.For scenarios not generated by the 4-parameter PK/PD functions, the estimates closely followed the true peak PE level, t p∕2 and AUC.There was more discrepancy between the t p estimates and the true t p , with absolute bias not more than 3 days in most scenarios.The largest discrepancy was −3.5 days (8.4 vs 11.9) in one of the scenarios under (, C 50 ,  a , γ, ) = (2, 0.4, 3, 3, 0.1).

DISCUSSION
Information on the trajectory of protective efficacy over time is important in the evaluation and deployment of preventive interventions against infectious diseases.Experiences from Covid-19 and malaria have highlighted the needs to understand the pattern from the initial rise of PE to rebound effects in the long term.As compared to analysis of single event, the use of the AG model for recurrent events mitigates against the time-varying bias towards the null value due to unobserved heterogeneity.It also offers more accuracy and precision to the modeling of waning efficacy and rebound effects as it has a larger number of events in the later part of the study period.Protective efficacy of medicinal products for prevention of infectious diseases usually wanes over time.Estimation and presentation of PE as a single value as if PE is constant over time is often insufficient and sometimes misleading. 39he shorter the trial duration, the larger the single PE value may appear to be.Yet, PE summarized as a single percentage is often compared between medicinal products and between trials, without taking into account differences in duration of follow-up or in the timing of exposure to infection, for example, malaria vaccines administered uniformly over calendar year or shortly before the malaria transmission season. 40Estimation of the PE trajectory is more informative.Graphical presentation of the results can help communicate clearly.In addition, based on the trajectory information, investigators can derive the AUC as a summary measure.This summary measure is more portable across trials and can facilitate comparison even if the trials have different follow-up durations, provided that each trial has a follow-up duration sufficient for the estimation.
We have embedded the PK/PD and other trajectories in the AG model.The AG model assumes that events are indistinguishable and ordered, such that a person cannot have two events at the same time.Same as other standard time-to-event analysis methods, it also assumes non-informative censoring.It is sometimes said that the AG model assumes event dependence, meaning the number of events a person had in earlier time intervals does not alter the number of events this person may have in later intervals.However, recent research has shown that this assumption is not required in the estimation of the overall effect of an intervention. 12,21As previously discussed, the AG model is suitable for estimation of the overall intervention effect, and the overall effect is important from a public health perspective. 12,21,22Hence the choice we made here.One alternative approach for modeling recurrent events is the frailty model.The choice of models depends on the specific research aim.If the aim is to evaluate intervention effects by strata defined by event order or to characterize the correlation between recurrent events, the frailty model is useful. 12,41,42Future research may investigate embedding the PK/PD trajectory in the frailty model for achieving other specific research aims, for example, see the discussion in the malaria vaccine context. 12s opposed to splines or harmonic oscillator, PK/PD models are well-known in biomedical research and the parameters have biomedical interpretation.In the case study of malaria prevention, this model out-performed the alternative 4-parameter, step and spline functions.It is possible that the cumulative dose history affects the trajectory after dosing.Therefore, we fitted more flexible models that allow dose-specific parameter values.In this trial, the timing of doses was planned to coincide with immunization visits.The time interval between doses were therefore quite wide.The effect of the previous dose might have largely dissipated before another dose was given.This may explain the lack of better fit in the more complex models with dose-specific parameter values in the present study.
In this study we use the same function g 1 (t) for different doses of intervention, with or without constraining the parameter values to be the same across doses.If the context is appropriate, it is possible to modify the approach to allow different functions for capturing the impact of different components of an intervention package.One such possible context is that the intervention package includes different products at different time points, such as malaria chemoprevention and malaria vaccine.The features of different products in the intervention package may require this approach.
In simulation evaluation, the maximum partial likelihood estimators for ln(C 50 ) and  was practically unbiased and their robust sandwich variance estimators were accurate in most scenarios.The estimators for the other two parameters had some bias and the confidence interval tended to have smaller coverage probability than the nominal level, especially when sample size was 1000 and number of doses was one.Their performance improved as sample size or number of events increased within realistic range.For examples, the malaria prevention trial we re-analyzed had a larger number of participants (≅2000), doses (4) and events (≅3000) than our simulation settings.However, for studies with smaller sample size, for example, phase II trials, the parameter estimates and their confidence intervals may be inaccurate.
Regardless of sample size, the method provided accurate estimates of three of the four features of PE trajectories, that is, peak PE, t p∕2 and AUC, even in simulation scenarios (b) and (c) where the PE trajectories in the data generating process was different from that in the analysis model.Despite the larger discrepancy in the estimation of t p in scenarios (b) and (c), the absolute difference was not more than 3 days in most scenarios.In the simulation, t p was about 8 to 16 days.Vaccine trials often define 14 days after completion of primary series as the analysis start time.This implies a biomedical expectation that t p is not very far from 14 days.A deviation from true t p by a few days is unlikely a concern in most infectious disease prevention efforts.In contrast, we expect that accurate estimation of peak PE (guiding the decision to delay the second dose in the Covid-19 pandemic), t p∕2 (guiding the timing of booster vaccine dose or re-administration of prophylactic drugs), and AUC (summary measure that is comparable across trials and products) are important.The proposed method works well in this regard.
A limitation of the proposed method is that the estimation of the model needs enough amount of observation time that the PE level is close to the long-term/rebound effect , either between doses or after the final dose.Otherwise the parameter estimation can be unstable.In the malaria prevention trial and in the simulation, over 50% and at least about 30% of the follow-up time, respectively, was close to this level.This issue should be considered during study planning.
Statistical analysis plans for vaccine trials, especially for per protocol analysis, often define the starting point for evaluation of efficacy as 14 or 28 days after completion of the primary vaccination, to give enough time for the body to develop immunity.For that setting, the monotonic function in Equation ( 1) is likely more appropriate than the PK/PD function.For long-term follow-up studies, the rebound-effect parameter, , may eventually taper off.Further development of the function will be needed to capture long-term pattern beyond typical clinical trial duration.On the other hand, studies of acute adverse reactions to drug or vaccine may not expect long-term rise in adverse reaction incidence.In that case, the  term may be omitted, simplifying the function to have three parameters only.
Another limitation is that the computation is time consuming.Estimation of the PH model for the malaria prevention data in Table 1 took about 1 min on a processor with base frequency 2.2 GHz.In contrast, the proposed model with shared parameter values for all four doses took about 7 h.Breslow et al. proposed a short-cut to reduce computation time for maximum partial likelihood estimation. 43This is equivalent to treating a cohort data set as a nested case-control study data set, keeping all the events and for each event randomly select a small number of time-matched controls.Investigators may consider this approach if they have very large data sets.
While our research was motivated by and applied to the estimation of PE trajectory for infectious disease prevention trials using recurrent event data, the PK/PD function is generic and has broad applicability.Some examples include (a) using the 4-parameter function to model changes in gut microbiota diversity index (a continuous variable) after a course of antibiotic by non-linear least square method and (b) using the 3-parameter version of the function aforementioned to conduct surveillance of transient rise in the incidence of acute adverse events after receiving a drug or vaccine by self-controlled case series method.
In conclusion, the proposed method has potential to improve understanding of the trajectory of protective efficacy of preventive medical interventions and hence facilitate informed clinical and policy decision making.Further research to evaluate its performance in the studies of other products or endpoints will be useful.

F I G U R E 2
Estimated time-varying protective efficacy of four doses of a malaria prevention drug combination, assuming the timing of dosing were at 0, 1, 6 and 9 months as scheduled.Solid line: G 1 (t); dotted reference line at PE = 0.
TA B L E 1 Simulation results of parameter estimation in scenarios that used the 4-parameter PK/PD function for data generation.Simulation results of features of PE trajectory.Boldfaced are true values of the trajectory features.
TA B L E 3Note: