Chapter 18

Sample Size Calculation for Observational Studies

Sin-Ho Jung
Taiyeong Lee
Elizabeth DeLong

Abstract

18.1 Introduction

18.2 Continuous Variables

18.3 Binary Variables

18.4 Two-Sample Log-Rank Test for Survival Data

18.5 Two-Sample Longitudinal Data

18.6 Discussion

Appendix: Asymptotic Distribution of Wilcoxon Rank Sum Test under H_a

References

 

Abstract

For a number of reasons, observational studies are currently being used to provide evidence to support medical decisions. Because these studies do not carry the same level of credibility conferred to randomized trials, appropriate attention to statistical issues is paramount. In particular, considerations of event rate, sample size, and power frequently occur. This chapter provides mechanisms for dealing with these considerations under a variety of scenarios. Additionally, the methods presented here are applicable more generally to clinical trials as well.

 

18.1 Introduction

The increasing availability of computerized clinical patient data has stimulated greater interest in research on observational data. Individual medical practices and hospitals are now attempting to draw inferences on the basis of data gathered on their own patients. Additionally, the proliferation of large registry databases that harvest clinical data from hundreds of sites has created opportunities to study questions regarding rare diseases and rare events. However, even studies using these large databases, when reduced to the appropriate study populations and endpoints, encounter issues of sample size and power.

Observational studies can be designed to address a variety of objectives. If the primary analysis of the study is simply to estimate some parameter, then one can use standard sample size calculations for requiring a confidence interval of at most a certain width. If the objective involves the comparison of cohorts, then in many cases the statistical analysis will involve some form of selection bias adjustment due to the non-randomized nature of observational research. Chapters 1 through 11 cover the issues of selection bias adjustment and standard analytical methods for comparing cohorts in these settings, including propensity score regression, stratification, matching, doubly robust adjustment, instrumental variables, local control, and various longitudinal analyses. Regarding the computation of sample sizes in study planning in these scenarios, such complex designs and analyses may require custom sample size and power estimate calculations through simulations. In some cases, the primary hypothesis may be simplified for the purpose of estimating power and sample size, and generalized SAS macros can perform the needed calculations.

In this chapter, we present a series of SAS macros for computing sample sizes for studies with various types of outcomes and comparisons (all for comparisons between two groups):

• t-test (Section 18.2.1)

• Wilcoxon Rank Sum test (Section 18.2.2)

• Two-Sample Tests on Binary outcomes (Section 18.3.1)

• Weighted Mantel-Haenszel test (Section 18.3.2)

• Log-Rank test for survival data (Section 18.4)

• Longitudinal Data—Continuous outcomes (Section 18.5.3.3)

• Longitudinal Data—Binary outcomes (Section 18.5.3.3)

The sample size macro for the Weighted Mantel-Haenszel test (provided in Section 18.3.1) is directly applicable for observational studies for which the analysis is based on the commonly utilized propensity score stratification approach. This propensity score stratification analysis methodology is described in detail in Chapter 2. In brief, propensity score stratification involves estimating a propensity score for each study participant, stratifying by this score, estimating the difference in outcomes between the two cohorts within each of the strata, and then combining the estimated cohort differences over strata.

Other sections in this chapter present sample size macros based on analysis methods without specific selection bias adjustments. However, these can be very useful for observational study planning when simple cohort comparisons are planned (for example, such as when causal inference is not the objective) or as a very easy-to-use initial calculation when other methods are being utilized such as propensity score matching. In this example, one would supplement the sample sizes computed in these macros with an estimate of the number of patients excluded from the analysis due to a lack of propensity score matched patients in the other cohort. In addition, most of the sample size calculation methods in this chapter can be used for randomized clinical trials as well.

 

18.2 Continuous Variables

In this section, we consider two two-sample tests for continuous variables: t-test and Wilcoxon rank test.

 

18.2.1 Two-Sample -test

Suppose that, for group k(=1,2), X_k1,..., X_k,nk are independent and identically distributed (IID) normal random variables with mean μ_k and variance σ². We want to test H₀ : μ₁ = μ₂ vs. $H_a:{\mu }_1\ne {\mu }_2$. ${\overline{X}}_k=n_k^{-1}{\sum }_{i=1}^{n_k}{X}_{ki}$ Let denote the sample mean for group k, and

$s_p^2=\frac{\sum _{i=1}^{n_1}{(X_{1i}-{\overline{X}}_1)}^2+\sum _{i=1}^{n_2}{(X_{2i}-{\overline{X}}_2)}^2}{n_1+n_2-2}$

denote the pooled sample variance. Then, under H₀,

$T=\frac{{\overline{X}}_1-{\overline{X}}_2}{s_p\sqrt{n_1^{-1}+n_2^{-1}}}$

follows the t-distribution with n₁ + n₂ − 2 degrees of freedom. Hence, given type I error probability α, we reject H₀ if $\left|T\right|>t_{n_1+n_2-2,1-\alpha /2}$, where t_v,γ is the 100γ-th percentile of the t-distribution with v degrees of freedom.

If n₁ and n₂ are large, then we do not require the normal distribution assumption, and the critical value from the t-test, $t_{n_1+n_2-2,1-\alpha /2}$ , can be approximated by the 100γ-th percentile of the standard normal distribution, z_1−α/2. How large the sample sizes should be for the large approximation depends on how close the distribution. of the observations for each group is to a normal distribution. To improve the normality, we often apply a transform to the raw data, such as a log-transformation for positive variables (Carroll and Ruppert, 1988). We derive the sample size formula based on the large sample approximation.

We consider a specific alternative hypothesis H_a : |μ₁ − μ₂|= ∆σ. Note that ∆ = σ^-1 (μ₁ − μ₂) denotes a standardized effect size under the alternative hypothesis. Without loss of generality, we assume that μ₁ > μ₂ under H_a. Then, under H_a,

$T\approx \frac{{\overline{X}}_1-{\overline{X}}_2-\sigma \Delta }{\sigma \sqrt{n_1^{-1}+n_2^{-1}}}+\frac{\Delta }{\sqrt{n_1^{-1}+n_2^{-1}}}$ $=Z+\Delta \sqrt{n{r}_1{r}_2},$

where Z : N(0, 1), r_k = n_k/n, and n = n₁ + n₂. For a given n, the power is calculated as

$1-\beta =\overline{\varphi }(z_{1-\alpha /2}-\Delta \sqrt{n{r}_1{r}_2}),$                                    (1)

where $\overline{\varphi }(z)=1-\varphi (z)$ and $\varphi (z)=P(Z\le z)$ is the cumulative distribution function of the standard normal distribution, N(0, 1). By solving (1) with respect to n, we obtain the required total sample size as

$n=\frac{{(z_{1-\alpha /2}+z_{1-\beta })}^2}{r_1{r}_2{\Delta }^2},$

or n_k = n х r_k for group k. For one-sided testing, replace α/2 with α. This sample size is based on a large sample approximation to a t-distribution, so that this formula underestimates the sample size for small sample sizes. If the final sample size by the formula based on normal approximation is smaller than 30, we may consider increasing it by about 10%.

In summary, the sample size calculation for a two-sample t-test is conducted as follows:

• In addition to type I error probability α and power 1 − β, specify

(a) $\Delta ={\sigma }^{-1}|{\mu }_1-{\mu }_2|$, standardized effect size

(b) the prevalence of group k $(r_1+r_2=1)$

• Calculate the sample size

$n=\frac{{(z_{1-\alpha /2}+z_{1-\beta })}^2}{r_1{r}_2{\Delta }^2}\mathrm{.}$

Example 18.1

Suppose that we want to detect a difference of 0.5σ in the population means with two-sided α= 0.05 and power 1 − β = 0.9. Then $z_{1-\alpha /2}=1.96$ and $z_{1-\beta }=1.282$. Assuming an equal proportion of two groups in the population (that is $r_1=r_2=0.5$), we obtain the required sample size . PROC POWER in SAS/STAT can directly handle this scenario, as shown in the following code:

Program 18.1 SAS Code for Two-Sample t-Test

proc power;
     twosamplemeans
     nfractional
     meandiff = 0.5              /* mean difference                */
     stddev  = 1                 /* sigma                          */
     groupweights = (0.5 0.5)    /* r1 r2                          */
     sides = 2                   /* 1: one sided 2: 2: two-sided   */
     power =0.9                  /* power                          */
     alpha =0.05                 /* alpha                          */
     ntotal=.;                   /* Sample Size                    */
run;

 

Output from Program 18.1

 

18.2.2 Wilcoxon Rank Sum Test

It is known that -tests are sensitive to outliers. The Wilcoxon rank sum test (WRST) has been widely used as a robust test. We assume that for group $k(=1,2)$, $X_{k1},...,X_{k,n_k}$ are IID random variables from a distribution with cumulative distribution function $F(x-{\theta }_k)= P (X_{ki}\le x)$. For $\Delta ={\theta }_1-{\theta }_2$, we want to test $H_0:\Delta =0$vs. $H_a:\Delta \ne 0$. Mann and Whitney (1947) propose to use $W=(n_1{n}_2{)}^{-1}{\sum }_{i=1}^{n_1}{\sum }_{j=1}^{n_2}I(X_{1i}>X_{2j})$ for testing the hypothesis. The expected value of W is $P (X_{1i}>X_{2j})$, the probability that a randomly chosen measurement from group 1 is greater than a randomly chosen measurement from group 2. Thus, the test statistic will be close to 1/2 if H₀ is true, and closer to 0 or 1 if H_a is true.

For large n(=n₁ + n₂), we reject if the absolute value of

$T=\frac{W-{\nu }_0}{{\sigma }_0}$

is larger than z_1−α/2, where v₀ = 1/2 and ${\sigma }_0^2=(n+1)/(12n_1{n}_2)$ are the mean and variance of W under H₀, respectively.

Let $f(x)=\partial F/\partial x$ denote the probability density function of F(x). The appendix shows that, under H_a, W has mean

${\nu }_a={\int }_{-\infty }^{\infty }F(x+\Delta )f(x)dx$

and variance

${\sigma }_a^2=\frac{{\sigma }_1^2}{n_1}+\frac{{\sigma }_2^2}{n_2},$

where

${\sigma }_1^2={\int }_{-\infty }^{\infty }{F}^2(x+\Delta )f(x)dx-{\{{\int }_{-\infty }^{\infty }F(x+\Delta )f(x)dx\}}^2,$

and

${\sigma }_2^2={\int }_{-\infty }^{\infty }{F}^2(x-\Delta )f(x)dx-{\{{\int }_{-\infty }^{\infty }F(x-\Delta )f(x)dx\}}^2\mathrm{.}$

Then,

$T=Z\frac{{\sigma }_a}{{\sigma }_0}+\frac{{\nu }_a-{\nu }_0}{{\sigma }_0},$

where $Z=(W-{\nu }_a)/{\sigma }_a$ is an asymptotically random variable. Given power 1 − β,

$1-\beta = P (|T|>z_{1-\alpha /2}|H_a)=P(Z\frac{{\sigma }_a}{{\sigma }_0}+\frac{|{\nu }_a-{\nu }_0|}{{\sigma }_0}>z_{1-\alpha /2}).$

$=\overline{\varphi }(\frac{{\sigma }_0}{{\sigma }_a}(z_{1-\alpha /2}-\frac{|{\nu }_a-{\nu }_0|}{{\sigma }_0})),$

where $\overline{\varphi }(z)= P (Z>z)$. By replacing n_k with nr_k in ${\sigma }_0^2$ and ${\sigma }_a^2$, and solving with respect to n, we obtain the required sample size

$n=\frac{1}{12r_1{r}_2}{\left\{\frac{z_{1-\alpha /2}+z_{1-\beta }\sqrt{12(r_2{\sigma }_1^2+r_1{\sigma }_2^2)}}{{\nu }_a-1/2}\right\}}^2\mathrm{.}$

Although the testing does not require the specification of the distribution function F(x), a sample size calculation does. The alternative hypothesis is specified by two distributions with a location shift. The assumption of a location shift model is not required in order to carry out the rank sum test, but it is needed if we want to make inferences about medians or means. It is, of course, needed in the sample size calculation as presented here.

While a t-test requires a finite variance of the observations, WRST does not. For example, if has a Cauchy distribution, we can use WRST, but not a t-test.

A sample size calculation for WRST assuming normal distributions can be described as follows:

• In addition to type I error probability α and power 1 − β, specify

(a) $\Delta ={\sigma }^{-1}({\mu }_1-{\mu }_2)$, standardized effect size

(b) r_k = the prevalence of group k $(r_1+r_2=1)$

• Calculate

${\nu }_a={\int }_{-\infty }^{\infty }\varphi (x+\Delta )\varphi (x)dx$

${\sigma }_1^2={\int }_{-\infty }^{\infty }{\varphi }^2(x+\Delta )\varphi (x)dx-{\{{\int }_{-\infty }^{\infty }\varphi (x+\Delta )\varphi (x)dx\}}^2$

and

${\sigma }_2^2={\int }_{-\infty }^{\infty }{\varphi }^2(x-\Delta )\varphi (x)dx-{\{{\int }_{-\infty }^{\infty }\varphi (x-\Delta )\varphi (x)dx\}}^2,$

where $\varphi (x)$ and $\varphi (x)$ are the probability density function and the cumulative distribution function of , respectively.

• Calculate the sample size

$n=\frac{1}{12r_1{r}_2}{\left\{\frac{z_{1-\alpha /2}+z_{1-\beta }\sqrt{12(r_2{\sigma }_1^2+r_1{\sigma }_2^2)}}{{\nu }_a-1/2}\right\}}^2\mathrm{.}$

Example 18.2

Suppose that we want to detect a difference of 0.5σ in the population means with two-sided α = 0.05 and power 1 − β = 0.9. Then $z_{1-\alpha /2}=1.96$ and $z_{1-\beta }=1.282$. Assuming an equal proportion of two groups in the population (that is, $r_1=r_2=0.5$) , we obtain the required sample size n = 177. Note that the two-sample t-test requires a slightly smaller n = 171, but, in case there are outliers in the final data, the Wilcoxon rank test may be more powerful.

Program 18.2 SAS Code for Wilcoxon Rank Sum Test

%macro  WilcoxonRankSumTest( delta=,    /* Standardized effect size     */
                             r1 =  ,    /* prevalence of group 1        */
                             r2 =  ,    /* prevalence of group 2        */
                             sides=,    /* number of test sides: 1 or 2 */
                             alpha=,    /* alpha                        */
                             power=     /* power                        */
                             );
proc iml;
  start CDFPDFDELTA(t);
      m=t+δ c=cdf('normal',m,0,1); p=pdf('normal',t,0,1); v=c*p;
return(v);
  finish;
  start CDF2PDFDELTA(t);
      m=t+δ c=cdf('normal',m,0,1)**2; p=pdf('normal',t,0,1); v=c*p;
return(v);
  finish;

  start CDFPDF_DELTA(t);
      m=t-δ c=cdf('normal',m,0,1); p=pdf('normal',t,0,1); v=c*p;
return(v);
  finish;
  start CDF2PDF_DELTA(t);
      m=t-δ c=cdf('normal',m,0,1)**2; p=pdf('normal',t,0,1); v=c*p;
return(v);
  finish;
  start s1(t);
      m=t+δ c=cdf('normal',m,0,1); p=pdf('normal',t,0,1); v=c*p;
return(v);
  finish;

  interval = .M || .P;
  call quad(nu1,"CDFPDFDELTA",   interval);
  call quad(s1, "CDF2PDFDELTA",  interval);
  call quad(nu2,"CDFPDF_DELTA",  interval);
  call quad(s2, "CDF2PDF_DELTA", interval);

  sigma1 = s1-nu1**2; sigma2 = s2-nu2**2;
  p1=1-&alpha/&sides; z_alpha = probit(p1);
  p2=&power; z_beta  = probit(p2);
  n =
(1/(12*&r1*&r2))*((z_alpha+z_beta*sqrt(12*(&r2*sigma1+&r1*sigma2)))/(nu1-
0.5))**2;
  delta=δr1 =&r1; r2=&r2;sides =&sides;alpha=αpower=&power;
  print  'Sample Size';
  print 'Wilcoxon Rank Sum Test for Mean Difference';
  print delta r1 r2 sides alpha power;
  print n;
quit;
run;
%mend WilcoxonRankSumTest;

/*---------------------- Run the macro for Example 2. -------------------*/

%WilcoxonRankSumTest(delta=0.5, r1 = 0.5, r2 = 0.5, sides=2, alpha=0.05, power=0.9);

 

Output from Program 18.2

 

18.3 Binary Variables

In this section, we investigate two-sample tests for binary outcome variables with or without stratification. We use a weighted Mantel-Haenszel test for stratified analysis to adjust for covariates that are unbalanced between two groups.

 

18.3.1 Two-Sample Test on a Binary Outcome

Suppose that, for group k(=1, 2), X_k denotes the number of responders from n_k independent subjects. If group k has a response probability p_k, then X_k is a binomial distribution with n_k independent trials and response probability p_k. We want to test $H_0:p_1=p_2$ vs. $H_a:p_1\ne p_2$. Let ${\widehat{p}}_k=X_k/n_k$ and $\widehat{p}=(X_1+X_2)/(n_1+n_2)$ denote the sample proportion for group k and the pooled data, respectively. Then, for large n₁ and n₂ (say, n₁,n₂ > 30) ,

$T=\frac{{\widehat{p}}_1-{\widehat{p}}_2}{\sqrt{\widehat{p}\widehat{q}(n_1^{-1}+n_2^{-1})}}$

follows the standard normal distribution under H₀, where $\widehat{q}=1-\widehat{p}$. Hence, given type I error probability α, we reject H₀ if $\left|T\right|>z_{1-\alpha /2}$. It is easy to show that T² is identical to the chi-squared test with 1 degree of freedom for a 2 х 2 table.

Let p₁, p₂ denote the response probability under a specific alternative hypothesis for sample size calculation, and $\Delta =p_1-p_2$. Then, for large $n(=n_1+n_2)$, $\widehat{p}$ converges to $\overline{p}=r_1{p}_1+r_2{p}_2$, where $r_k=n_k/n$. Under H_a,

$T\approx \frac{{\widehat{p}}_1-{\widehat{p}}_2-\Delta }{\sqrt{\overline{p}\overline{q}(n_1^{-1}+n_2^{-1})}}+\frac{\Delta }{\sqrt{\overline{p}\overline{q}(n_1^{-1}+n_2^{-1})}}$

$=Z+\Delta \sqrt{\frac{n{r}_1{r}_2}{\overline{p}\overline{q}}},$

where Z : Z(0, 1) and $\overline{q}=1-\overline{p}$. Hence, given n, the power is calculated as

$1-\beta =\overline{\varphi }(z_{1-\alpha /2}-\Delta \sqrt{n{r}_1{r}_2}),$                                       (2)

By solving (2) with respect to n, we obtain the required total sample size as

$n=\frac{{(z_{1-\alpha /2}+z_{1-\beta })}^2}{r_1{r}_2{\Delta }^2},$

In summary, the sample size calculation for two-sample binomial proportions is conducted as follows:

• In addition to type I error probability α and power 1 − β, specify

(a) p₁, p₂ = binomial proportions under H_a

(b) r_k the prevalence of group k (r₁ + r₂ = 1)

• Calculate the sample size

$n=\frac{{(z_{1-\alpha /2}+z_{1-\beta })}^2\overline{p}\overline{q}}{r_1{r}_2{(p_1-p_2)}^2},$

where $\overline{p}=r_1{p}_1+r_2{p}_2$and $\overline{q}=1-\overline{p}$.

Example 18.3

Suppose that we want to detect a difference of p₁ = 0.4 and p₂ = 0.5 with two-sided α = 0.05 and power 1 − β = 0.9. Assuming an equal proportion of two groups in the population (that is, r₁ = r₂ = 0.5), we obtain the required sample size n = 10.38.

Program 18.3 SAS Code for Two-Sample Binomial Proportions

proc power;
     twosamplefreq  test=pchi
     groupproportions =(0.4 0.5) /* binomial proportions to be detected  */
     groupweights = (1 1)        /* (w1, w2):Weight of two groups,(r1,r2)*/
     sides = 2                   /* 1: one sided 2: 2: two-sided         */
     power =0.9                  /* power                                */
     alpha =0.05                 /* alpha                                */
     ntotal=.;                   /* sample size                         */
    run;

 

Output from Program 18.3

 

18.3.2 Weighted Mantel-Haenszel Test

In an observational study, patients are not randomly assigned to treatment groups with equal probability. Instead, the probability of assignment varies from patient to patient, possibly depending on the patient’s baseline covariates. This often results in non-comparable treatment groups due to imbalance of the baseline covariates and consequently invalidates the standard methods commonly employed in data analysis. To overcome this problem, Rosenbaum and Rubin (1983, 1984) developed the propensity score method.

Suppose that there are j strata formed by propensity score stratification. Let n denote the total sample size and n_j the sample size in stratum j ${\sum }_{j=1}^J{n}_j=n$. The data on each subject comprise the response variable x = 1 for response and 0 for no response and j and k for the stratum and treatment group, respectively, to which the subject is assigned ($1\le j\le J;k=1,2$). Frequency data in stratum j can be described as follows:

	Group
Response	1	2	Total
Yes	xj11	xj12	xj1
No	xj21	xj22	xj2
Total	nj1	nj2	nj

Let $O_j=x_{j11}$, $E_j=n_{j1}{x}_{j1}/n_j$, and

$V_j=\frac{n_{j1}{n}_{j2}{x}_{j1}{x}_{j2}}{n_j^2(n_j-1)}\mathrm{.}$

Then, the weighted Mantel-Haenszel (WMH) test is given by

$T=\frac{\sum _{j=1}^J{\widehat{w}}_j(O_j-E_j)}{\sqrt{\sum _{j=1}^J{\widehat{w}}_j^2{V}_j}},$

where the weights ${\widehat{w}}_j$ converge to a constant w_j as $n\to \infty$. The weights are ${\widehat{w}}_j=1$ for the original Mantel-Haenszel test and ${\widehat{w}}_j={\widehat{q}}_j=x_{j2}/n_j$ for the statistic proposed by Gart (1985).

Let a_j = n_j/n denote the allocation proportion for stratum $({\sum }_{j=1}^J{a}_j=1)$, and $b_{jk}=n_{jk}/n_j$ denote the allocation proportion for group k within stratum j $(b_{j1}+b_{j2}=1)$. Let p_jk denote the response probability for group k in stratum j and $q_{jk}=1-p_{jk}$. Under $H_0:p_{j1}=p_{j2},1\le j\le J$, T is approximately N(0, 1). The optimal weights maximizing the power depend on the allocation proportions $\{(a_j,b_{1j},b_{2j}),j=1,...,J\}$ and effect sizes $(p_{j1}-p_{j2},1,...,J)$ under H₁.

In order to calculate the power of WMH, we have to derive the asymptotic distribution of ${\sum }_{j=1}^J{\widehat{w}}_j(O_j-E_j)$ and the limit of ${\sum }_{j=1}^J{\widehat{w}}_j^2{V}_j$ under H₁. We assume that the success probabilities satisfy for under . Note that $(p_{jk},1\le j\le J,j=1,2)$ satisfy $p_{j2}{q}_{j1}/(p_{j1}{q}_{j2})=\varphi$ for $\varphi \ne 1$ under H₁. Note that a constant odds ratio across strata holds if there exists no interaction between the treatment and the propensity score when the binary response is regressed on the treatment indicator and the propensity score using a logistic regression. The following derivations are based on H₁. It can be verified that

   $O_j-E_j=\frac{n_{j1}{n}_{j2}}{n_j}({\hat{p}}_{j1}-{\hat{p}}_{j2})$

   $=\frac{n_{j1}{n}_{j2}}{n_j}({\hat{p}}_{j1}-p_{j1}-{\hat{p}}_{j2}+p_{j2})+\frac{n_{j1}{n}_{j2}}{n_j}(p_{j1}-p_{j2})$

   $=n{a}_j{b}_{j1}{b}_{j2}({\hat{p}}_{j1}-p_{j1}-{\hat{p}}_{j2}+p_{j2})+n{a}_j{b}_{j1}{b}_{j2}(p_{j1}-p_{j2}).$

Thus, under H₁, ${\sum }_{j=1}^J{\hat{w}}_j(O_j-E_j)$ is approximately normal with mean nΔ and variance $n{\sigma }_1^2$, where

   $\Delta =\sum _{j=1}^J{w}_j{a}_j{b}_{j1}{b}_{j2}(p_{j1}-p_{j2})$

   $=(1-\Phi )\sum _{j=1}^J{w}_j{a}_j{b}_{j1}{b}_{j2}\frac{p_{j1}{q}_{j1}}{q_{j1}+\Phi p_{j1}}$

and

   ${\sigma }_1^2=n^{-1}\sum _{j=1}^J{w}_j^2\frac{n_{j1}^2{n}_{j2}^2}{n_j^2}(\frac{p_{j1}{q}_{j1}}{n_{j1}}+\frac{p_{j2}{q}_{j2}}{n_{j2}})$

   $=\sum _{j=1}^J{w}_j^2{a}_j{b}_{j1}{b}_{j2}(b_{j2}{p}_{j1}{q}_{j1}+b_{j1}{p}_{j2}{q}_{j2}).$

Also under H₁, we have

$$\sum _{j=1}^J{w}_j^2{V}_j=n{\sigma }_0^2+o_p(n),$$

where

$${\sigma }_0^2=\sum _{j=1}^J{w}_j^2{a}_j{b}_{j1}{b}_{j2}(b_{j1}{p}_{j1}+b_{j2}{p}_{j2})(b_{j1}{q}_{j1}+b_{j2}{q}_{j2}).$$

Hence, the power of WMH is given as

$$1-\beta = P (|T|>z_{1-\alpha /2}|H_1)$$

$$= P (\frac{{\sigma }_1}{{\sigma }_0}Z+\sqrt{n}\frac{\left|\Delta \right|}{{\sigma }_0}>z_{1-\alpha /2})$$

$$=\overline{\Phi }(\frac{{\sigma }_0}{{\sigma }_1}{z}_{1-\alpha /2}-\sqrt{n}\frac{\left|\Delta \right|}{{\sigma }_1}),$$

where Z is a standard normal random variable and $\overline{\Phi }(z)=P(Z>z)$ . Thus, the sample size required for achieving a desired power of 1 − β can be obtained as

$$n=\frac{{({\sigma }_0{z}_{1-\alpha /2}+{\sigma }_1{z}_{1-\beta })}^2}{{\Delta }^2}\mathrm{.}$$                        (3)

Following Jung, Chow, and Chi (2007), the sample size calculation for the weighted Mantel-Haenszel test can be carried out as follows:

1. Specify the input variables

- Type I and II error probabilities (α β).

- Success probabilities for group 1 p₁₁,...p_J1, and the odds ratio Φ under H₁. Note that $p_{j2}=\Phi p_{j1}/(q_{j1}+\Phi p_{j1})$ .

- Incidence rates for the strata, (a_j, j = 1,...,J). (Yue [2007] proposes to use $a_j\approx 1/J$.

- Allocation probability for group 1 within each stratum, $(b_{j1},j=1,...,J)$.

2. Calculate n by

$$n=\frac{{({\sigma }_0{z}_{1-\alpha /2}+{\sigma }_1{z}_{1-\beta })}^2}{{\Delta }^2},$$

    where

$\Delta =\sum _{j=1}^J{a}_j{b}_{j1}{b}_{j2}(p_{j1}-p_{j2})$

${\sigma }_1^2=\sum _{j=1}^J{a}_j{b}_{j1}{b}_{j2}(b_{j2}{p}_{j1}{q}_{j1}+b_{j1}{p}_{j2}{q}_{j2})$

${\sigma }_0^2=\sum _{j=1}^J{a}_j{b}_{j1}{b}_{j2}(b_{j1}{p}_{j1}+b_{j2}{p}_{j2})(b_{j1}{q}_{j1}+b_{j2}{q}_{j2}).$

Example 18.4

Suppose that we want to compare the response probabilities between control (k = 1) (and experimental (k = 2) (groups. We consider partitioning the combined data into J = 5 strata, and the allocation proportions are projected as and $(a_1,a_2,a_3,a_4,a_5)=(.15,.15,.2,.25,.25)$ and $(b_{11},b_{21},b_{31},b_{41},b_{51})=(.4,.4,.5,.6,.6)$. Also, suppose that the response probabilities for the control group are given as $(p_{11},p_{21},p_{31},p_{41},p_{51})=(.5,.6,.7,.8,.9)$, and we want to calculate the sample size required for a power of 1 − β = 0.8 to detect an odds ratio of $\Phi =2$ using two-sided α = 0.05. For $\Phi =2$, the success probabilities for the experimental group are given as $(p_{12},p_{22},p_{32},p_{42},p_{52})=(.6667,.7500,.8235,.8889,.9474)$. Under these settings, by (3), we need n = 447 for Mantel-Haenszel.

Program 18.4 SAS Code for Weighted Mantel-Haenszel Test with Strata

%macro WMHTestwithStrata ( J = , /* number of strata                         */
                          inA= , /* incidence rates for the strata           */
                          inB= , /* allocation probability for control group */
                          inP1=, /* success probability for control group    */
                          phi =, /* odds ratio under H1                      */
                          power=, /* power                                   */
                          alpha=, /* alpha                                   */
                          sides=  /* 1: one-sided test 2: Two-sided test     */
                         );

 proc iml;
   %let K = 2; /* two groups */
   A=&inA; B =&inB;P1=&inP1;
   P2 =J(&J,1,0);Q1 =J(&J,1,0); Q2 =J(&J,1,0);
   do j=1 to &J;
      Q1[j]=1-P1[j];
      P2[j]=&phi*P1[j]/(Q1[j]+&phi*P1[j]);
      Q2[j]=1-P2[j];
   end;
   z_p1=1-&alpha/&sides; z_alpha = probit(z_p1);
   z_p2=&power; z_beta = probit(z_p2);
   delta = 0; s0_sq= 0; s1_sq= 0;

   do j=1 to &J;
    delta = delta+A[j]*B[j]*(1-B[j])*(P1[j]-P2[j]);
    s1_sq= s1_sq+A[j]*B[j]*(1-B[j])*((1-B[j])*P1[j]*Q1[j]+B[j]*P2[j]*Q2[j]);
    s0_sq= s0_sq+A[j]*B[j]*(1-B[j])*(B[j]*P1[j]+(1-B[j])*P2[j])*(B[j]*Q1[j]+(1-B[j])*Q2[j]);
   end;
   n = (1/(delta**2))*((sqrt(s0_sq)*z_alpha+sqrt(s1_sq)*z_beta)**2);

   print 'Sample Size Calculation';
   print 'Weighted Mantel-Haenszel Test with Strata';
   alpha = α power = &power; phi= φ sides = &Sides;
   print  alpha power phi sides ;
   print A B P1 P2;
   print delta s0_sq s1_sq;
   print n;
  quit;
  run;
 %mend WMHTestwithStrata;

 /*---------------- Run the macro for Example 18.4. ---------------------*/

 %WMHTestwithStrata(
   J = 5 ,                                  /* number of strata        */
   inA= %str({0.15, 0.15, 0.2, 0.25, 0.25}), /*Incidence Rates for Strata*/
   inB= %str({0.4 , 0.4, 0.5, 0.6, 0.6}),  /* Allocation Probability  */
   inP1=%str({0.5 , 0.6, 0.7, 0.8, 0.9}),  /* Success Probability     */
   phi = 2,                                 /* odds Ratio Under H1     */
   power =0.8,                              /* Power                   */
   alpha = 0.05,                            /* Alpha                   */
   sides = 2                                /* Two-sided test          */
  );

 

Output from Program 18.4

 

18.4 Two-Sample Log-Rank Test for Survival Data

Suppose that n_k subjects are accrued to a study from group k(=1,2). For subject i(=1,...,n_k) in group k(=1,2), let T_ki denote the survival variable (that is, time to an event of interest), with marginal cumulative hazard function ${\Lambda }_k(t)$. We want to test the null hypothesis,

$H_0:{\Lambda }_1(t)={\Lambda }_2(t) \hspace{0.17em}for\hspace{0.17em}all\hspace{0.17em} t\ge 0$

against the alternative hypothesis,

$H_a:{\Lambda }_1(t)\ne {\Lambda }_2(t) \hspace{0.17em}for\hspace{0.17em}some\hspace{0.17em} t\ge 0.$

Because a cumulative hazard function uniquely determines the distribution, H₀ implies that $(T_{1i},i=1,...,n_1)$ and $(T_{2i},i=1,...,n_2)$ have the same distributions.

Due to loss to follow up or termination of study before all subjects experience events, survival times are censored for some subjects. Let C_ki denote the censoring time for the subject i in group k. Then, we observe $\{(X_{ki},{\delta }_{ki}),1\le i\le n_k,k=1,2\}$, where $X_{ki}=min(T_{ki},C_{ki})$ and ${\delta }_{ki}=I(T_{ki}\le C_{ki})$. We assume that C_ki and T_1i are independent within each group.

The log-rank test (Peto and Peto, 1972) is given by

$W={\int }_0^{\infty }\frac{Y_1(t)Y_2(t)}{Y_1(t)+Y_2(t)}\{d{\hat{\Lambda }}_1(t)-d{\hat{\Lambda }}_2(t)\},$

where ${\hat{\Lambda }}_k(t)={\int }_0^t{Y}_k{(s)}^{-1}d{N}_k(s)$ is the Nelson-Aalen estimator (Nelson, 1969; Aalen, 1978) of ${\Lambda }_k(t)$, $N_k(t)={\sum }_{i=1}^{n_k}{N}_{ki}(t)$, $N_{ki}(t)={\delta }_{ki}I(X_{ki}\le t)$, $Y_k(t)={\sum }_{i=1}^{n_k}{Y}_{ki}(t)$, and $Y_{ki}(t)=I(X_{ki}\ge t)$. For large n with and $n_k/n=r_k$, $W/\hat{\sigma }$ is approximately normal with mean 0 and variance 1 under H₀, where

${\hat{\sigma }}^2={\int }_0^{\infty }\frac{Y_1(t)Y_2(t)}{\{Y_1(t)+Y_2(t)\}}\{d{N}_1(t)+d{N}_2(t)\},$

Refer to Fleming and Harrington (1991), for example. We reject H₀ if the absolute value of $W/\hat{\Sigma }$ is larger than $z_{1-\alpha /2}$ .

For sample size and power calculation, we assume a proportional hazards model $\Delta ={\lambda }_2(t)/{\lambda }_1(t)$. The power of the log-rank test depends on the number of events, rather than the number of patients. By Rubinstein and colleagues (1981), the number of events D₁ and D₂ in the two groups, required for a two-sided α test to detect a hazard ratio of Δ(<1) with power 1 − β, is given as

${(\frac{log\Delta }{z_{1-\alpha /2}+z_{1-\beta }})}^2=D_1^{-1}+D_2^{-1}\mathrm{.}$                          (4)

In order to calculate D_k, we need to specify the survival distributions under H_a and the common censoring distribution. Suppose that the subjects are accrued at a constant rate for a period of A and followed for an additional period of B after the completion of accrual. The total study period is A + B. Then, assuming no loss to follow up, the censoring distribution is uniform between B and A + B, as in the following:

$G(t)= P (C_{ki}\ge t)=\{\begin{array}{ll}1\hfill & \hspace{0.17em}if\hspace{0.17em} t<B\hfill \\ 1+(B-t)/A\hfill & \hspace{0.17em}if\hspace{0.17em} B\le t\le A+B\hfill \\ 0\hfill & \hspace{0.17em}otherwise\mathrm{.} \hfill \end{array}$

One of the most popular survival models in sample size calculation is the exponential distribution,

$S_k(t)= P (T_{ki}\ge t)=exp(-{\lambda }_{k}t)$

for group k. Then, the probability that a subject in arm k experiences an event is given as

$d_k= P (T_{ki}\le C_i)={\int }_0^{\infty }{S}_k(t)dG(t)=1-\frac{exp(-{\lambda }_{k}B)\{1-exp(-{\lambda }_{k}A)\}}{A{\lambda }_k},$

so that we have $D_k=n_k{d}_k=n{r}_k{d}_k$. By plugging this in (4) and solving with respect to n, we obtain

$n=(\frac{1}{r_1{d}_1}+\frac{1}{r_2{d}_2})(\frac{z_{1-\alpha /2}+z_{1-\beta }}{\log \Delta }{)}^2\mathrm{.}$

We assumed no loss to follow up in calculating d_k, but this assumption can be easily loosened by assuming a distribution for the time to loss to follow up.

In summary, a sample size is calculated as follows:

• In addition to type I error probability α and power 1 − β, specify

(a) ${\lambda }_1,{\lambda }_2=$ hazard rates under H_a (Δ = λ₂/λ₁)

(b) r_k = the prevalence of group k (r₁ + r₂ = 1)

(c) A = accrual period, B = additional follow-up period

• Calculate the probability of an event for a subject in group k(=1,2)

$d_k=1-\frac{\exp (-{\lambda }_{k}B)\{1-\exp (-{\lambda }_{k}A)\}}{A{\lambda }_k},$

• Calculate the sample size

$n=(\frac{1}{r_1{d}_1}+\frac{1}{r_2{d}_2})(\frac{z_{1-\alpha /2}+z_{1-\beta }}{\log \Delta }{)}^2\mathrm{.}$

 

Example 18.5

Suppose that the control group (k = 1) has a median survival time of ${\mu }_1=3$ years, and the experimental group will be considered acceptable if it extends the median survival by 50% (that is, ${\mu }_2=4.5$ years). Under the exponential survival model, the annual hazard rates are ${\lambda }_1=0.231$ and ${\lambda }_2=0.154$ (Δ = 0.667). Also, suppose that (r₁,r₂) = (0.3,0.7), and patients were uniformly accrued for A = 3 years and the final analysis is conducted B = 2 years after the completion of accrual. Then, we have (d₁,d₂) = (0.5455,0.4115), n = 613 and (n₁ = 184,n₂ = 429), and (D₁,D₂) = (100,177). PROC POWER in SAS/STAT(Program 18.5.2) also provides sample sizes for the log-rank test based on Lakatos (1988), which, unlike our formula, calculates the limit of the variance estimator of the log-rank test under the null hypothesis. As demonstrated here, it gives a slightly smaller sample size, n = 593, under the design setting.

Program 18.5.1 SAS Code for Two-Sample Log-Rank Tests for Survival Data

%macro SS_TwoSmpleLogRank(
          Accrual =,   /* accrual period                          */
          Follow  =,   /* additional follow-up period             */
          inR     =,   /* group allocation proportion(a1, a2)     */
          inLambda=,   /* hazard rates under the alternative      */
          alpha   =,   /* alpha                                   */
          power   =,   /* power                                   */
          sides   =    /* 1: One-sided test 2: Two-sided test     */
                     );
 proc iml;
  %let K = 2;  r = &inR; Group = J(&K,1,0); lambda = &inLambda;
  delta = lambda[2]/lambda[1]; d_prob = J(&K, 1,0);
  n = J(&K, 1,0); D = J(&K, 1,0);
  do i=1 to &K;
    group[i] = i;
    d_prob[i]=1-(exp(-lambda[i]*&Follow))*(1-exp(-lambda[i]*&Accrual))
              /(&Accrual*lambda[i]);
  end;
  z_p1=1-&alpha/&sides; z_alpha = probit(z_p1);
  z_p2=&power; z_beta = probit(z_p2);
  Total =
int((1/(r[1]*d_prob[1])+1/(r[2]*d_prob[2]))*((z_alpha+z_beta)/log(delta))**
2)+1;
  do i=1 to &K; n[i]= r[i]*Total; D[i]= d_prob[i]*n[i]; end;
  print ' Sample Size Calculation';
  print 'for Two-Sample Log-Rank Test for Survival Data';
  Accrual=&Accrual; Follow_Up= &Follow;
  print 'Accrual Period : ' &Accrual;
  print ' Follow Up period : ' &Follow;
  print Group lambda delta r d_prob;
  print Group n D;
  print Total;
  quit;
 run;
%mend SS_TwoSmpleLogRank;

/*-------------- Run the macro for Example 18.5.1 ----------------------------*/

%SS_TwoSmpleLogRank(
  Accrual =3,                  /* accrual period                         */
  Follow  =2,                  /* additional follow-up period            */
  inR     =%str({0.3, 0.7}),   /* group allocation proportion(r1, r2)    */
  inLambda=%str({0.231, 0.154}),/* hazard rates under the alternative    */
  alpha   =0.05,               /* alpha                                  */
  power   =0.9,                /* power                                  */
  sides   =2                   /* two-sided test                         */
      );

 

Output from Program 18.5.1

 

Program 18.5.2 SAS Code (PROC POWER) for Example 18.5

      proc power;
          twosamplesurvival test=logrank
          nfractional
          groupmedsurvtimes = 3 | 4.5
          accrualtime = 3
          followuptime = 2
          groupweights = (0.3,0.7)
          alpha = 0.05
          sides=2
          power = 0.9
          ntotal = .
         ;
      run;

 

Output from Program 18.5.2

 

18.5 Two-Sample Longitudinal Data

In a longitudinal study, we measure the outcome of research interest repeatedly from each subject over a time period. When there are two groups of subjects, one popular primary objective is the between-group comparison of change rate in the expected value of a response variable over time. Typically, the repeated measurements within each subject are correlated. Because of the robustness to possible misspecification of the correlation structure among the repeated measurements, the generalized estimating equation (GEE) method has been one of the most popular methods to fit the regression models and to test on the change rates (Liang and Zeger, 1986). In this section, we discuss sample size estimation methods for such testing. We consider cases where the response variable is continuous or dichotomous.

 

18.5.1 Generalized Estimating Equations

Suppose that there are n_k subjects in treatment group k(=1,2), n₁,n₂ = n. Let, for group k, $N_k={\sum }_{i=1}^{n_k}{m}_{ki}$ denote the total number of observations and r_k = n_k/n the allocation proportion (r₁ + r₂ = 1). For subject i(i = 1,...,n_k) in group k, let y_kij denote the outcome variable at measurement time t_kij (j = 1,..., m_ki) with ${\mu }_{kij}= E (y_{kij})$ that is expressed as

$g({\mu }_{kij})=a_k+b_k{t}_{kij},$

where g(.) is a known link function. In other words, we assume that there exists a link function g(.) linearizing the trajectory in response over time. To simplify the discussions, we use the identity link $g(\mu )=\mu$ for continuous outcome variables and the logit link $g(\mu )=\text{log{}\mu \text{/(1-}\mu \text{)}}$ for binary outcome variables, but the generalization to the use of other links is simple.

The coefficient b_k represents the change rate per unit time in mean response if y is continuous and the change rate in log-odds if y is binary. The measurement times may vary subject by subject due to missing measurements, patients’ visits for measurements at unscheduled times, loss to follow up, or other causes. In this chapter, we assume that any missing data is missing completely at random (Rubin, 1976).

The repeated measurements $(y_{kij},1\le j\le m_{ki})$ within each subject tend to be correlated. However, the true correlation structure is usually unknown or of secondary interest. Liang and Zeger (1986) proposed a consistent estimator, called the GEE estimator, based on a working correlation structure. Using either the identity or logit link, the GEE estimator $({\hat{a}}_k,{\hat{b}}_k)$ based on the working independence structure solves U_k(a,b) = 0, where

$U_k(a,b)=\frac{1}{{\sqrt{n}}_k}\sum _{i=1}^{n_k}\sum _{j=1}^{m_{ki}}\{y_{kij}-{\mu }_{kij}(a,b)\}\left(\begin{array}{c}1\\ t_{kij}\end{array}\right)$

and ${\mu }_{kij}(a,b)=g^{-1}(a+b{t}_{kij})$.

18.5.1.1 Continuous Outcome Variable Case

In the continuous outcome variable case, we have a closed form solution to U_k (a,b) = 0

${\hat{b}}_k=\frac{\sum _{i=1}^{n_k}\sum _{j=1}^{m_{ki}}(t_{kij}-{\overline{t}}_k)y_{kij}}{\sum _{i=1}^{n_k}\sum _{j=1}^{m_{ki}}{(t_{kij}-{\overline{t}}_k)}^2}$

and ${\hat{a}}_k={\overline{y}}_k-{\hat{b}}_k{\overline{t}}_k$, where ${\overline{t}}_k=N_k^{-1}{\sum }_{i=1}^{n_k}{\sum }_{j=1}^{m_{ki}}{t}_{kij}$ and ${\overline{y}}_k=N_k^{-1}{\sum }_{i=1}^{n_k}{\sum }_{j=1}^{m_{ki}}{y}_{kij}$.

By Liang and Zeger (1986), as $n\to \infty$, ${\sqrt{n}}_k({\hat{b}}_k-b_k)\to N(0,v_k)$ in distribution. Here v_k is consistently estimated by

${\hat{v}}_k=\frac{\sum _{i=1}^{n_k}{\{\sum _{j=1}^{m_{ki}}(t_{kij}-{\overline{t}}_k){\hat{\varepsilon }}_{kij}\}}^2}{{\left\{\sum _{i=1}^{n_k}\sum _{j=1}^{m_{ki}}{(t_{kij}-{\overline{t}}_k)}^2\right\}}^2}$

where ${\hat{\varepsilon }}_{kij}=y_{kij}-{\hat{a}}_k-{\hat{b}}_k{t}_{kij}$.

18.5.1.2 Binary Outcome Variable Case

Let $p_{kij}={\mu }_{kij}$ in the binary outcome variable case. We have to solve the equation using a numerical method, such as the Newton-Raphson algorithm: at the l-th iteration,

$\left(\begin{array}{c}{\hat{a}}_k^{(l)}\\ {\hat{b}}_k^{(l)}\end{array}\right)=\left(\begin{array}{c}{\hat{a}}_k^{(l-1)}\\ {\hat{b}}_k^{(l-1)}\end{array}\right)+n_k^{-1/2}{A}_k^{-1}({\hat{a}}_k^{(l-1)},{\hat{b}}_k^{(l-1)})U_k({\hat{a}}_k^{(l-1)},{\hat{b}}_k^{(l-1)}),$

where

$A_k(a,b)=-n_k^{-1/2}\frac{\partial U_k(a,b)}{\partial (a,b)}=\frac{1}{n_k}\sum _{i=1}^{n_k}\sum _{j=1}^{m_{ki}}{p}_{kij}(a,b)q_{kij}(a,b)\left(\begin{array}{cc}1& t_{kij}\\ t_{kij}& t_{kij}^2\end{array}\right)$

$p_{kij}(a,b)=g^{-1}(a+b{t}_{kij})=\frac{e^{a+b{t}_{kij}}}{1+e^{a+b{t}_{kij}}}$                 (5)

and $q_{kij}(a,b)=1-p_{kij}(a,b)$.

By Liang and Zeger (1986), for large n, ${\sqrt{n}}_k{({\hat{a}}_k-a_k,{\hat{b}}_k-b_k)}^T$ is asymptotically normal with mean 0 and variance V_k that can be consistently estimated by

$${\hat{V}}_k=A_k^{-1}({\hat{a}}_k,{\hat{b}}_k){\hat{\Sigma }}_k{A}_k^{-1}({\hat{a}}_k,{\hat{b}}_k)$$

where

$${\hat{\Sigma }}_k=\frac{1}{n_k}\sum _{i=1}^{n_k}{\left\{\sum _{j=1}^{m_{ki}}{\hat{\varepsilon }}_{kij}\left(\begin{array}{c}1\\ t_{kij}\end{array}\right)\right\}}^{\otimes 2},$$

${\hat{\varepsilon }}_{kij}=y_{kij}-p_{kij}({\hat{a}}_k,{\hat{b}}_k)$ and $c^{\otimes 2}=c{c}^T$ for a vector . Let ${\hat{v}}_k$ be the (2,2)-component of ${\hat{V}}_k$.

 

18.5.2 Sample Size Calculation

Suppose that we want to test the rate of change between two groups (that is, H₀: b₁ = b₂) . Based on the asymptotic results from the previous section, we can reject H₀: b₁ = b₂ in favor of H₀: b₁ ≠ b₂ when

$\left|\frac{{\hat{b}}_1-{\hat{b}}_2}{\sqrt{{\hat{v}}_1/n_1+{\hat{v}}_2/n_2}}\right|>z_{1-\alpha /2},$                   (6)

Where Z_1−α/2 is the 100(1−α/2) percentile of a standard normal distribution.

In this section, we derive a sample size formula for the two-sided α test (6) to detect $H_1:|b_2-b_1|=\Delta (>0)$ with power 1−β. When designing a study, we usually schedule fixed visit times $t_1<\cdots <t_m$for m repeated measurements from each subject. We often set t₁ = 0 for the baseline measurement time. When the study is conducted, however, the subjects may skip some visit times due to various reasons, which results in missing values, or they may not follow the visit schedule correctly so that the observed visit times may be variable. Jung and Ahn (2003) show through simulations in a continuous outcome variable case that the sample size formula based on fixed measurement times is very accurate even when the observed measurement times are widely distributed around the scheduled times.

By simple algebra, we can derive a sample size formula to detect the specified difference$|b_2-b_1|=\Delta$ with power 1−β,

$n=\frac{{(z_{1-\alpha /2}+z_{1-\beta })}^2(v_1/r_1+v_2/r_2)}{{\Delta }^2},$          (7)

where $v_k={lim}_{n\to \infty }{\hat{v}}_k$. The expression of v_k is slightly different between continuous and binary outcome variable cases, as shown in the following subsections.

In order to allow for missing values, let δ_j denote the proportion of patients with observations at t_j and δ_j the proportion of patients with observations at both t_j and t_j’. Also, let . ${\rho }_{j{j}^{\prime }}= corr (y_{kij},y_{ki{j}^{\prime }})$We assume that the missing pattern and correlation structure are common in two treatment groups. We discuss specific models for missing pattern and correlation structure in Section 18.5.3.

18.5.2.1 Continuous Variable Case

We assume that the continuous outcome variable has a constant variance $var (y_{kij})={\sigma }^2$over time, which is common between two treatment groups. Jung and Ahn (2003) show that is v₁ = v₂ = v expressed as

$$v=\frac{{\sigma }^2(s^2+c)}{s^4},$$

where

$$s^2=\sum _{j=1}^m{\delta }_j{(t_j-\tau )}^2$$$$c=\sum {\sum }_{j\ne j^{\prime }}{\delta }_{j{j}^{\prime }}{\rho }_{j{j}^{\prime }}(t_j-\tau )(t_{j^{\prime }}-\tau )$$

and

$$\tau =\frac{\sum _{j=1}^m{\delta }_j{t}_j}{\sum _{j=1}^m{\delta }_j}\mathrm{.}$$

Hence, (7) is simplified to

$n=\frac{v{(z_{1-\alpha /2}+z_{1-\beta })}^2}{{\Delta }^2{r}_1{r}_2}\mathrm{.}$                    (8)

Note that we do not have to specify the true values for $\{(a_k,b_k),k=1,2\}$but only the difference $|b_2-b_1|=\Delta$in sample size calculation. The sample size is proportional to the variance of measurement error and decreases as r₁ approaches 1/2.

18.5.2.2 Binary Variable Case

Let P_kj = a_k + b_kt_k denote the success probabilities under H₁. Jung and Ahn (2005) show that

$$v_k=\frac{s_k^2+c_k}{s_k^4},$$

where

$$s_k^2=\sum _{j=1}^m{\delta }_j{p}_{kj}{q}_{kj}{(t_j-{\tau }_k)}^2$$$$c_k=\sum {\sum }_{j\ne j^{\prime }}{\delta }_{j{j}^{\prime }}{\rho }_{j{j}^{\prime }}\sqrt{p_{kj}{q}_{kj}{p}_{k{j}^{\prime }}{q}_{k{j}^{\prime }}}(t_j-{\tau }_k)(t_{j^{\prime }}-{\tau }_k),$$

and

$${\tau }_k=\frac{\sum _{j=1}^m{\delta }_j{p}_{kj}{q}_{kj}{t}_j}{\sum _{j=1}^m{\delta }_j{p}_{kj}{q}_{kj}}\mathrm{.}$$

As shown here, the sample size formula under the binary outcome variable case depends on the probabilities (P_kj, j = 1, . . . , m) under H₁, so that we need to specify all regression parameters $\{(a_k,b_k),k=1,2\}$. Let Δ = b₁−b₂ denote the difference in slope between two groups under H₁. If pilot data exist for the control group (group 1), then we may use the estimates as the parameter values, a₁ and b₁. If t = 0 is the baseline, the intercepts in the two groups are set the same (that is, a₁ = a₂) . By setting b₂ = b₁ + Δ, we can specify all regression coefficients under H₁.

If no pilot data are available, we may specify the binary probabilities at the baseline, p₁₁, and at the end of follow up, p_1m , for the control group. Then we obtain (a₁b₁) by

$b_1=\frac{g(p_{1m})-g(p_{11})}{t_m-t_1},$               (9)

$$a_1=g(p_{11})-b_1{t}_1=g(p_{11}).$$

And we set a₂ = a₁(since p₁₁ = p₂₁ at the baseline) and b₂ = b₁ + Δ.

 

18.5.3 Modelling Missing Pattern and Correlation Structure

Calculation of (or) requires projection of the missing probabilities and the true correlation structure.

18.5.3.1 Missing Pattern

In order to specify δ_jj, we need to estimate the missing pattern. If missing at time t_j is independent of missing at time t_j for each patient, then we have δ_jj’ = δ_jj, and we call this type independent missing.

In some studies, subjects missing at a measurement time may be missing at all following measurement times, as in the labor pain study discussed in Example 18.6. This type of missing is called monotone missing. In this case, we have δ_jj’ = δ_j for ($j<j^{\prime }({\delta }_1\ge \cdots \ge {\delta }_m)$). In monotone missing cases, one may want to specify the proportion of patients who will contribute exactly the first j observations, say £j. Then, noting that ${\eta }_j={\Delta }_j-{\Delta }_{j+1}$for j = 1,...,m−1 and μ_m = δ_m , we can obtain δ_j recursively starting from δ_m. Note also that ${\sum }_j{\eta }_j={\Delta }_1$ , which equals 1 if all patients have measurements at the first measurement time t₁.

In summary, we consider two missing patterns as candidates to the true one: (a) independent missing, where ${\delta }_{j{j}^{\prime }}={\delta }_j{\delta }_{j^{\prime }}$ , and (b) monotone missing, where ${\delta }_{j{j}^{\prime }}={\delta }_{j^{\prime }}$ for $j<j^{\prime }$ (${\delta }_1\ge \cdots \ge {\delta }_m$ ).

18.5.3.2 Correlation Structure

Now, we specify the true correlation structure ${\rho }_{j{j}^{\prime }}$ . A reasonable model for ${\varepsilon }_{ij}=y_{ij}-g^{-1}(a+b{t}_{ij})$ may be

$${\varepsilon }_{ij}=u_i+e_{ij},$$

where u_i is a subject-specific error term with variance ${\sigma }_u^2$ and e_ij is a serially correlated within-subject error term with variance ${\sigma }_e^2$ and correlation coefficients ${\tilde{\rho }}_{j{j}^{\prime }}$. Assuming that u_i and e_ij are independent, we have

$$cov ({\varepsilon }_{ij},{\varepsilon }_{i{j}^{\prime }})={\sigma }_a^2+{\sigma }_e^2{\tilde{\rho }}_{j{j}^{\prime }}\mathrm{.}$$

Often, the variation between subjects (${\sigma }_u^2$) is much larger than that within subjects (${\sigma }_e^2$). In this case, we have

$$cov ({\varepsilon }_{ij},{\varepsilon }_{i{j}^{\prime }})\approx {\sigma }_u^2$$

and an exchangeable correlation structure, ${\rho }_{j{j}^{\prime }}=\rho$ for $j\ne j^{\prime }$ , may be a reasonable approximation to the true one.

On the other hand, if the variation within the subject (${\sigma }_e^2$ ) dominates over the variation between subjects (${\sigma }_u^2$ ), then we will have

$$cov ({\varepsilon }_{ij},{\varepsilon }_{i{j}^{\prime }})\approx {\sigma }_e^2{\tilde{\rho }}_{j{j}^{\prime }},$$

so that a serial correlation structure may be a reasonable approximation to the true one. One of the most popular serial correlation structures, especially when measurement times are not equidistant, is a continuous autocorrelation model with order 1, AR(1), for which ${\rho }_{j{j}^{\prime }}={\rho }^{|t_j-t_{j^{\prime }}|}$ .

We consider two correlation structures as candidate approximations to the true one: (i) exchangeable, where ${\rho }_{j{j}^{\prime }}=\rho$ , and (ii) AR(1), where ${\rho }_{j{j}^{\prime }}={\rho }^{|t_j-t_{j^{\prime }}|}$ .

18.5.3.3 Examples

For sample size calculation, the following parameters need to be specified commonly in both continuous and discrete outcome variable cases:

• Type I and II errors α and β, respectively.

• The size of difference to detect, $|b_2-b_1|=\Delta$ .

• Allocation proportion for group k, r_k.

• Correlation structure and the associated correlation parameter ρ, (that is, (i) ${\rho }_{j{j}^{\prime }}=\rho$ for exchangeable or (ii) ${\rho }_{j{j}^{\prime }}={\rho }^{|t_j-t_{j^{\prime }}|}$ for AR(1)).

• Proportion of patients with an observation at t_j, δ_j.

• Missing pattern: (a) independent missing (${\delta }_{j{j}^{\prime }}={\delta }_j{\delta }_{j^{\prime }}$ ) or (b) monotone missing (${\delta }_{j{j}^{\prime }}={\delta }_{j^{\prime }}$ for $j<j^{\prime }$ ) .

In addition, we need to specify ${\sigma }^2= var ({\varepsilon }_{ij})$ in the continuous outcome variable case and the regression coefficients $\{(a_k,b_k),k=1,2\}$ under H₁ in the binary outcome variable case.

We demonstrate our sample size formula with real longitudinal studies.

Example 18.6 Continuous Outcome Case

In a study on labor pain (Davis, 1991), 83 women in labor were assigned to either a pain medication group (43 women) or a placebo group (40 women). At 30-minute intervals, the self-reported amount of pain was marked on a 100mm line, where 0 = no pain and 100 = extreme pain. The maximum number of measurements for each woman was , but there were numerous missing values at later measurement times with monotone missing pattern. A simple approach to such a study objective might be to estimate and compare the slopes of pain scores over time in the two treatment groups. In this study, the outcome variable is continuous. Suppose that we want to design a new study on labor pain based on the data reported by Davis (1991). As in the original study, we assume monotone missing. In this study, the measurement times were equispaced, so that we set t_j = j − 1 (j = 1,…6) for convenience. From the data, we obtained σ² = 815.84. Suppose we want to detect a difference of σ in mean pain score between two groups at t₆. So, we project Δ = σ/(t₆ − t₁ = 5.71 in a new study. We consider a balanced design (that is, r₁ = r₂ = 1/2) . Also, from the data, the proportion of observed measurements are

$$({\delta }_1,{\delta }_2,{\delta }_3,{\delta }_4,{\delta }_5,{\delta }_6)=(1,.90,.78,.67,.54,.41).$$

From these results τ = 2.02, we obtain and s² = 11.42.

Suppose that we want to detect a difference of Δ = 5.71 with 80% of power (z_1−β = 0.84) using a two-sided α = .05 (z_1−α/2 = 1.96) test. Under an exchangeable correlation structure, we obtain ρ = .64 and c = −3.13 from the data, so that we have v = 815.84×(11.42−3.13)/11.42² = 51.84. Hence, from (8), the required sample size is calculated as

$$n=[\frac{51.84\times {(1.96+0.84)}^2}{{5.71}^2}]+1=50.$$

where [x] is the largest integer not exceeding x.

Under AR(1), we obtain from the data ρ = .80 and c = 2.31, so that the required sample size for detecting Δ = 5.71 with 80% of power using a two-sided α = .05 test is given as n = 83. The sample size under AR(1) is larger than that under an exchangeable correlation, by about 66%, in this example.

Program 18.6 SAS Code for Two-Group Comparision of Repeated Continuous Measurement

%macro SS_RepeatedContinuousMeasurement(
               missingPattern = , /* 1: independent missing 2: monotone missing     */
               corrStructure  = , /* 1: compound symetric , 2: AR(1)                */
               rho            = , /* associated correlation parameter               */
               m              = , /* number of measurement time points              */
               sigma_sq       = , /* variance                                       */
               inR            = , /* group allocation proportion(r1, r2)            */
               inDelta        = , /* proportion of observed measurements            */
               alpha          = , /* alpha                                          */
               power          = , /* power                                          */
               sides          = , /* 1: one-sided test  2: two-sided test           */
               print    = 0       /* 0:default 1: detail                            */
                               );
proc iml;
  %let K = 2; r = &inR; delta = &inDelta; t = J(&m,1,0);
  do j=1 to &m; t[j] = j-1; end;
  %if &inDelta eq %then %do; delta=J(&m,1,0);do j=1 to &m; delta[j]=1-(j-1)/20; end; %end;
  d = sqrt(&sigma_sq)/(t[&m]-t[1]);
  start g(p); gp = log(p/(1-p)); return (gp); finish g;
  start prob(a,b,t); p = 1/(1+exp(-a-b*t)); return (p); finish prob;
  start rho(i,j,r,c);
    if c=1 then do; /* CS */ if i=j then rho_ij =1;else rho_ij = r; end;
    else do; /* AR(1) */ dist = abs(i-j); rho_ij = r**dist; end;
    return (rho_ij);
  finish rho;
  tau_num =0;  tau_denum = 0;
  do j=1 to &m; tau_num = tau_num+delta[j]*t[j]; tau_denum = tau_denum+delta[j]; end;
  tau= tau_num/tau_denum; s_sq=0;
  do j=1 to &m; s_sq = s_sq+delta[j]*((t[j]-tau)**2); end;
  c=0;
  do i=1 to &m;
     do j=1 to &m;
           if i ^= j then do;
              if &missingPattern = 1 then do; delta_ij = delta[i]*delta[j]; end;
          else do; if j > i then max_ij=j; else max_ij=i; delta_ij=delta[max_ij]; end;
          c = c+delta_ij*rho(i,j,&rho,&corrStructure)*(t[i]-tau)*(t[j]-tau);
                end;
      end;
   end;
   v = &sigma_sq*(s_sq+c)/(s_sq**2);
   z_p1=1-&alpha/&sides; z_alpha = probit(z_p1); z_p2=&power;  z_beta  = probit(z_p2);
   n = int((v*(z_alpha+z_beta)**2)/(d**2*r[1]*r[2]))+1;
   print 'Sample Size Calculation for a Two-Group Comparision';
   print ' of Repeated Continous Measurements';
   alpha  =&alpha; power=&power; rho=&rho; sides=&sides;
   sigma_sq=&sigma_sq;
   print alpha power rho sides;
   if &missingPattern = 1 then do;print ' Missing Pattern: Independent '; end;
   else if &missingPattern = 2 then do; print ' Missing Pattern: Monotone '; end;
   if &corrStructure=1 then do; print ' Correlation Structure : Compound Symetric ';end;
   else if &corrStructure = 2 then do; print ' Correlation Structure : AR(1) '; end;
   %if &print = 1 %then %do; print delta; print tau s_sq c v ; %end;
   print d sigma_sq;
   print n;
   quit;
run;
%mend SS_RepeatedContinuousMeasurement;

/*--Run the macro for Example 18.6 when Correlation Structure = Compound Symetric --*/

%SS_RepeatedContinuousMeasurement(
       missingPattern = 2,              /* monotone missing                      */
       corrStructure  = 1,              /* compound symetric                     */
       rho            = 0.64,           /* associated correlation parameter      */
       m              = 6,              /* number of measurement time points     */
       sigma_sq       = 815.84,         /* variance                              */
       inR          = %str({0.5, 0.5}), /* group allocation proportion(r1, r2)   */
       inDelta      = %str({1, 0.9, 0.78, 0.67, 0.54, 0.41}),
                                        /* proportion of observed measurements   */
       alpha        = 0.05,             /* alpha                                 */
       power        = 0.8,              /* power                                 */
       sides         = 2,               /* two-sided test                        */
       print         = 1
        );

/*--- Run the macro for Example 18.6 when Correlation Structure = AR(1) ---------*/

%SS_RepeatedContinuousMeasurement(
     missingPattern  = 2,              /*  monotone missing                      */
     corrStructure  = 2,               /* AR(1)                                  */
     rho             = 0.8,            /*  associated correlation parameter      */
     m              = 6,               /* number of measurement time points      */
     sigma_sq        = 815.84,         /* variance                               */
     inR            = %str({0.5, 0.5}),/* group allocation proportion(r1, r2)    */
     inDelta        = %str({1, 0.9, 0.78, 0.67, 0.54, 0.41}),
                                       /* proportion of observed measurements    */
     alpha          = 0.05,            /* alpha                                  */
     power          = 0.8,             /* power                                  */
     sides           = 2,              /* two-sided test                         */
     print = 1
        );

Output from Program 18.6

Example 18.7 Binary Outcome Case

The Genetics vs. Environment In Scleroderma Outcome Study (GENISOS) is an observational study designed as a collaboration of the University of Texas-Houston Health Science Center with the University of Texas Medical Branch at Galveston and the University of Texas-San Antonio Health Science Center (Reveille et al., 2001). Scleroderma, or systemic sclerosis, is a multisystem disease of unknown etiology characterized by cutaneous and visceral fibrosis, small blood vessel damage, and autoimmune features (Medsger, 1997). The study subjects are regularly followed to check for the occurrence of pulmonary fibrosis. In this case, the outcome is a binary variable.

Suppose that we want to develop a study to examine the effect of a new drug in preventing the occurrence of pulmonary fibrosis in subjects with scleroderma compared to no intervention. The parameter values specified here for sample size calculation are approximated by the estimates from the current data set of GENISOS.

We want to estimate the sample size for the new study using α = .05 and 1−β = .8 . As in GENISOS, presence or absence of pulmonary fibrosis is assessed at baseline and at months 6, 12, 18, 24, and 30. Because the measurement times are equidistant, we set t_j = j − 1 j ( j = 1,…,6 ) for the m = 6 time points. The within-group correlation structure of the repeated measurements conforms to AR(1) with the adjacent correlation equal to ρ = .8 (that is, ${\rho }_{j{j}^{\prime }}{=.8}^{|j-j^{\prime }|}$ ). We consider assigning an equal number of scleroderma patients in each of two groups (that is, r₁ = r₂ =1/2).

Approximately 75% of scleroderma patients do not have pulmonary fibrosis at the baseline in the ongoing GENISOS. We project that the proportion of subjects without pulmonary fibrosis is 75% at baseline (that is, p_1,1 = .75) and 50% at 30 months (that is, p_1,j = .50) in a placebo group. We assume that a new therapy will prevent or delay further occurrence of pulmonary fibrosis. That is, the proportion of subjects without pulmonary fibrosis will remain 75% during the 30-month study in a new therapy group (in other words, p_2,1 = p_2,6 = .75 ). From these values and (9), we obtain

$$b_1=\frac{g(.5)-g(.75)}{5-0}=-0.220,$$

and a₁ = g(.75) = 1.099. Similarly, we obtain (a₂, b₂)=(1.099, 0), so that Δ = 0 − (−0.220) = 0.220. By (5), the probabilities of no pulmonary fibrosis at the 6 time points are obtained as (.750,.707,.659,.608,.555,.500) for the placebo group and (.750,.750,.750,.750,.750,.750) for the treatment group.

The proportions of observed measurements are expected to be

$$({\delta }_1,{\delta }_2,{\delta }_3,{\delta }_4,{\delta }_5,{\delta }_6)=(1,.95,.90,.85,.80,.75).$$

Suppose that we expect independent missing. Then, we obtain ${\delta }_{j{j}^{\prime }}$ from the specified δ_j values using ${\delta }_{j{j}^{\prime }}={\delta }_j{\delta }_{j^{\prime }}$ . Now we have all the parameter values required, and we obtain v₁ = 0.305 and v₁ = 0.353. Finally, from (7), we obtain

$$n=[\frac{{(1.96+0.84)}^2(0.305/0.5+0.353/0.5)}{{0.220}^2}]+1=215.$$

If we assume monotone missing ( ${\delta }_{j{j}^{\prime }}={\delta }_{j\vee j^{\prime }}$ ), we obtain v₁ = 0.324 , , v₂ = 0.308 and n = 229 .

Program 18.7 SAS Code for Two-Group Comparision of Repeated Binary Measurement

%macro SS_RepeatedBinaryMeasurement(
       missingPattern =   , /* 1: independent missing 2: monotone missing       */
       corrStructure =    , /* 1: compound symetric , 2: AR(1)                  */
       rho           =    , /*  associated correlation parameter                */
       m             =    , /* number of measurement time points                */
       inP1          =    , /* proportion of subjects for control group         */
       inP2          =    , /* proportion of subjects                           */
       inR           =    , /* group allocation proportion(r1, r2)              */
       inDelta =          , /* proportion of observed measurements              */
       alpha   =          , /* alpha                                            */
       power   =          , /* power                                            */
       sides    =           /* 1: one-sided test  2: two-sided test             */
     );
  proc iml;
   start g(p); gp = log(p/(1-p)); return (gp); finish g;
   start prob(a,b,t); p = 1/(1+exp(-a-b*t)); return (p); finish prob;
   start rho(i,j,r,c);
    if c=1 then do; /* CS */ if i=j then rho_ij =1; else rho_ij = r; end;
    else if c = 2 then do; /* AR(1) */ dist = abs(i-j); rho_ij = r**dist; end;
    return (rho_ij);
   finish rho;
   %let K = 2;
   P1=&inP1;  P2=&inP2; r =&inR; a=J(&K,1,0);  b=J(&K,1,0); t=J(&m,1,0);
   tau = J(&K,1,0);  s_sq = J(&K,1,0); v = J(&K,1,0); c_sq = J(&K,1,0);
   do j=1 to &m;  t[j] = j-1; end;
   %if &inDelta eq %then %do;
       delta = J(&m,1,0);
       do j=1 to &m; delta[j] = 1-(j-1)/20; end;
   %end;
   b[1] = (g(P1[&m])-g(P1[1]))/(t[&m]-t[1]); a[1] = g(P1[1])-B[1]*t[1];
   b[2] = (g(P2[&m])-g(P2[1]))/(t[&m]-t[1]); a[2] = g(P2[1])-B[1]*t[1];
   d = b[2]-b[1];
   do j=1 to &m;
      T[j] = j-1; P1[j]= prob(a[1], b[1], t[j]); P2[j]= prob(a[2], b[2], t[j]);
   end;
   tau1_num =0; tau1_denum = 0; tau2_num =0; tau2_denum = 0;
   do j=1 to &m;
     tau1_num = tau1_num+delta[j]*P1[j]*(1-P1[j])*t[j];
     tau1_denum = tau1_denum+delta[j]*P1[j]*(1-P1[j]);

     tau2_num = tau2_num+delta[j]*P2[j]*(1-P2[j])*t[j];
     tau2_denum = tau2_denum+delta[j]*P2[j]*(1-P2[j]);
   end;
   tau[1]= tau1_num/tau1_denum; tau[2]= tau2_num/tau2_denum;
   s1=0;s2=0;
   do j=1 to &m;
      s1 = s1+delta[j]*P1[j]*(1-P1[j])*((t[j]-tau[1])**2);
      s2 = s2+delta[j]*P2[j]*(1-P2[j])*((t[j]-tau[2])**2);
   end;
   s_sq[1]=s1; s_sq[2]=s2; c1=0; c2=0; c12=0;
   do i=1 to &m;
      do j=1 to &m;
           if i ^= j then do;
              if &missingPattern = 1 then do;
             delta_ij = delta[i]*delta[j];
              end;else do;
                 if j > i then max_ij = j; else max_ij = i; delta_ij = delta[max_ij];
          end;
          c1=c1+delta_ij*rho(i,j,&rho,&corrStructure)*
             sqrt(P1[i]*(1-P1[i])*P1[j]*(1-P1[j]))*((t[i]-tau[1])*(t[j]-tau[1]));
          c2=c2+delta_ij*rho(i,j,&rho,&corrStructure)*
             sqrt(P2[i]*(1-P2[i])*P2[j]*(1-P2[j]))*((t[i]-tau[2])*(t[j]-tau[2]));
         end;
      end;
   end;
   c_sq[1]=c1; c_sq[2]=c2;
   do k =1 to &K; v[k] = (s_sq[k]+c_sq[k])/(s_sq[k]**2); end;
   z_p1=1-&alpha/&sides; z_alpha = probit(z_p1);
   z_p2=&power; z_beta  = probit(z_p2);
   n = int((z_alpha+z_beta)**2*(v[1]/r[1]+v[2]/r[2])/d**2)+1;

   print 'Sample Size Calculation for a Two-Group Comparision';
   print ' of Repeated Binary Measurements';
   alpha  =&alpha; power=&power; rho=&rho; sides=&sides;
   print alpha power rho sides;
   if &missingPattern = 1 then do; print ' Missing Pattern: Independent '; end;
   else if &missingPattern = 2 then do;print ' Missing Pattern: Monotone '; end;
   if &corrStructure = 1 then do; print ' Correlation Structure : Compound Symetric ';
   end; else if &corrStructure = 2 then do; print ' Correlation Structure : AR(1) ';
   end;
   print P1 P2 delta v ; print d; print n;
   quit;
run;
%mend  SS_RepeatedBinaryMeasurement;

/*--------------- Run the macro for Example 18.7 ------------------------------*/

%SS_RepeatedBinaryMeasurement(
     missingPattern = 1,          /* 1: independent missing 2: monotone missing     */
     corrStructure  = 2,          /* 1: compound symetric , 2: AR(1)                */
     rho     = 0.8,               /*  associated correlation parameter              */
     m       = 6,                 /* number of measurement time points              */
    inP1    = %str({0.75, 0, 0, 0, 0, 0.5}),
                                  /* proportion of subjects for control group       */
    inP2    = %str({0.75, 0.75, 0.75, 0.75, 0.75, 0.75}), 
                                  /* proportion of subjects                         */
     inR     = %str({0.5, 0.5}),  /* group allocation proportion(r1, r2)            */
     inDelta = ,      /* delta[j] = 1-(j-1)/20, proportion of observed measurements */
     alpha   = 0.05,              /* alpha                                          */
     power   = 0.8,               /* power                                          */
     sides   = 2                  /* 1: one-sided test  2: two-sided test           */
     );

%SS_RepeatedBinaryMeasurement(
    missingPattern = 2,          /*  monotone missing                               */
    corrStructure = 2,           /* AR(1)                                           */
     rho     = 0.8,                 /* associated correlation parameter             */
     m       = 6,                   /* number of measurement time points            */
    inP1   = %str({0.75, 0, 0, 0, 0, 0.5}),   

                                    /* proportion of subjects for control group     */
    inP2    = %str({0.75, 0.75, 0.75, 0.75, 0.75, 0.75}), 
                                 /* proportion of subjects                          */
     inR     = %str({0.5, 0.5}),    /* group allocation proportion(r1, r2)          */
     inDelta = ,      /* delta[j] = 1-(j-1)/20, proportion of observed measurements */
     alpha   = 0.05,                /* alpha                                        */
     power   = 0.8,                 /* power                                        */
     sides   = 2                    /* 1: one-sided test  2: two-sided test         */
     );

 

Output from Program 18.7

 

18.6 Discussion

In designing an observational study, we should choose a statistical testing method that will be used when analyzing the final data. For an accurate sample size calculation, the sample size formula should be directly derived from the power function of the chosen testing method.

There are two types of parameters included in sample size formulas: primary and nuisance parameters. The primary parameters are those shown in the statistical hypotheses (for example, means, binomial proportions, and hazard rates), and the nuisance parameters (for example, prevalence rates, accrual period in survival data, correlation coefficients, and missing type and probabilities in longitudinal data) are those of no or secondary interest. Because both types of parameters determine the final sample size, it is important to specify the parameter values as accurately as possible. If a database exists, we usually estimate the values of nuisance parameters and the primary parameter value of the control group from the database. Although r_k may be chosen by the prevalence of a disease, we may decide to use different allocation proportions from the prevalence rates in the natural population depending on the relative cost to accrue subjects in two groups. The primary parameter value for the case group is usually chosen based on the clinical significance of the associated intervention on the outcome variable, which is measured in terms of the difference or ratio of the primary parameter values between the two groups. When there is uncertainty in some of these parameters, we may conduct a sensitivity analysis on the power to demonstrate that the power would be adequate under a range of scenarios.

 

Appendix: Asymptotic Distribution of Wilcoxon Rank Sum Test under H_α

Let ${\tilde{X}}_{ki}=X_{ki}-{\theta }_k$. Then ${\tilde{X}}_{11},\dots ,{\tilde{X}}_{1n_1},{\tilde{X}}_{21},\dots ,{\tilde{X}}_{2n_2}$ are IID with cumulative density function (CDF) F(x). Noting that

$$W=\frac{1}{n_1{n}_2}\sum _{i=1}^{n_1}\sum _{j=1}^{n_2}I({\tilde{X}}_{1i}>{\tilde{X}}_{2j}-\Delta ),$$

we have $W={\int }_{-\infty }^{\infty }{\widehat{F}}_2(x+\Delta )d{\widehat{F}}_1(x)$, where ${\widehat{F}}_k(x)=n_k^{-1}{\sum }_{i=1}^{n_k}I({\tilde{X}}_{ki}\le x)$ is the empirical CDF of ${\tilde{X}}_{k1},\dots ,{\tilde{X}}_{k{n}_k}$ that uniformly converges to F(x) as $n_k\to \infty$ . Let

${\nu }_a= P ({\tilde{X}}_{1i}>{\tilde{X}}_{2j}-\Delta )={\int }_{-\infty }^{\infty }F(x+\Delta )dF(x)$ . For a large n, $W-{\nu }_a$ is expressed as

$${\int }_{-\infty }^{\infty }{\widehat{F}}_2(x+\Delta )d{\widehat{F}}_1(x)-{\int }_{-\infty }^{\infty }F(x+\Delta )dF(x)$$

$$={\int }_{-\infty }^{\infty }F(x+\Delta )d\{{\widehat{F}}_1(x)-F(x)\}+{\int }_{-\infty }^{\infty }\{{\widehat{F}}_2(x+\Delta )-F(x+\Delta )\}dF(x)+o_p(n^{-1}),$$

where $o_p(n^{-1})={\int }_{-\infty }^{\infty }\{{\widehat{F}}_2(x+\Delta )-F(x+\delta )\}d\{{\widehat{F}}_1(x)-F(x)\}$ is negligible for a large n.

Because

$${\int }_{-\infty }^{\infty }F(x+\Delta ){\widehat{F}}_1(x)=\frac{1}{n_1}\sum _{i=1}^{n_1}F({\tilde{X}}_{1i}+\Delta )$$

and

$${\int }_{-\infty }^{\infty }F(x+\Delta ){\widehat{F}}_1(x)=\frac{1}{n_1}\sum _{i=1}^{n_1}F({\tilde{X}}_{1i}+\Delta )$$

we have $var (W)={\sigma }_1^2/n_1+{\sigma }_2^2/n_2$, where

$${\sigma }_1^2= var \{F({\tilde{X}}_{1i}+\Delta )\}$$

$${\sigma }_2^2= var \{\overline{F}({\tilde{X}}_{2i}-\Delta )\}= var \{F({\tilde{X}}_{2i}-\Delta )\}.$$

 

References

Aalen, O. 1978. “Nonparametric inference for a family of counting processes.” Annals of Statistics 6: 701–726.

Carroll, R. J., and D. Ruppert. 1988. Transformation and Weighting in Regression. New York: Chapman and Hall/CRC.

Davis, C. S. 1991. “Semi-parametric and non-parametric methods for the analysis of repeated measurements with applications to clinical trials.” Statistics in Medicine 10: 1959–1980.

Fleming, T. R., and D. P. Harrington. 1991. Counting Processes and Survival Analysis. New York: John Wiley & Sons, Inc.

Gart, J. J. 1985. “Approximate tests and interval estimation of the common relative risk in the combination of 2 × 2 tables.” Biometrika >72: 673–677.

Jung, S. H., and C. Ahn. 2003. “Sample size estimation for GEE method for comparing slopes in repeated measurements data.” Statistics in Medicine >22: 1305–1315.

Jung, S. H., and C. W. Ahn. 2005. “Sample size for a two-group comparison of repeated binary measurements using GEE.” Statistics in Medicine 24(17): 2583–2596.

Jung, S. H., S. C. Chow, and E. Chi. 2007. “A note on sample size calculation based on propensity analysis in non-randomized trials.” Journal of Biopharmaceutical Statistics >17: 35–41.

Lakatos, E. 1988. “Sample sizes based on the log-rank statistic in complex clinical trials.” Biometrics 44: 229–241.

Liang, K. Y., and S. L. Zeger. 1986. “Longitudinal data analysis using generalized linear models.” Biometrika >73: 13–22.

Mann, H. B., and D. R. Whitney. 1947. “On a test of whether one of two random variables is stochastically larger than the other.” The Annals of Mathematical Statistics 18: 50-60.

Medsger, Jr., T. A. 1997. “Systemic sclerosis (scleroderma): clinical aspects. In  Arthritis and allied conditions: a textbook of rheumatology. 13th ed. Edited by W. J. Koopman. Baltimore, MD: Williams &  Wilkins, pp 1433–1464.

Nelson, W. 1969. “Hazard plotting for incomplete failure data.” Journal of Quality Technology 1(1): 27–52.

Peto, R., and J. Peto. 1972. “Asymptotically efficient rank invariant test procedures (with discussion).” Journal of the Royal Statistical Society, Series A: General 135: 185–206.

Reveille, J. D., M. Fischbach, T. McNearney, et al. 2001. “Systemic sclerosis in 3 US ethnic groups: a comparison of clinical, sociodemographic, serologic, and  immunogenetic determinants.” Seminars in Arthritis and Rheumatism >30(5): 332–346.

Rosenbaum, P. R., and D. B. Rubin. 1983. “The central role of the propensity score in observational studies for causal effects.” Biometrika >70: 41–55.

Rosenbaum, P. R., and D. B. Rubin. 1984. “Reducing bias in observational studies using subclassification on the propensity score.” Journal of the American Statistical Association 79: 516–524.

Rubin, D. B. 1976. “Inference and missing data.” Biometrika >63: 581–590.

Rubinstein, L. V., M. H. Gail, and T. J. Santner. 1981. “Planning the duration of a comparative clinical trial with loss to follow-up and a period of continued observation.” Journal of Chronic Diseases >34(9/10): 469–479.

SAS Institute Inc. 2009. SAS/STAT9.2 User’s Guide, Second Edition. Cary, NC: SAS Institute Inc.

Yue, L. Q. 2007. “Statistical and regulatory issues with the application of propensity score analysis to non-randomized medical device clinical studies.” Journal of Biopharmaceutical Statistics 17: 1–13.