In clinical trials with repeated measurements, comparing treatments based on the time-averaged difference (TAD), defined as the difference in the average of longitudinally measured responses between treatment groups, is often considered a meaningful metric for the treatment effect. Overall and Doyle (1994) suggested that TAD can provide a more powerful evaluation of treatment efficacy than the rate of changes over time when the treatment effect has rapid onset and repeated measurements are obtained across an extended period after the maximum effect has been achieved. Many sample size formulas have been developed for the inference of TAD between two treatment groups (Overall and Doyle, 1994; Diggle et al., 2013; Zhang and Ahn, 2012; Lou et al., 2017a). However, sample size calculation for the comparison of time-averaged responses among multiple treatment groups has received less attention in the literature. Randomized trials with multiple treatment arms are widely used in practice (Parmar et al., 2014). They increase the chance of finding an effective treatment as well as reduce the cost and time requirement of clinical trials by testing more treatments simultaneously. For multi-arm trials with continuous outcomes, Zhang and Ahn (2013) presented the sample size formula to compare the time-averaged responses based on the generalized estimating equation (GEE) method (Liang and Zeger, 1986). This sample size formula took into account arbitrary missing patterns, various correlation structures, and unbalanced randomization. Lou et al. (2017b) further extended the sample size approach to multi-arm trials with repeatedly measured count outcomes.

In this paper, we investigate sample size calculation for the comparison of time-averaged responses among K ≥ 3 groups where a binary outcome is repeatedly measured over the study period. The paper is arranged as follows. In Section 2, we briefly review the GEE method for the inference of TAD in multi-arm trials with a repeatedly measured binary outcome. In Section 3, we derive a closed-form sample size formula for the assessment of TAD among treatment groups over time using the GEE method. We demonstrate that this formula is flexible to account for various missing patterns and correlation structures. In Section 4, we conduct simulation studies to assess the performance of the proposed sample size formula under various practical settings. In Section 5, we illustrate the proposed method using a clinical trial example. Section 6 concludes with the discussion.

Suppose in a clinical trial a total of n subjects are enrolled and randomly assigned to one of K treatment groups. Each subject is scheduled to obtain J measurements over the study period. Let Y_{ki j} be the binary response measurement obtained at time t_j( j = 1, …, J) from subject i(i = 1, …, n_k) of the kth treatment group (k = 1, …, K), where n_k denotes the number of subjects assigned to the kth treatment group. We use r_k = n_k/n to denote the proportion of subjects assigned to the kth treatment group. To make inference about the TAD among the K treatment groups, we model Y_{ki j} with the following logistic model:

We have E(Y_{ki j}) = p_k = e^b_k/(1 + e^b_k). That is, b = (b₁, …, b_K)′ represents the time-averaged response on the log-odds scale for the K treatment groups. As a binary variable, it is obvious that Var(Y_{ki j}) = p_k(1− p_k). Furthermore, we model the within-subject correlation among measurements obtained from the same subject by corr(Y_{ki j}, Y_{ki j′}) = ρ_{j j′} with ρ_{j j} = 1. We assume Y_{ki j}’s to be independent across subjects.

The GEE estimator of b, denoted by b̂ = (b̂₁, …,b̂_K)′, can be obtained from the following equation:

which is derived based on an independent working correlation structure. Liang and Zeger (1986) showed that as n → ∞, $\sqrt{n}(\widehat{\mathit{b}}-\mathit{b})$ approximately has a normal distribution with mean vector 0 and variance matrix V. We can consistently estimate V by ${\mathit{V}}_n=\mathit{W}{\mathit{A}}_n^{-1}(\widehat{\mathit{b}}){\mathbf{\Sigma }}_n{\mathit{A}}_n^{-1}(\widehat{\mathit{b}})\mathit{W}$, where W is a diagonal matrix with diagonal elements being ($1/\sqrt{r_1},\dots ,1/\sqrt{r_K}$),

and

Here ε̂_{ki j} = y_{ki j} − e^b̂_k/(1 − e^b̂_k) denotes the residual.

To compare the time-averaged responses among K treatment groups, the null hypotheses of interest is H₀ : b₁ = · · · = b_K. The test statistic is

where C = (c₁, …, c_K)′ is a vector denoting a contrast of treatment effects with ${\sum }_{k=1}^K{c}_k=0$. For example, one reasonable specification would be C = (−1, 1/(K − 1), …, 1/(K − 1))′. The null hypothesis is rejected if |Z| > z_1−α/2, where z_1−α/2 is the 100(1−α/2)th percentile of the standard normal distribution.

As n → ∞, let A and ∑ denote the limits of the A_n and ∑_n, respectively. Then V_n converges to V = WA⁻¹(b̂)∑A⁻¹(b̂)W. Given true treatment effects θ = (θ₁, …, θ_K), the required sample size to reject the null hypothesis with power 1 − γ at significance level α is

Missing data are frequently encountered in clinical trials with repeated measurements. We now show that a closed-form extension of (3.1) can be obtained to account for missing data. Let Δ_{ki j} be the missing indicator, which takes value 0/1 for a missed/observed outcome measurement. Then the generalized expressions of A_n and ∑_n that accommodate missing data are

and

We assume that the missing probabilities only depend on time. In addition, we define δ_j = E(Δ_{ki j}) to be the probability of a subject having a measurement at time t_j and δ_{j j′} = E(Δ_{ki j}Δ_{ki j′}) to be the probability of a subject having measurements simultaneously at t_j and t_j′. Note that δ_{j j} = δ_j. We can use probabilities δ_j and δ_{j j′} to describe various types of missing patterns. Then, as n → ∞, we have

and

The general variance V can be expressed as

Plugging V and C = (−1, 1/(K − 1), …, 1(K − 1))′ into Equation (3.1), the generalized sample size formula that accommodates various correlation structures and missing data patterns is

To assess the performance of the proposed sample size method, we conduct simulation under different parameter settings. Suppose subjects are randomly assigned to one of four treatment groups (K = 4). Each subject is planned to obtain T = 6 measurements at scheduled times (t_j = j − 1 for j = 1, …, 6). We investigate three missing patterns: independent missing (IM), monotone missing (MM), and mixed missing (MIX). Under IM, missing measurements occur independently over time with δ_{j j′} = δ_jδ_j′ for j ≠ j′. Under MM, a subject missing a measurement at time t_j will miss all the following measurements, such that δ_{j j′} = δ_j′ for j ≤ j′. IM represents scenarios where subjects miss clinical visits due to random reasons. MM represents scenarios where subjects randomly drop out during the study period. In real clinical trials, it is likely that subject dropout and missing visits both occur, which implies a mixture of IM and MM, denoted by MIX. Let (${\delta }_1^{(\text{IM})},\dots ,{\delta }_J^{(\text{IM})}$) and (${\delta }_1^{(\text{MM})},\dots ,{\delta }_J^{(\text{MM})}$) be the marginal probabilities of obtaining a measurement under the IM and MM patterns, respectively. Moreover, let ${\delta }_{j{j}^{\prime }}^{(\text{IM})}$ and ${\delta }_{j{j}^{\prime }}^{(\text{MM})}$ be the corresponding joint probabilities under these two patterns. Suppose the proportion of subjects who follow IM and MM patterns are w and 1 − w, then for mixed missing, we have

In simulation studies, we set w = 0.5 for the MIX pattern. For marginal observation probabilities, we assume δ^(IM) = δ^(MM) = δ = (δ₁, …, δ_J), where four sets of values are explored:

These settings correspond to four scenarios. We assume no missing data at t₁. δ₁ indicates complete data throughout the study. δ₂ indicates that the missing probabilities have a linear trend with 5% increase in dropout rate at each subsequent time point. δ₃ indicates an accelerating trend missing probabilities toward the end of study. δ₄ indicates a trend opposite to that of δ₃. Note that δ₂ − δ₄ have the same dropout rate (25%) at the end of study.

The simulation study also explores two correlation structures: compound symmetry (CS, ρ_{j j′} = ρ for j ≠ j′) and auto-regressive (AR(1), ρ_{j j′} = ρ^{|t_j−t_j′|}) with ρ = 0.3, 0.5. We set the Type I error and power at α = 0.05 and 1 − γ = 0.8, respectively. For simplicity, we assume the balanced design with r₁ = · · · = r_K = 1/K. However, the proposed sample size formula is applicable to any unbalance design. We consider two alternative hypotheses. The first alternative hypothesis is that the first group is control group. The others are treatment groups with same treatment effect. Specifically, H₁ : θ₁ = 0, θ₂ = θ₃ = θ₄ = 0.5. The second alternative hypothesis is that the treatment groups are ordered by treatment effects. Specifically, H₁ : θ₁ = 0, θ_k = θ₁ +(k − 1) Δ₀ for k = 2, 3, 4 and Δ₀ = 0.25. For each combination of design parameters (missing pattern, marginal observation probability δ, correlation structure, correlation ρ, alternative hypothesis H₁), the simulation is conducted as:

• Calculate the required sample size (n) based on Equation (3.2).

• Generate n samples under null hypothesis and alternative hypothesis separately. Each sample contains n multivariate Bernoulli random variable Y_ki = (Y_ki1, …, Y_kiJ)′ with correlation parameter ρ_{j j′} under the assumed correlation structure. The correlated binary vectors are generated by the method of Emrich and Piedmonte (1991).

• Create incomplete data sets. Generate missing indicators based on the specified missing pattern and the marginal observation probabilities.

• Calculate b̂, A_n, ∑_n, and W, and the test statistic Z based on Equation (2.1).

• Repeat Steps 2–4 for L = 10,000 times. The empirical Type I error and empirical power are estimated by ${\sum }_{l=1}^{L}I\hspace{0.17em}(|Z|>z_{1-\alpha /2})/L$ under the null hypothesis and alternative hypothesis, respectively.

Tables 1 and 2 summarize the required sample sizes and their corresponding empirical Type I errors and empirical powers under different combinations of simulation parameters. Table 1 is obtained under alternative hypotheses H₁ : θ₁ = 0, θ₂ = θ₃ = θ₄ = 0.5, and Table 2 under H₁ : θ₁ = 0, θ₂ = 0.25, θ₃ = 0.5, θ₄ = 0.75. We have several observations: (1) The empirical powers and Type I errors are close to the nominal levels, which are 0.8 and 0.05, respectively. Therefore, the proposed sample size formula has a good performance over a wide range of design parameter settings; (2) The sample sizes under the AR(1) are always smaller than those under the CS correlation structure, with other design parameters being the same; (3) The sample size increases with the correlation (ρ), which is obvious from the sample size formula (3.2); (4) Given the same marginal observation probabilities, the order of the required sample sizes under different missing pattern is MM > MIX > IM. Under MM missing data tend to concentrate in successive observations from a few subjects, while under IM missing data are randomly distributed across subjects over the study period. The MM missing pattern leads to greater information loss (hence greater sample size requirement) than IM, with MIX lies in between. (5) Despite the same dropout rate at the end of study, the order of required sample sizes is δ₄ > δ₂ > δ₃, because the overall proportion of missing values is the greatest under δ₄.

PASS sample size software manual (Pass14, 2015) illustrated the sample size estimation to test proportions in a repeated measurement design between two treatment groups. Here, we show sample size estimation to test proportions among three treatment groups. An investigator wants to design a study that compares the efficacy of a prophylactic treatment for the common cold with two active drugs and a placebo. The null hypothesis is that there is no difference in the proportion of patients who get sick among three treatment groups. Patients will be randomly assigned to one of three treatment groups with an equal probability, and followed monthly from September to April (beginning in October, hence J = 7) to determine the patient’s disease status (present or absent). The study investigated if there is an overall difference in the proportion of patients who get sick among three treatment groups. A baseline of 60% disease rate for the common cold is estimated based on previous studies. It is expected that the disease rate will continue to be 60% in the placebo group. Suppose that a clinically meaningful difference in efficacy is 30% relative reduction in disease rate in two active medication groups, which corresponds to a disease rate of 42%. The hypothesis of interest is then H₀ : b₁ = b₂ = b₃ versus H₁ : b₁ = logit(60%) = 0.4055, b₂ = b₃ = logit(42%) = −0.3228. We would like to calculate the sample size that can detect the difference in treatment effect with Type I error α = 0.05 and power 1 − γ = 0.8 under a balanced design. We set the measurement times at t_j = j−1( j = 1, …, 7). The observation probabilities are assumed to follow a linear trend with a 30% dropout at the end of study, δ = (1, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70). Under the AR(1) correlation structure with ρ = 0.5, the sample sizes required under the IM, MM, and MIX (assuming a balanced mixture of IM and MM) patterns are 104, 110, and 107, respectively. However, the required sample sizes under the IM, MM, and MIX are 165, 175, and 170 under the CS correlation structure.

In this study, we derived a sample size formula to compare the time-averaged responses of repeated binary outcomes among K groups. This proposed sample size formula can accommodate arbitrary correlation structures, missing patterns, marginal observation probabilities, and unbalanced experimental designs. We develop a sample size formula based on the GEE method because: (1) It has been widely used to analyze data from clinical trials with longitudinal/repeated observations; (2) It is robust to the misspecification of correlation structure (Liang and Zeger, 1986; Jung and Ahn, 2003); (3) It is flexible to accommodate missing data (Zeger et al., 1988). Our simulation studies show that the empirical powers and the empirical Type I errors are very close to the nominal levels under a wide range of design parameters. When we compare the time-averaged responses among K groups, a larger correlation is always associated with a larger sample size, which is obvious in Equation (3.2) and discussed in Lou et al. (2017a).