Comparing statistical methods for detecting and estimating waning efficacy of rotavirus vaccines in developing countries

ABSTRACT

Introduction

Vaccination has significantly reduced morbidity and mortality resulting from rotavirus infection worldwide. However, rotavirus vaccine efficacy (VE) appears to wane over the first 2 years since vaccination, particularly in developing countries. Statistical methods for detecting VE waning and estimating its rate have been used in a few studies, but comparisons of methods for evaluating VE waning have not yet been performed. In this work we present and compare three methods – Durham’s method, Tian’s method, and time-dependent covariate (TDC) method – based on generalizations of the Cox proportional hazard model.

Methods

We developed a new stochastic agent-based simulation model to generate data from a hypothetical rotavirus vaccine trial where the protective efficacy of the vaccine may vary over time. Input parameters to the simulation model were obtained from studies on rotavirus infections in four developing countries. We applied each of the methods to four simulated datasets and compared the type-1 error probabilities and the powers of the resulting statistical tests. We also compared estimated and true values of VE over time.

Results

Durham’s method had the highest power of detecting true VE waning of the three methods. This method also provided quite accurate estimates of VE in each period and of the per-period drop in VE.

Conclusions

Durham’s method is somewhat more powerful than the other two Cox proportional hazards model-based methods for detecting VE waning and provides more information about the temporal behavior of VE.

KEYWORDS:

• waning vaccine efficacy
• rotavirus vaccine
• hazard of infection
• Cox regression model
• agent-based simulations

1. Introduction

Rotavirus (RV) infections cause over 200,000 deaths annually world-wide, 85% of which occur in low-income countries of Asia and Africa.¹ RV vaccines have been used globally since 2006, and by 2020 over 100 countries have included rotavirus vaccination in their national immunization programs. There are four RV vaccines pre-qualified by the World Health Organization and available for use worldwide: GlaxoSmithKline’s (GSK’s) monovalent vaccine (Rotarix), Merck’s pentavalent vaccine (RotaTeq), Bharat Biotech’s monovalent vaccine (Rotavac), and Serum Institute of India’s pentavalent vaccine (Rotasiil).

While vaccination has significantly reduced RV-associated morbidity and mortality, there are two serious challenges related to RV vaccines in low-income settings. First, RV vaccines are less effective in low-income settings than in high-income settings.^2,3 Second, the effectiveness of RV vaccine appears to decline faster with increasing age in low-income countries than in medium- and high-income countries.² Clark et al. used a Bayesian hierarchical Poisson meta-regression model to estimate the pooled cumulative vaccine efficacy (VE) and its waning by mortality strata.² They found that, compared to low mortality settings, RV vaccine immunogenicity is lower and wanes more rapidly in high mortality settings. Rogawski et al. reached similar conclusions; they explain that waning of VE over time may be because unvaccinated children are more likely than vaccinated children to become infected and acquire natural immunity to subsequent infection (after the first infection).³ Therefore, the difference in protection between vaccinated and unvaccinated children seems to decrease over time. When VE is estimated from the cumulative incidence (attack rates) in vaccinated and unvaccinated children, and even when only the first infection of each child is considered, the estimated VE may decrease over time while the true VE remains fixed.⁴

It is very important to detect RV VE waning and estimate the rate of decrease in VE over time. Parents should be aware of the long-term protective effect of any vaccines given to their children, and health authorities should use information about VE waning when designing or modifying vaccination schedules. Over the years, a few methods have been proposed and used to evaluate VE waning for various vaccines. However, we could not find any study that compares these methods or recommends which method to evaluate VE waning should be preferred.

The main goal of this work is to compare different statistical methods to detect true waning of RV VE and to estimate the magnitude of age-related changes in VE. These comparisons are performed using new stochastic agent-based simulation software that generates data like those observed in real-life RV VE studies. Input parameters are based on RV studies conducted in developing countries. We focus on VE against the first occurrence of rotavirus diarrhea (RVD) to avoid confounding of natural immunity acquired after the first infection with the protective effect of the vaccine.³ To our best knowledge, this is the first article comparing different statistical methods for estimation of waning VE.

2. Materials and methods

2.1. Estimation of vaccine efficacy over time

Vaccine efficacy is usually estimated as one minus the ratio of the attack rates (cumulative incidence rates) over a specific time interval between vaccinated and unvaccinated persons. Since this a summary estimate of VE over a time interval, it does not reflect the magnitude of VE at a given time point. Moreover, this estimate of VE may decrease over time even when there is no temporal change of the true VE.⁴ Therefore, we recommend to estimate VE at a given time point as one minus the ratio of the hazard of the outcome of interest (first occurrence of RVD in our case) in a vaccinated and an unvaccinated person at this time point. Note that the hazard at a given time point is defined as the instantaneous risk of the outcome in persons who have not experienced the outcome (i.e. are still susceptible) just before that time point.

2.2. The simulation model and input data

Our stochastic agent-based simulation model assumes that a cohort of two thousand 60-day-old infants is followed up until they reach the age of 360 days. We assume that a pre-set proportion of the children in the cohort received all required doses of the rotavirus vaccine 14 days prior to the onset of the study. Vaccination is done at random. The study duration (300 days) is divided into 10 periods of 30 days each. The input parameters for each period are as follows: (1) $\lambda$, which denotes the daily hazard of an unvaccinated susceptible child to contract her/his first RVD occurrence, and (2) θ, which is defined such that the daily hazard of RVD in a vaccinated susceptible child is $\lambda \cdot \theta$. Hence, the true vaccine efficacy is VE = 1- θ. The values of these daily parameters are fixed during each period but may vary from one period to the next.

The input values for the daily hazards ($\lambda$’s) of the first occurrence of RVD in unvaccinated children were obtained from MAL-ED, a multi-country birth cohort study on enteric infections and malnutrition in 8 countries from April 2009 to February 2014.⁵ Healthy children were enrolled within 17 days of birth and were followed until 24 months of age. We used data from four countries where RV vaccination was not yet included in the national immunization program: Bangladesh, India, Nepal and Pakistan. For the θ ’s we assumed that at the beginning of the study VE was 60%,² and we considered 3 levels of vaccine efficacy waning: (1) no waning (i.e. θ was fixed over all periods), (2) moderate waning, where we let the absolute VE decrease 3% in each 30-day period, and (3) severe waning, where the absolute VE decreased 6% per period.

For each combination of country and level of waning we conducted 100 simulations. The output from each simulation was an ‘outcomes file’ which included each child’s vaccination status and the age (in days) at the first RVD episode (if any). Three methods (described below) were used to test the null hypothesis that the efficacy of the RV vaccine does not change with time and to estimate the rate of waning. These methods were applied to the outcomes file from each simulation, and the proportion of simulations where the null hypothesis of no waning was not accepted was determined. When there is no waning, this proportion estimates the type I error probability of the statistical test. When there is waning, this proportion is an estimate of the power of the test. In addition, one of the methods was used to estimate the value of VE in each period and the rate of change in VE from one period to the next.

A copy of the simulation software is available from the authors upon request.

2.3. The Cox Proportional hazard model and its generalizations

Since we defined VE as one minus the hazard ratio comparing vaccinated and unvaccinated children, the hypothesis that VE does not change with age is equivalent to the hypothesis of a fixed (time-independent) hazard ratio, or equivalently, the proportional hazard hypothesis. The most common estimate of a fixed hazard ratio is based on Cox’s proportional hazard model. Let t denote the time (in days) from vaccination and let V denote the binary vaccination status. Under the Cox model, the hazard of a child of vaccination status V at time t is written as: $\lambda \left(t|V\right)={\lambda }_0\left(t\right)\cdot Exp\left\{\beta \cdot V\right\}$, where λ₀(t) is the hazard of an unvaccinated child at time t (the baseline hazard), and β is a coefficient that can be estimated from the data. Then the fixed hazard ratio is HR = Exp(β) and the fixed VE is estimated as one minus the estimate of Exp(β). Therefore, changes in VE over time violate the proportional hazard assumption. To allow for time-variable VE, the Cox model presented above has to be replaced by a more general model, such as: $\lambda \left(t|V\right)={\lambda }_0\left(t\right)\cdot Exp\left\{\beta \left(t\right)\cdot V\right\}$. In other words, the coefficient β is allowed to vary over time. Then VE at time t is VE(t) = 1−Exp{β(t)}, and the null hypothesis of no waning is equivalent to β(t) ≡ β for all time points t. We consider three methods for testing this hypothesis.

2.4. Durham’s method

This method is based on smoothing scaled residuals from the proportional hazard model.⁶ It consists of four steps. First, an ordinary proportional hazard model is fitted. Second, Schoenfeld residuals⁷ are calculated. The residual is the difference between the covariate at the failure time and the expected value of the covariate at this same time. Third, these residuals are expanded about the regression coefficients in a Taylor series and added to the coefficient from the ordinary proportional hazard model. Fourth, after smoothing the estimated proportional hazard coefficient and the scaled Schoenfeld residuals, we can get the estimated hazard ratios as function of time. This allows the estimation of vaccine efficacy at each time point.

2.5. Tian’s method

In this method, a kernel-weighted partial likelihood approach is used to obtain a confidence band of the time-dependent coefficient $\beta \left(t\right)$.⁸ At each time point, the estimate is obtained by maximizing a smooth concave function of a p x 1 vector of parameters, where p is the dimension of the vector of covariates. The $\left(1-\alpha \right)$ confidence bands for the time-dependent coefficient $\beta \left(t\right)$ can be obtained by

$\hat{\beta }\left(t\right)\pm c_{\alpha }\hat{w}{\left(t\right)}^{-1}$, where t is in the time interval of interest $\left(b_1\le t\le b_2\right)$, $c_{\alpha }$ is the $100\left(1-\alpha \right)$th percentile of the approximate distribution of $\hat{\beta }\left(t\right),$and $\hat{w}\left(t\right)$ is a positive weighting function. These confidence bands can also be used to test the hypothesis of no VE waning (i.e., that the proportional hazards assumption is met).

2.6 Time-dependent covariates (TDC) method

In this method, the following Cox proportional hazard model with a time-dependent covariate is fit: $\lambda \left(\mathrm{t}|\mathrm{V}\right)={\lambda }_0\left(\mathrm{t}\right)\cdot \mathrm{E}\mathrm{x}\mathrm{p}\left\{\beta \cdot \mathrm{V}+\gamma \cdot \mathrm{V}\cdot \mathrm{g}\left(\mathrm{t}\right)\right\},\hspace{0.17em}$where $g\left(t\right)$ is an arbitrary function of time; usually the function $g\left(t\right)=\mathrm{l}\mathrm{o}\mathrm{g}\left(t\right)$ is used.⁹ This is a proportional hazard model with two covariates: the time-fixed binary vaccination status V and the time-dependent covariate $V\cdot g\left(t\right).$ Under this model, vaccine efficacy at time t is: $VE\left(t\right)=1-exp\left\{\beta +\gamma \cdot g\left(t\right)\right\}.$ Therefore, the hypothesis of no vaccine efficacy waning is equivalent to H₀: $\gamma =0$.

For each of the methods listed above (Durham’s, Tian’s, and TDC) we determined the proportion of simulations where the null hypothesis of no VE waning was rejected under each of the three waning scenarios.

3. Results

3.1. Comparison of proportion rejections of the null hypothesis of no VE waning between the three methods

Table 1 presents the proportion of simulations where the null hypothesis of no VE waning was rejected under each of the three waning scenarios (Section 2.2) for each of the three methods (Durham’s, Tian’s, and TDC). In the no waning scenario, this proportion should be around or below the nominal significance level of 0.05. This was satisfied with all three methods in the simulations that used the hazard rates from Bangladesh, India and Nepal. In the simulations that used the hazard rates from Pakistan, the proportions of rejections in the no waning scenario exceeded 0.05 with all 3 methods (0.08, 0.07, 0.08, respectively).

Table 1. Proportion of simulations where the hypothesis of no VE waning was rejected

	No waning			Moderate waning			Severe waning
	Durham	Tian	TDC	Durham	Tian	TDC	Durham	Tian	TDC
Bangladesh	0.01	0.03	0.02	0.23	0.20	0.18	0.63	0.55	0.55
India	0.04	0.02	0.03	0.18	0.14	0.13	0.25	0.21	0.22
Nepal	0.04	0.06	0.04	0.17	0.15	0.08	0.36	0.31	0.34
Pakistan	0.08	0.07	0.08	0.23	0.21	0.19	0.54	0.48	0.48

Next, we look at the results for the moderate and severe waning scenarios, where the true VE is assumed to decrease by 3 and 6% in each 30-days period, respectively. Under these scenarios, where the proportion of rejections estimates the power of the test of no waning, the proportion of rejections with Durham’s method was somewhat higher than with the other two methods based on the data from all four countries. For example, for the simulations based on Bangladesh’s hazard rates, under the severe waning scenario the proportions of rejections with Durham’s, Tian’s and the TDC methods were 0.63, 0.55 and 0.55, respectively.

3.2. Accuracy of the VE estimates

Of the three methods described earlier, only Durham’s method provides estimates of VE in each period. The means of estimated VEs over 100 simulations for each period using Durham’s method are shown in Table 2 for the simulations based on the Bangladesh hazard rates. The estimated VEs were usually close to the true VEs used in the simulations. For example, the absolute biases (estimated VE minus true VE) in period 3 for the no waning, moderate waning and severe waning scenarios were 0.4%, −0.5% and −1.3%, respectively. In period 6, the biases were −0.1%, −4.2%, and 1.8%, respectively, and in period 9 the biases were 0.0%, −5.0%, and −2.3%, respectively. Similar results (not shown) were observed for the simulations based on the hazard rates from the other three countries.

Table 2. Comparing estimated VE via Durham’s method with true VE from simulations based on Bangladesh’s hazard rates

	No waning		Moderate waning		Severe waning
Period	Estimated VE(%)	True VE(%)	Estimated VE(%)	True VE(%)	Estimated VE(%)	True VE(%)
1	60.4	60.0	57.9	60.0	53.5	60.0
2	60.3	60.0	56.1	57.0	50.7	54.0
3	60.4	60.0	53.5	54.0	46.7	48.0
4	60.1	60.0	51.3	51.0	41.9	42.0
5	60.0	60.0	48.9	48.0	36.3	36.0
6	59.9	60.0	45.8	45.0	31.8	30.0
7	60.1	60.0	41.7	42.0	25.4	24.0
8	59.9	60.0	36.2	39.0	18.1	18.0
9	60.0	60.0	31.0	36.0	9.7	12.0
10	60.9	60.0	26.8	33.0	1.7	6.0

3.3. Accuracy of estimated drop in VE per period

We compared the average estimated per-period drop in VE using Durham’s method to the true per-period drop used in the input to the simulation software. The results are shown in Table 3. The average change in VEs per period was close to 0 in the no waning scenario. For the moderate and severe waning scenarios, the estimated changes in VE per period were close to the true value of −0.03 and −0.06, respectively.

Table 3. Estimated change in absolute VE per period (Δ) using Durham’s method

	No waning $\mathrm{\Delta }=0\mathrm{\%}$	Moderate waning $\mathrm{\Delta }=-3\mathrm{\%}$	Severe waning $\mathrm{\Delta }=-6\mathrm{\%}$
Bangladesh	0.1%	−3.5%	−5.8%
India	0.1%	−2.7%	−5.4%
Nepal	0.2%	−2.8%	−4.9%
Pakistan	0.8%	−2.9%	−6.1%

4. Discussion

In our comparison of three different methods based on generalizations of the Cox proportional hazard model to evaluate waning of RV VE, we found that all three methods can detect drops in VE over time when applied to simulated data with true decreasing VE. Durham’s method is somewhat more likely to reject the null hypothesis of no VE waning, i.e., it has the highest power of detecting true VE waning. Durham’s method also provides quite accurate estimates of VE in each period, and of the per-period drop in VE. Tian’s and the TDC method do not provide estimates of VE in each period and cannot be used to estimate the rate of VE waning. Output from Durham’s method can also be used to plot the estimated VEs as a function of time since vaccination. These plots provide the user with visual information about the temporal behavior of VE.

Based on these comparisons, we conclude that Durham’s method is somewhat more powerful than the other two Cox proportional hazards model-based methods in terms of detecting VE waning and provides more information about the temporal behavior of VE. However, Durham’s method has a few limitations. First, though Durham’s method has somewhat higher power for rejecting the hypothesis of no waning, compared to the other two methods, the power is still a bit low given the relatively large sample size of 2000. Even in the severe waning scenario, when VE declines by 6% in each 30-day period, the power we observed in this simulation study ranged between 0.25 and 0.63. Under moderate waning, where the true decline in VE was 3% per period, the power ranged between 0.17 and 0.23. Hence, future research will be needed to develop more powerful methods for detecting VE waning. Second, this method is based on the assumption that the hazard of infection follows the (generalized) Cox regression model, i.e., the hazard can be written as $\lambda \left(t|V\right)={\lambda }_0\left(t\right)\cdot Exp\left\{\beta \left(t\right)\cdot V\right\}$, where $t$ is the time since vaccination and $V$ is the binary vaccination status. This assumption about the hazard function may not always hold. Third, this method assumes that vaccination status does not vary over time. In contrast, Tian’s and the TDC method allows persons to receive the vaccine at different times or ages. More research is needed to develop methods that can be applied to situations where study participants are vaccinated at different times points. In addition, it may be possible to develop more general methods that allow researchers to use information on temporal changes in the prevalence of infection in the population by modeling the hazard of infection as a function of the prevalence.

Our study has a few limitations. First, our model assumes that all vaccinated study participants received the vaccine at the same time. Second, we assumed a randomized vaccine trial or a cohort study where one does not adjust for possible confounders. Third, we only used input data from developing countries. Finally, we used only a single value of true VE = 60% (we repeated some of the simulations with a higher VE of 80% and observed similar results). We plan to address these issues in the future. We also plan to use data on other vaccine-preventable diseases, such as influenza, pneumococcal diseases, and COVID-19.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Additional information

Funding

This work was supported by the Centers for Disease Controls and Prevention (CDC) via an IPA [19IPA1912112]. The findings and conclusions in this report are those of the author(s) and do not necessarily represent the official position of the Centers for Disease Control and Prevention.