Comparison of three sample size calculation methods for non-inferiority vaccine trials with multiple continuous co-primary endpoints

ABSTRACT

Clinical trials that study immunogenicity of combination vaccines often have less power than desirable. To make up for the reduction in statistical power at the study level, researchers have to increase the study sample size. To study immunogenicity variables, we used the geometric mean concentration of immune response after vaccination as immunologic endpoint and compared 3 sample size calculation methods: the “Inflation factors” method, the “Incrementing” method, and the “Bonferroni correction” method when there are multiple continuous co-primary endpoints. The parameters were set according to the actual situation of the use of combination vaccines and the simulation results were used as reference. The present study demonstrates that all 3 methods are applicable when the effect size of each endpoint is similar and the endpoints are at most weakly correlated, but when there is a true difference in effect sizes among endpoints, the “Incrementing” method has the best performance.

Keywords:

• co-primary endpoints
• sample size calculation
• overall power
• combination vaccine
• intersection-union

Introduction

Combination vaccines^1,2, which are composed of 2 or more inactivated organisms, live organisms or purified antigens, have been extensively used in recent decades as they offer protection against multiple diseases or several subtypes of bacteria or viruses and help simplify the current immunization schedule³. Owing to the multiple components in these vaccine products, they require multiple co-primary immunogenicity endpoints; for example, pneumococcal vaccination often requires 13 or 23 co-primary immunogenicity endpoints. Other multi-strain vaccines include those for human papillomavirus, rotavirus, influenza, meningococcus and poliomyletitis. Accordingly, the quality of such products might be determined by several parameters⁴. To demonstrate the efficacy of these vaccines⁵, regulatory agencies often require all endpoints to be significant at the prespecified significance level. Therefore, the hypotheses regarding each serotype or component are referred to as co-primary endpoints. When the purpose of the trial is to ascertain the effectiveness of a vaccine on all of the designated co-primary variables, there is no need to adjust type I error rate. However, as described in “ICH Harmonised Tripartite Guideline for Statistical Principles Clinical Trials E9”⁶ and “Committee for Human Medicinal Products Guideline on Multiplicity Issues in Clinical Trials”⁷, the influence of co-primary endpoints on type II error rate should be thoroughly considered^8,9. To compensate for the loss of study-wide statistical power, researchers should increase the sample size. The natural question is, “what is the appropriate sample size?”

The Center for Drug Evaluation and Research (CDER) Guidance for Industry for the Evaluation of Combination Vaccines for Preventable Diseases¹ recommends that the sample size required for each endpoint should be calculated separately and the largest sample size should be selected as the sample size for the overall trial. This method of calculating sample size is extensively accepted and used currently. By using this method, the loss of overall power may not be so severe when the endpoints are correlated or widely differ in effect size and variability estimates. When the endpoints are entirely independent, the overall power to attain success will be reduced to the the number of co-primary endpoints power of the prespecified study-wide statistical power. Where there are 4 mutually independent endpoints, if the type II error rate is 0.2 for each endpoint, then the overall power will be 0.8⁴ =0.41. Thus, the likelihood of the study failing to support a conclusion of a favorable vaccine effect when such an effect exists would be 59.04%. Data extracted from previously published studies suggest that the high and low correlation levels between co-primary endpoints can range from 0.2 to 0.8 in clinical trials⁹. In vaccine trials, since the range of correlations between immunogenicity endpoints is still not clear, detailed discussion based on a wide range of correlation is worthwhile. Therefore, we investigated the actual data of phase II and phase III vaccine clinical trials for the purpose of summarizing the relevant patterns between endpoints and revealing the impact that multiple co-primary immunogenicity endpoints have on sample size and overall power.

Inflation of the Type II error rate in trials based on multiple co-primary endpoints has been studied^8,10–15. Xiong et al.¹⁰ explored sample size calculation in a study with 2 normally distributed continuous endpoints. Sozu et al.¹¹ extended the sample size calculation formula to a study with 3 correlated endpoints^8,12. Sugimoto et al.¹⁵ took correlations between endpoints and the ratios of effect size in multiple co-primary outcomes into consideration. These previous reports mostly focused on superiority trials with 2 or 3 co-primary endpoints. For non-inferiority clinical trials of combination vaccines when there are more than 2 endpoints, there are some other methods^16,17 which might be more practical. Julious et al.¹⁷ provided a method of calculating sample size required for the entire study by multiplying the sample size for a single endpoint by the inflation factor. This method extends the results provided by Sankoh that explains both the number of comparisons and the correlations between endpoints for the correction of Type II error rate from five endpoints to seven endpoints. Varga et al.¹⁶ presents another method which allocates more power to the groups that may have larger standardized effect sizes based on prior information. The premise of this method is the independence of endpoints. Besides these, “Bonferroni correction” is still the most popular approach used for correcting the error rate.

As far as we know, there is no previous report on sample size calculation for trials with multiple co-primary endpoints considering the nature of immunogenicity variables. In our study, we used the geometric mean concentration (GMC) of immune response after vaccination as the immunologic endpoint and discussed the following 3 methods for sample size calculation: the “Inflation factors” method, the “Incrementing” method, and the “Bonferroni correction” method.

Methods

Considering a non-inferiority vaccine trial, the objective is to determine whether the efficacy of the experimental vaccine is not worse than that of the active comparator concerning the null hypothesis. When GMC is the outcome measured in a randomized clinical vaccine trial, then the parameter of interest is the ratio of geometric means. The trial is considered to have shown benefit for all components of the vaccine if the lower limits of the confidence interval for the ratio of the experimental vaccine group to the control vaccine group is greater than the pre-specified non-inferiority margin for all endpoints.

We start with some notation. Suppose there are K (≥ 2) co-primary outcomes in a study with n_T subjects in the experimental group and n_C subjects in the active control group. Let V denote the post-vaccination humoral and cellular immunogenicity value, ln(V_i) (i = 1. K) are normally distributed with mean μ_i, variance σ_i², and GMC =$e^{{\mu }_i}$. The means μ_Ti, μ_Ci, the variances σ²_Ti, σ²_Ci, and the correlation coefficients ρ_Tij, ρ_Cij (i = 1. K, j = 1. K, i ≠ j) among endpoints are specified in advance. The vector of log-transformed immunogenicity response Y_j = (Y_j1, … Y_jK)^T (j = 1. n) obeys K-variate normal distributions^8,18 with mean vector E[Y] = μ = (μ₁. μ_K)^T, and common covariance matrix Σ, i.e.,

$Y_j\sim {\mathrm{N}}_K\left(\mu ,\mathbf{\Sigma }\right),$

where

$\mathbf{\Sigma }=\left[\begin{array}{ccc}{\sigma }_1^2& \begin{array}{ccc}\cdots & {\rho }_{1j}{\sigma }_1{\sigma }_j& \cdots \end{array}& {\rho }_{1K}{\sigma }_1{\sigma }_K\\ \begin{array}{c}⋮\\ {\rho }_{i1}{\sigma }_i{\sigma }_1\\ ⋮\end{array}& \ddots & \begin{array}{c}⋮\\ {\rho }_{iK}{\sigma }_i{\sigma }_K\\ ⋮\end{array}\\ {\rho }_{K1}{\sigma }_K{\sigma }_1& \begin{array}{ccc}\cdots & {\rho }_{Kj}{\sigma }_K{\sigma }_j& \cdots \end{array}& {\sigma }_K^2\end{array}\right].$

For a single endpoint, the null (H₀) and alternative (H₁) hypotheses take the following forms:

H₀: The experimental vaccine is inferior to the control vaccine (μ_T-μ_C≤-d).

H₁: The experimental vaccine is non-inferior to the control vaccine (μ_T-μ_C >-d).

When defining the null and alternative hypotheses, d is known as the non-inferiority limit. The lower limit of the confidence interval for the difference between μ_T and μ_C should be greater than -d, which is the same as requiring the lower limit of the confidence interval for the ratio of the geometric means $e^{{\mu }_T}$ and $e^{{\mu }_C}$ greater than λ, where d = ln(1/λ). A commonly used margin^19-22 for the ratio is ln(1/2), which means that the experimental vaccine be at least half as efficacious as the control vaccine.

We first introduce the sample size calculation procedures we used for the non-inferiority trials and then describe the design schemes that assume the endpoints of interest to be following lognormal distribution. All the calculations and simulations were conducted using SAS (ver 9.4). The related SAS code to implement the “Incrementing” method is provided in Supplementary 3.

Sample size calculation procedure

For a trial with a single endpoint, many methods can be used to calculate the sample size for non-inferiority parallel group trials with normal data^23,24. Among them, the most commonly used formula is proposed by Julious et al.²³ as,

(1) $n_T=\frac{\left(r+1\right){(Z_{1-\beta }+Z_{1-\alpha /2})}^2{\sigma }^2}{r{\left(\left({\mu }_T-{\mu }_C\right)-d\right)}^2}$(1)

Here, r is the allocation ratio between groups (n_T = rn_C), μ is the mean, σ² is the variances, d is the non-inferiority limit, and Z_1-α/2 and Z_1-β are the corresponding values for type I and type II error risks α and β, respectively, assuming a standard normal distribution.

From formula (1), if we set the allocation ratio to 1 for maximum study power, the factors that determine the sample size can be summarized as type I and II error rates, non-inferiority limit, and the effect sizes. Generally, the non-inferiority limit and the type I and II error rates are fixed. Factors that can have a direct effect on sample size are $\left({\mu }_T-{\mu }_C\right)-d$ and the variance $\sigma$². For trials with multiple co-primary endpoints, denote by c_ij the sample size ratio of the sample size calculated by formula (1) for endpoint i to that for endpoint j (i = 1… K, j = 1… K, i ≠ j), which is,

(2) $c_{ij}=\frac{n_{Ti}}{n_{Tj}}$(2)

Here, we are interested in performing hypothesis tests on the difference in GMC between the experimental arm and the control arm. Assuming equal variance and identical samples of the 2 arms, let diff = μ_T-μ_C. According to formula (1), c_ij could be simplified as,

(3) $c_{ij}=\frac{{\left(dif{f}_j-d_j\right)}^2}{{\left(dif{f}_i-d_i\right)}^2}$(3)

Because different study endpoints may choose different non-inferiority thresholds, the choice of d_j and d_i must always be justified on both clinical and statistical grounds.

To calculate the sample size for trials with 2 co-primary endpoints, the “Jumping method”²⁵ is frequently used. With this method, the conjunctive power is (1-β₁)*(1-β₂), if 1-β₁ and 1-β₂ represent the 2 power values assigned to the 2 endpoints. The sample size is iterated from small to large and the power of the 2 endpoints is calculated using the same sample size. The smallest sample size that makes the conjunctive power meet expectations is the sample size of the entire study. For trials with more than 2 co-primary endpoints, an improved method, the “Incrementing” method can be adopted¹⁶. With the assumption that the endpoints are mutually independent, this method assigns the power of $\sqrt[k]{1-\beta }$ to each endpoint and calculates the corresponding sample sizes. The smallest sample size among endpoints is used as the initial value for the iteration and the power for each endpoint is recalculated. The sample size is increased by one and the power for each endpoint is recalculated if the overall power does not reach the expected value. The iteration stops when the overall power is no less than 1-β,

(4) $1-\beta =\prod \left({1-\beta }_{\mathrm{i}}\right)\left(\mathrm{i}=1,...,\mathrm{K}\right)$(4)

This design assigns a larger power to an endpoint that has a larger effect size, which provides an efficient way of using the sample size. Another popular approach is to use Bonferroni corrections²⁶. In order to control the family-wise type II error rate, the type II error rate for each endpoint is corrected as β/K. Although this method is conservative especially when the comparisons are correlated, it is simple and practical. The third approach of calculating the sample size for trials with multiple co-primary endpoints is provided by Julious et al.^17,23. When the multiple co-primary endpoints are potentially correlated, the sample size for a single endpoint in clinical trials with continuous co-primary endpoints is increased by using the “Inflation factors” method. This method is an extension of the result of Sankoh²⁷, which considers the number of comparisons and the correlations between endpoints for the type II error rate as follows,

(5) ${\beta }_k=1-\left(\sqrt[t]{1-c\beta }\right)$(5)

where,

(6) $t=k^{1-\rho }$(6)

and k is the number of comparisons, ρ is the correlation coefficient between the comparisons, β is the prespecified type II error rate, β_k is the value to calculate the sample size for each arm as described in methods for single endpoint, and the values for c are described by Julious, where K ≤ 7. Thus, the increase in sample size for must-win trials compared with a single endpoint trial can be described as,

(7) $\mathrm{I}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{F}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r}=\frac{{\left(Z_{1-\alpha /2}+Z_{1-{\beta }_k}\right)}^2}{{\left(Z_{1-\alpha /2}+Z_{1-\beta }\right)}^2}$(7)

The “Inflation factor” method uses a single correlation coefficient instead of the correlation coefficient matrix between the endpoints. This assumption of equal correlation may not be true; however, it represents a pragmatic solution to the problem if the correlations are not known to differ greatly. Supplementary Table 2 shows that the coefficient of variation of the correlations between the endpoints is no more than 33%. In practice, it might be difficult to set a detailed correlation matrix before the study; therefore, a minimum estimate of correlation between endpoints should be chosen to be conservative.

Table 1. Sample size estimated to acquire an overall power of 80% or 90% at significant level α = 0.05 (2-tailed) by simulation. The proportion of sample size decreases as the correlation coefficient increases compared to r = 0 were provided in parentheses. For example, when c_ij = 1, r = 0.3 and K = 2, the proportion of sample size decreases was 100%×(45–44)/45 = 2.22%.

	Sample Size
	80% Overall Power				90% Overall Power
Correlation Coefficient	2 endpoints	5 endpoints	10 endpoints	20 endpoints	2 endpoints	5 endpoints	10 endpoints	20 endpoints
cij = 1
0	45	58	66	78	54	69	78	88
0.3	44(2.22%)	54(6.9%)	62(6.06%)	74(5.13%)	53(1.85%)	67(2.9%)	72(7.69%)	84(4.55%)
0.6	42(6.67%)	51(12.96%)	57(14.52%)	67(14.86%)	51(5.66%)	63(8.96%)	68(13.89%)	79(10.71%)
0.9	38(15.56%)	43(29.41%)	46(35.09%)	49(43.28%)	48(11.76%)	54(23.81%)	60(26.47%)	63(31.65%)
cij = 1.5
0	55	66	72	82	67	78	85	92
0.3	54(1.82%)	63(4.55%)	68(5.56%)	77(6.1%)	67(0%)	77(1.28%)	81(4.71%)	90(2.17%)
0.6	53(3.7%)	57(14.29%)	63(13.24%)	69(16.88%)	66(1.49%)	74(5.19%)	77(9.88%)	83(10%)
0.9	50(9.43%)	52(24.56%)	54(28.57%)	54(40.58%)	65(3.03%)	70(10.81%)	69(20.78%)	71(25.3%)
cij = 3
0	95	95	104	106	127	129	138	132
0.3	95(0%)	95(0%)	102(1.92%)	104(1.89%)	127(0%)	129(0%)	138(0%)	132(0%)
0.6	94(1.05%)	94(1.05%)	100(3.92%)	102(3.85%)	127(0%)	129(0%)	138(0%)	132(0%)
0.9	94(1.06%)	94(1.06%)	100(4%)	101(4.9%)	127(0%)	129(0%)	138(0%)	132(0%)

The sample size shown in the table is the required sample size for a single group.

Table 2. Sample size calculated by these three methods and by simulation to acquire an overall power of 80% or 90%. The number of endpoints were 2, 5, 10 and 20, the correlation between endpoints were 0, 0.3, 0.6 and 0.9 and the c_ij were 1, 1.5 and 3, respectively. The ratio of the difference between sample sizes calculated by these three methods and the sample size estimated by simulation compared to the sample size estimated by simulation were provided in parentheses.

		Sample Size
		80% Overall Power				90% Overall Power
Number of Endpoints	Correlation	Simulation	Incrementing	Bonferroni	Inflation Factors	Simulation	Incrementing	Bonferroni	Inflation Factors
cij = 1
2	0	45	43(−4.44%)	43(−4.44%)	42(−6.67%)	54	54(0%)	54(0%)	53(−1.85%)
	0.3	44	43(−2.27%)	43(−2.27%)	41(−6.82%)	53	54(1.89%)	54(1.89%)	52(−1.89%)
	0.6	42	43(2.38%)	43(2.38%)	41(−2.38%)	51	54(5.88%)	54(5.88%)	51(0%)
	0.9	38	43(13.16%)	43(13.16%)	37(−2.63%)	48	54(12.5%)	54(12.5%)	48(0%)
5	0	58	57(−1.72%)	57(−1.72%)	55(−5.17%)	69	67(−2.9%)	67(−2.9%)	65(−5.8%)
	0.3	54	57(5.56%)	57(5.56%)	54(0%)	67	67(0%)	67(0%)	64(−4.48%)
	0.6	51	57(11.76%)	57(11.76%)	51(0%)	63	67(6.35%)	67(6.35%)	61(−3.17%)
	0.9	43	57(32.56%)	57(32.56%)	45(4.65%)	54	67(24.07%)	67(24.07%)	55(1.85%)
10	0	66	66(0%)	67(1.52%)		78	76(−2.56%)	76(−2.56%)
	0.3	62	66(6.45%)	67(8.06%)		72	76(5.56%)	76(5.56%)
	0.6	57	66(15.79%)	67(17.54%)		68	76(11.76%)	76(11.76%)
	0.9	46	66(43.48%)	67(45.65%)		60	76(26.67%)	76(26.67%)
20	0	78	76(−2.56%)	76(−2.56%)		88	86(−2.27%)	85(−3.41%)
	0.3	74	76(2.7%)	76(2.7%)		84	86(2.38%)	85(1.19%)
	0.6	67	76(13.43%)	76(13.43%)		79	86(8.86%)	85(7.59%)
	0.9	49	76(55.1%)	76(55.1%)		63	86(36.51%)	85(34.92%)
cij = 1.5
2	0	55	55(0%)	66(20%)	64(16.36%)	67	71(5.97%)	81(20.9%)	81(20.9%)
	0.3	54	55(1.85%)	66(22.22%)	63(16.67%)	67	71(5.97%)	81(20.9%)	80(19.4%)
	0.6	53	55(3.77%)	66(24.53%)	62(16.98%)	66	71(7.58%)	81(22.73%)	79(19.7%)
	0.9	50	55(10%)	66(32%)	57(14%)	65	71(9.23%)	81(24.62%)	74(13.85%)
5	0	66	64(−3.03%)	86(30.3%)	84(27.27%)	78	77(−1.28%)	101(29.49%)	100(28.21%)
	0.3	63	64(1.59%)	86(36.51%)	83(31.75%)	77	77(0%)	101(31.17%)	98(27.27%)
	0.6	57	64(12.28%)	86(50.88%)	78(36.84%)	74	77(4.05%)	101(36.49%)	94(27.03%)
	0.9	52	64(23.08%)	86(65.38%)	69(32.69%)	70	77(10%)	101(44.29%)	84(20%)
10	0	72	71(−1.39%)	101(40.28%)		85	83(−2.35%)	115(35.29%)
	0.3	68	71(4.41%)	101(48.53%)		81	83(2.47%)	115(41.98%)
	0.6	63	71(12.7%)	101(60.32%)		77	83(7.79%)	115(49.35%)
	0.9	54	71(31.48%)	101(87.04%)		69	83(20.29%)	115(66.67%)
20	0	82	79(−3.66%)	115(40.24%)		92	90(−2.17%)	129(40.22%)
	0.3	77	79(2.6%)	115(49.35%)		90	90(0%)	129(43.33%)
	0.6	69	79(14.49%)	115(66.67%)		83	90(8.43%)	129(55.42%)
	0.9	54	79(46.3%)	115(112.96%)		71	90(26.76%)	129(81.69%)
cij = 3
2	0	95	98(3.16%)	129(35.79%)	126(32.63%)	127	130(2.36%)	160(25.98%)	159(25.2%)
	0.3	95	98(3.16%)	129(35.79%)	124(30.53%)	127	130(2.36%)	160(25.98%)	156(22.83%)
	0.6	94	98(4.26%)	129(37.23%)	122(29.79%)	127	130(2.36%)	160(25.98%)	154(21.26%)
	0.9	94	98(4.26%)	129(37.23%)	112(19.15%)	127	130(2.36%)	160(25.98%)	144(13.39%)
5	0	95	99(4.21%)	169(77.89%)	165(73.68%)	129	130(0.78%)	198(53.49%)	196(51.94%)
	0.3	95	99(4.21%)	169(77.89%)	162(70.53%)	129	130(0.78%)	198(53.49%)	192(48.84%)
	0.6	94	99(5.32%)	169(79.79%)	154(63.83%)	129	130(0.78%)	198(53.49%)	184(42.64%)
	0.9	94	99(5.32%)	169(79.79%)	134(42.55%)	129	130(0.78%)	198(53.49%)	165(27.91%)
10	0	104	100(−3.85%)	198(90.38%)		138	130(−5.8%)	226(63.77%)
	0.3	102	100(−1.96%)	198(94.12%)		138	130(−5.8%)	226(63.77%)
	0.6	100	100(0%)	198(98%)		138	130(−5.8%)	226(63.77%)
	0.9	100	100(0%)	198(98%)		138	130(−5.8%)	226(63.77%)
20	0	106	103(−2.83%)	226(113.21%)		132	131(−0.76%)	253(91.67%)
	0.3	104	103(−0.96%)	226(117.31%)		132	131(−0.76%)	253(91.67%)
	0.6	102	103(0.98%)	226(121.57%)		132	131(−0.76%)	253(91.67%)
	0.9	101	103(1.98%)	226(123.76%)		132	131(−0.76%)	253(91.67%)

The sample size shown in the table is the required sample size for a single group.

Table 3. Factors considered in the sample size calculation by these mehtods.

Methods	Type I and Type II error rates	Numbers of endpoints	Differences in Effect Size among Endpoints	Correlation among Endpoints
Incrementing	Yes	Yes	Yes	No
Bonferroni	Yes	Yes	No	No
Inflation Factors	Yes	Yes	No	Yes

Simulation strategy and methods

In order to validate the methods mentioned above and to gain further insight into statistical power and correlation, the sample size calculated was compared to the sample size obtained by simulation, used as a reference. For simplicity, the allocation ratio r is set to 1 by assuming n_T = n_c, which is the most common case.

Simulation methods

The following steps were involved in the simulation of sample size:

Step 1: The log-transformed immunogenicity response mean vector E[Y] = μ = (μ₁…μ_k)^T, common covariance matrix Σ, and an initial sample size n is used to generate K sets of normally distributed data of 2 groups of subjects.

Step 2: Based on the simulation data, the hypothesis of non-inferiority is further tested by ANOVA and the lower limit of the 95% confidence of the difference between 2 groups is obtained.

Step 3: If all lower limits of the 95% confidence of the K co-primary endpoints are greater than the pre-defined non-inferiority limit, then the trial is marked as “success”, otherwise, the trial is labeled as “failure”.

Step 4: Steps 1–3 are repeated 1000 times and the number m of “successes” is counted.

Step 5: For the sample size as n, the overall power is m/1000.

Step 6: The sample size is slightly increased and steps 1–5 are repeated until the maximum sample size is reached.

In this study, the maximum value of the sample size calculated by formula (1) for each endpoint when the significance level is α and the predetermined power of each study is 1−β was used as an initial value in simulation. The sample size calculated with Bonferroni corrected 1-β_K (β_K = β/K) was used to determine the maximum sample size required to achieve the required overall power 1-β. Note that we have to repeatedly generate simulation data and testing hypothesis in order to understand the performance of these methods in different situations; at least 100 repeats for each situation should be used to endure accuracy. We have used 1000 simulations to achieve appropriate precision.

Design scheme

In order to capture the nature of immunogenicity variables and design the parameters for simulation, we analyzed the raw data of 4 combination vaccine trials. GMC, standard deviation of GMC of each immunologic endpoint in the experimental group and the control group, and the Pearson correlation coefficient among endpoints from these trials were calculated. Considering the practicality of vaccine clinical trials, the non-inferiority limit was set to 0.693 (ln(1/2)). The mean μ_Ci was set to 3.5, and the variance σ²_Ci was 1. Two-sided test was used and the family-wise type I error rate (α) was 0.05. According to the number of endpoints (K), the correlation coefficient (ρ), the sample size ratio (c_ij), and the type II error rate (β), there were 96 (4*4*3*2) ways of setting the parameters in this simulation study. The numbers of endpoints were taken as 2, 5, 10, and 20, corresponding to the number of actual combinations or multivalent vaccine immunogenic outcomes. The correlation coefficient between 2 different endpoints varied from 0, 0.3, and 0.6 to 0.9 as no correlation, low correlation, moderate correlation or high correlation, respectively. We assumed that c_ij = 1, 1.5, or 3. c_ij = 1 represents the case where the power for each endpoint is equal. When c_ij is greater than 1, the effect size of one endpoint out of K is less than that of the other endpoints. When c_ij = 1, μ_Ti = μ_Ci, diff_i = 0. When c_ij = 1.5, for i = 1,…, K-1, μ_Ti = μ_Ci, diff_i = 0; for i = K, μ_TK = 3.37, μ_CK = 3.5, diff_K = -0.127. When c_ij = 3, for i = 1… K-1, μ_Ti = μ_Ci, diff_i = 0; for i = K, μ_TK = 3.21, μ_CK = 3.5, diff_K = -0.293. The family-wise type II error rate β values were set to 0.1 and 0.2.

We performed a simulation study for all of the 96 combinations of parameters described above and plotted the behavior of the overall power of different settings with increasing sample size. The “Incrementing”, “Bonferroni correction”, and the “Inflation factors” methods were used to calculate the sample size for each of these 96 cases. The results were compared with the simulation results. A deviation of more than 5% might suggest that the method is inappropriate for this situation.

Results

Behavior of overall power

We illustrated the behavior of overall power for 2, 5, 10, and 20 co-primary endpoints. Parameter c_ij was set to 1, 1.5, and 3. The correlation coefficient was set to 0, 0.3, 0.6, and 0.9. Two-sided tests were used, and the significance level α was 0.05. Parameter μ_Ci was set to 3.5, and the variance σ²_Ci was equal to 1. Figure 1 shows the overall power per group obtained by performing 1000 simulations at different sample sizes. Table 1 illustrates the sample size calculated to acquire an overall power of 80% or 90%.

Figure 1. The behavior of overall power with the increasement of sample size of different settings for c_ij and correlation coefficient (r). The type I error rate was 0.05. The numbers of endpoints were 2, 5, 10 and 20 in each set. The range of the horizontal coordinate axis X were 20–95, 35–110 and 65–140, respectively. The sample size shown in the figure is the required sample size for a single group.

The graph clearly shows that the greater the number of co-primary endpoints, the larger is the sample size required to reach the desired overall power. Figure 1 also shows that the larger the correlation coefficient, the closer are the curves of different co-primary endpoints when c_ij is fixed. For example, when c_ij = 1, to acquire an overall power of 80%, trials having 20 co-primary endpoints required 33 more cases per group than trials with 2 endpoints when r = 0; however, it can be as low as 11 when r = 0.9. Moreover, as the ratio of the sample size between different endpoints increases, the curve of overall power with more co-primary endpoints approaches the curve with less co-primary endpoints. Numerically, for the cases where the endpoints are independent from each other (r = 0), 33,27 or 11 more subjects were required for trials with 20 co-primary endpoints than that for trials with 2 co-primary endpoints to obtain an overall power of 0.8 when c_ij was 1,1.5 or 3, respectively.

Comparison of methods

We used numerical examples to compare the 3 sample size calculation methods mentioned in the previous section. Table 2 shows the sample sizes calculated using the “Incrementing” method, “Bonferroni correction” method, “Inflation factors” method, and sample size obtained by simulation. Table 3 shows the factors considered in the sample size calculation by these methods. Considering the situation where the ratio of the sample size at each endpoint is the same (c_ij = 1), the difference between the sample size calculated by the “Inflation factors” method and the reference estimate obtained via simulation did not exceed ± 7%. The sample size obtained using the “Incrementing” method was similar to that obtained using the “Bonferroni correction” method, which did not exceed ± 9% compared with the simulation results when r ≤ 0.3. In particular, the sample sizes obtained using the “Incrementing” method and the sample size obtained using the “Bonferroni correction” method were above the simulation results when K > 2 and r > 0.3. Next, considering the case where the theoretical sample size ratio between one endpoint and the other endpoints is 1.5 (c_ij = 1.5), the sample size calculated by the “Bonferroni correction” method and the “Inflation factors” method was at least 14% more than the simulation result. Meanwhile, when the correlation between the endpoints was small (≤ 0.3), the sample size calculated by the “Incrementing” method did not exceed ± 6% of the sample size acquired by simulation. As the number of co-primary endpoints and the correlation between endpoints increased, the sample size calculated by the “Incrementing” method gradually exceeded the sample size required calculated by simulation. Finally, when the theoretical sample size ratio between one endpoint and the other endpoints was 3 (c_ij = 3), the results obtained using the “Incrementing” method were within ± 6% of the simulation results. However, the sample sizes calculated using the “Inflation factors” method or the “Bonferroni correction” method were at least 19.15% or 35.79% more than that calculated by simulation in order to achieve an overall power of 0.8.

Overall, all 3 methods were applicable when the effect size of each endpoint was similar and there was no or weak correlation between each endpoint (r ≤ 0.3); when the correlation between endpoints is strong (r ≥ 0.6), the “Inflation factors” method considers the influence of correlation between endpoints. When the effect size of individual endpoints is significantly smaller than that of other endpoints (c_ij ≥ 1.5), the “Incrementing” method allocates more power to this endpoint, and thus, requires a smaller sample size to achieve the desired overall power. Moreover, when the gap of effect size between the endpoints is relatively large (c_ij ≥ 1.5), the impact of the correlation on the overall power is almost negligible. In an intermediate case, the sample size under each method can be calculated separately, and the largest sample size can be considered.

Discussion

This article describes 3 sample size calculation methods that are appropriate for vaccine clinical studies having multiple continuous co-primary endpoints. The inflation of type II error rate in combination vaccine trials with multiple co-primary endpoints has been recognized^7,28–32. Calculating the sample size required for each endpoint separately and selecting the maximum value as the sample size for the entire trial often underestimate the demands of sample size¹. A robust method is required for sample size calculation of vaccine trials when there are multiple co-primary endpoints. Most methods in the literature are limited to trials with either 2 or 3 co-primary endpoints or superiority trials^12,14,15,17; however, combination vaccines could have as many as 23 co-primary endpoints³³ and are often designed to be non-inferiority trials especially for generic products. Other methods, such as the “Bonferroni correction” method, although easy to implement in practice, tend to generate conservative sample sizes. Therefore, we have examined 3 ways of calculating the sample size for combination vaccine trials. We first introduced c_ij to indicate the ratio of sample sizes required for each endpoint separately in non-inferiority clinical trials. As c_ij could be simplified as$\frac{{\left(dif{f}_j-d\right)}^2}{{\left(dif{f}_i-d\right)}^2}$, when it takes a value of 1, the effect sizes are expected to be similar across endpoints and the “Inflation factors” method is relatively straightforward. When c_ij was set to 1.5 or 3, the sample size calculated by the “Incrementing” method was closer to the simulation result than that calculated by the other methods. In addition, it is noteworthy that the simulation study was based on the assumption that the effect size of a single endpoint is less than that of the other endpoints, which is a theoretical situation. In fact, the immunogenicity results among endpoints will be moderate; therefore, the results of other approaches will be less conservative.

Based on the assumption that the effect sizes among the endpoints are roughly equal, numerous studies have been conducted; however, our analysis of the actual data found that there is a gap in immunogenicity results among the serotypes or components of combination vaccines. c_ij is not equal to 1 at most times. Supplementary Table 1 shows that the c_ij between strain A1 and strain A3 stimulated to 5.80 in the influenza A subunit vaccine trial; the c_ij between serotype M9V and M18C in a 23-valent pneumococcal polysaccharide vaccine trial was 2.44; the maximum value of c_ij among the immunogenicity endpoints of the Adsorption Free Pertussis Combined Vaccine trial was 2.3. For such trials, the “Incrementing” method can assign more statistical power to the endpoints having a larger effect size, thereby allowing the endpoint with a smaller effect size to have a lower power to achieve the desired overall power and significantly reduce the sample size. Although in the simulation study, we analyzed the behavior of sample size when the correlation among endpoints reaches 0.9, it is a rare situation in actual clinical trials. Correlation coefficients between each serotypes or components of combination vaccines do not exceed 0.6 according to Supplementary Table 2. Due to the presence of cross-immunization³⁴, the body could develop immunity against another pathogen when vaccinated against a particular pathogen, resulting in a positive correlation among immunogenicity results. However, some bacteria or viruses do not have the cross-protection effect owing to the difference in evoking mechanism. Besides, there are many fluctuations during clinical trial conduct, and often the assumptions of treatment response are off the mark. In some cases, the correlation coefficient between the endpoints is unknown. The sample size calculation is based on the assumption that the endpoints are independent from each other. Considering the relatively small impact, for the trials with 2 or 5 co-primary endpoints, the influence of relevance on research results is negligible. However, when there are as many as 10 or more endpoints, a high correlation could cause the sample size calculated by the “Incrementing” method to slightly exceed the sample size required.

This study focused on the influence of the number of endpoints, the difference in immunological effect between the experimental group and the control group, and the correlation among the endpoints on the overall power for vaccine trials with multiple co-primary endpoints. For conservative considerations, the endpoints of each serotype of the control vaccine had the same level of immune response, and therefore, the non-inferiority margins for the different endpoints were the same. When the immune response levels of the various endpoints of the control group are different, we believe that different non-inferiority thresholds should be selected for different species. Although this does not affect the application of the 3 methods described in this paper, it affects the calculation of the parameter c_ij. For example, if we calculate the non-inferiority margin based on the mean difference³⁵, the non-inferiority margin for effect retention is δ = -ε(C-P), where 1-ε is the desirable proportion of the control effect to be retained, C and P refer to the parameters of benefit for the control vaccine and the placebo, respectively. Thus, $c_{ij}=\frac{{\left(dif{f}_j-\delta \right)}^2}{{\left(dif{f}_i-\delta \right)}^2}$. The calculation of cij needs to fully consider the difference in immune response between the experimental group and the control group and the non-inferiority margin.

In this paper, 3 convenient methods were compared. Generally, all 3 methods were acceptable when the effect sizes among all endpoints were similar. The “Inflation factors” method was the only method that could take correlation into consideration. When there is a valid difference in effect sizes among endpoints, the “Inflation factors” method and the “Bonferroni correction” method would be conservative. Our work has been restricted to continuous endpoints in a non-inferiority clinical trial. However, this work serves as the basis for designing other types of randomized clinical trials with multiple co-primary endpoints, including superiority clinical trials. To focus on the effects of correlation on different methods, we only changed the log-transformed immunogenicity response for one endpoint in different settings; the picture on the ground is far more complicated. Moreover, the sample size calculation methods for multiple dichotomous co-primary endpoints still need further discussion since positive rates are often used as the primary endpoint in vaccine trials.

Disclosure of potential conflicts of interest

No potential conflict of interest was reported by the authors.

Supplemental data

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by “the Fundamental Research Funds for the Central Universities” and “Postgraduate Research＆Practice Innovation Program of Jiangsu Province” under Grant KYCX17_0185;the Fundamental Research Funds for the Central Universities and Postgraduate Research * Practice Innovation Program of Jiangsu Province [KYCX17_0185];