Introduction
Determining the sample size that is necessary to achieve a certain target power
is a fundamental step in the planning phase of every clinical trial.
The sample size depends on the type 1 error rate \(\alpha\), the type 2 error rate
\(\beta\), the effect size \(\Delta\) and so-called nuisance parameters,
which are parameters that affect the distribution of the test
statistic but are not of interest in the test problem.
While the type 1 and type 2 error rates are usually predetermined and
the minimal clinically important effect size is known, there is often uncertainty
about the magnitude of the nuisance parameters, such as the variance of the data
\(\sigma^2\) for tests with continuous outcomes or the overall response rate \(p\)
for tests with binary outcomes.
Consider, as an example, the meta analysis by Nakata et al. (2018) that compares
minimally invasive preservation with splenectomy during distal pancreatectomy.
Among others, the overall morbidity of the two groups is compared and data
from 13 studies are reported (cf. Figure 2(c) in Nakata et al. (2018)).
Within these 13 studies, overall morbidity rates pooled over both groups
between 0.10 and 0.62 are reported.
This illustrates the high uncertainty about the “true” overall morbidity rate,
which is the nuisance parameter in this setting.
In these cases, an internal pilot study design with blinded sample size recalculation
can be used. In such a design, the nuisance parameter is estimated in a blinded way
(i.e., without using information about the group assignment of the patients)
after a certain number of outcome data is available, and the sample size is
recalculated using this information (Wittes and Brittain 1990).
While in principle blinded sample size recalculation could be done
without any a priori sample size calculation, it is still advisable to calculate
an initial sample size based on the best guess for the nuisance parameter available
in the planning phase and to determine when to recalculate the sample size
based on this initial calculation.
This is done to avoid conducting the recalculation too early
(so that there is still a great uncertainty about the magnitude of the nuisance
parameter when recalculation is performed), or too late
(so that there may be no room for adjusting the sample size any longer as the
recalculated sample size is already exceeded).
Using this method to recalculate the sample size is an attractive option because
the cost in terms of additional sample size is very small (depending on the outcome)
and in most scenarios the type 1 error
rate is unaffected by the blinded sample size recalculation.
Hence, whenever there is uncertainty about the value of a nuisance parameter and
the logistics of the trial allow it, blinded sample size recalculation can be used.
Meanwhile, this is even recommended by regulatory authorities.
For instance, the Committee for Medical Products for Human Use (CHMP) (2006) states that “(w)henever possible, methods for blinded
sample size reassessment (…) should be used”.
Methods to reassess the sample size in a blinded manner in an internal
pilot study design have been developed for a variety of outcomes.
Based on the early work by Stein (1945), Wittes and Brittain (1990) introduced the internal
pilot study design for continuous outcomes.
Their work was extended in different manners by different authors
(cf. among others Birkett and Day (1994), Denne and Jennison (1999), and Kieser and Friede (2000)).
In all these papers, the main task is to re-estimate the variance of a
continuous outcome in a blinded way.
These ideas can be applied to binary outcomes as well where the re-estimated
nuisance parameter is the overall response rate over both treatment arms.
Associated methods were, for instance, presented by Gould (1992) and Friede and Kieser (2004) for
superiority trials and by Friede, Mitchell, and Müller-Velten (2007) for non-inferiority trials.
However, despite the clear benefit of a blinded sample size recalculation and a
great number of publications on that topic, blindrecalc is to the
knowledge of the authors the first R-package on CRAN,
and thus a freely available software,
that helps with the planning of a clinical trial with such a design
by computing the operating characteristics and the distribution of the total
sample size of the study.
The package can be used for pre-planned and midcourse implemented
blinded sample size reassessments in order to evaluate the potential scenarios
the blinded sample size re-estimation may imply.
For continuous outcomes, we implemented the \(t\)-test for superiority trials
and the shifted \(t\)-test for non-inferiority trials.
For binary outcomes, we implemented the chi-squared test for superiority trials
and the Farrington-Manning test for non-inferiority trials.
The structure of the paper is as follows:
In the Statistical methods section, we explain the general way of proceeding
when conducting a trial with an internal pilot study and how to obtain a blinded
estimate of the nuisance parameter for continuous and binary outcomes.
The structure of the package is introduced in Package structure.
We demonstrate how blindrecalc can be utilized to plan a trial with an
internal pilot study design and blinded sample size recalculation in Usage and example.
In Development principles, we outline the principles of the development process
and how we ensure the quality of our code.
Finally, a brief Conclusion complements this paper.
Statistical methods
The general procedure for planning and conducting a trial with a blinded sample size
recalculation is as follows: At first, an initial sample size \(n_{init}\) is calculated
by using a best guess for the value of the nuisance parameter.
The sample size for the first stage of the trial, \(n_1\),
is then calculated as a fraction of \(n_{init}\), e.g.,
0.25, 0.5 or 0.75. After \(n_1\) observations are available, the total sample size
\(n_{rec}\) is recalculated in a blinded way based on the available data. The final
total sample size \(n\) is then determined as:
\[ n = \min(\max(n_1, n_{rec}), n_{max}),\]
where \(n_{max}\) is a prespecified maximal sample size. This is called the unrestricted
design with upper boundary.
A restricted design would use \(n_{init}\) as a lower boundary (Wittes and Brittain 1990).
Often \(n_{max}\) is set to a multiple of \(n_{init}\). The
special case of \(n_{max} = \infty\) results in the unrestricted design (Birkett and Day 1994).
After the final total sample size is calculated, the \(n_2 = n - n_1\) observations for
the second stage are gathered. Finally, the specified statistical test can be conducted
with the data of all \(n\) patients.
In the following, we shortly introduce the implemented tests and how to obtain
a blinded estimate of the nuisance parameter in each case.
Continuous outcomes
Assume a clinical two-arm trial with normally distributed outcomes where a higher
value is deemed to be favorable, with mean values \(\mu_E\) (experimental group) and
\(\mu_C\) (control group) and common unknown variance \(\sigma^2\).
The outcome of interest is the mean difference \(\Delta := \mu_E - \mu_C\).
By introducing a non-inferiority margin \(\delta>0\), the test problem is given by
\[H_0: \Delta \leq - \delta \text{ vs. } H_1: \Delta > -\delta.\]
The null hypotheses can be tested by a shifted \(t\)-test taking the non-inferiority
margin \(\delta\) into account.
Note that the special case \(\delta=0\) corresponds to the standard \(t\)-test
for superiority.
The approximate total sample size for a one-sided \(t\)-test to detect a mean difference
of \(\Delta = \Delta^* > -\delta\) with a power of \(1-\beta\) while controlling
the type 1 error rate at level \(\alpha\) equals
\[n=n_E+n_C=\frac{(1+r)^2}{r} \frac{(z_{1-\alpha/2} + z_{1-\beta})^2 \sigma^2}{(\Delta^*+\delta)^2}.\]
Here, \(r\) refers to the allocation ratio of the sample sizes between the experimental
and the control group, i.e., \(r = n_E / n_C\), and \(z_{1-q}\) denotes the
\(1-q\) quantile of the standard normal distribution.
In this framework, the nuisance parameter is the unknown variance \(\sigma^2\).
Due to potential uncertainty on the value of \(\sigma^2\), it seems
appropriate to re-estimate it in a blinded interim analysis to ensure that
the desired power level is met.
Inserting \(\widehat{\sigma}^2\) observed mid-course
into the above sample size formula may lead to a more reasonable sample size
for the respective trial than sticking to the value assumed in the planning stage.
There exist different methods for estimating \(\sigma^2\) in a blinded manner
(Zucker et al. 1999).
In blindrecalc, the one-sample variance estimator is implemented.
It is defined as
\[
\widehat{\sigma}^2 := \frac{1}{n_1-1} \sum_{j \in \{E, C \}} \sum_{k=1}^{n_{1,j}}
(x_{j,k} - \bar{x} )^2,
\]
where \(x_{j,k}\) is the outcome of patient \(k\) in group \(j\), \(n_{1,j}\) denotes
the first-stage sample size in group \(j\), i.e., \(n_1 = n_{1,E} + n_{1,C}\), and
\(\bar{x}\) equals the mean over all \(n_1\) observations.
Since the patient’s group allocation is not considered when computing
\(\widehat{\sigma}^2\), blinding is maintained when using this variance estimator.
In the superiority case, i.e., if \(\delta=0\), blinded sample size reassessment
can be performed without relevant type 1 error rate inflation (Kieser and Friede 2003).
In the non-inferiority case, however, the type 1 error rate may be inflated
by the internal pilot study design and a correction of the applied significance level
may become necessary to protect the nominal type 1 error rate at level
\(\alpha\) (Friede and Kieser 2003).
In particular, this inflation arises if the sample size recalculation
is performed too early, i.e., if \(n_1\) is chosen too small.
Interestingly, the cost of this procedure in terms of sample size is quite low.
Since the one-sample variance estimate slightly overestimates the variance,
an increase in sample size arises.
However, this increase amounts to only 8 patients with \(\alpha = 0.025\) and
\(\beta = 0.2\) or 12 patients with the same significance level and \(\beta = 0.1\)
(Friede and Kieser 2001).
In return, the sample size recalculation procedure implies that the trial’s
power meets the target value \(1-\beta\) for a wide range of values of \(\sigma^2\).
Lu (2016) gives closed formulas for the exact distribution of the test statistic
of the two-sample \(t\)-test in this setting.
This allows the simulation of error probabilities and of the sample size distribution
in an acceptable amount of time.
The proposals made by Lu (2016) are implemented in blindrecalc.
Thus, the design characteristics for continuous outcomes presented in
blindrecalc are obtained by simulation and not by exact computation.
This is the case for for binary outcomes that are presented in the following.
Binary outcomes
In a superiority trial with binary outcomes where a higher response probability
is assumed to be favorable, the one-sided null and alternative hypothesis are
\[H_0: p_E \leq p_C \text{ vs. } H_1: p_E > p_C,\]
where \(p_E\) and \(p_C\) denote the event probabilities in the experimental
and the control group, respectively.
While several tests exist for this test problem, the widely used chi-squared test
is implemented in blindrecalc.
The sample size for this test can be approximated with the formula (Kieser 2020):
\[n = \frac{1 + r}{r} \frac{\left(z_{1 - \alpha / 2}
\sqrt{(1 + r) \cdot p_0 \cdot (1 - p_0)} + z_{1 - \beta}
\sqrt{r \cdot p_{C,A} \cdot (1 - p_{C,A}) +
p_{E,A} \cdot (1 - p_{E,A})}\right)^2}
{\Delta^2}.\]
Again, \(r\) denotes to the allocation ratio of the sample sizes, and
\(z_{1 - \alpha / 2}\) and \(z_{1 - \beta}\) are the \(1 - \alpha / 2\) and the
\(1 - \beta\) quantiles of the standard normal distribution.
Furthermore, \(p_{C,A}\) and \(p_{E,A}\) are the
response probabilities in the control and the experimental group under
the assumed alternative, \(p_0\) is the overall response probability,
i.e., \(p_0 = (p_{C,A} + r \cdot p_{E,A}) / (1 + r)\), and \(\Delta\) is the effect
under the alternative, i.e. \(\Delta = p_{E,A} - p_{C,A}\).
The nuisance parameter here is \(p_0\), which can be estimated
in a blinded way after \(n_1\) observations with
\[\hat{p}_0 = \frac{X_{1,E} + X_{1,C}}{n_{1,E} + n_{1,C}},\]
where \(X_{1,E}\) and \(X_{2,C}\) denote the number of observed events in the
experimental and the control group and \(n_{1,E}\) and \(n_{1,C}\) represent
the first-stage sample sizes in the two groups. Blinded estimates of the
event rates in each group can then be obtained by
\(\hat{p}_{C,A} = \hat{p}_0 - \Delta \cdot r / (1 + r)\) and
\(\hat{p}_{E,A} = \hat{p}_0 + \Delta / (1 + r)\).
These estimates are used to recalculate the sample size.
The benefit of the blinded recalculation is that the desired power can be maintained,
even if the initially assumed value for \(p_0\) was wrong.
It is well known that the chi-squared test in a fixed design does not
maintain the nominal significance level, hence the same can be expected for a
chi-squared test with a blinded sample size recalculation.
In fact, Friede and Kieser (2004) showed that the actual levels of the test with and without
recalculating the sample size are very close.
In a non-inferiority trial, the null and alternative hypothesis are
\[H_0: p_E - p_C \leq - \delta \text{, } H_1: p_E - p_C > -\delta,\]
where \(\delta > 0\) is the fixed non-inferiority margin.
The most commonly used test for this problem was proposed by Farrington and Manning (1990).
An approximate sample size formula for this test is
\[n = \frac{1 + r}{r} \cdot \frac{\left(z_{1 - \alpha / 2}
\sqrt{r \cdot \tilde{p}_C \cdot (1 - \tilde{p}_C) +
\tilde{p}_E \cdot (1 - \tilde{p}_E)} + z_{1 - \beta}
\sqrt{r \cdot p_{C,A} \cdot (1 - p_{C,A}) +
p_{E,A} \cdot (1 - p_{E,A})} \right)^2}{(\Delta + \delta)^2},\]
where \(\tilde{p}_{E}\) and \(\tilde{p}_{C}\) are large sample approximations
of the restricted maximum likelihood estimators under the null hypothesis
restriction \(p_E - p_c = - \delta\) (see the Appendix of Farrington and Manning (1990) for the computation).
The same formulas as for the chi-squared test can be used to estimate
\(\hat{p}_E\) and \(\hat{p}_C\) in a blinded way, and these estimates have to be
used to obtain blinded estimates \(\hat{\tilde{p}}_E\) and \(\hat{\tilde{p}}_C\)
for the restricted maximum likelihood estimates.
Plugging these estimates into the sample size formula gives the re-estimated sample size.
Like the chi-squared test, the Farrington-Manning test is also no exact test
and can exceed the nominal significance level.
Friede, Mitchell, and Müller-Velten (2007) showed that in general no further inflation of the type 1 error rate
is caused by blinded re-estimation of the sample size.
Nevertheless, it is possible for the chi-squared test as well as for the
Farrington-Manning test to choose the nominal significance level smaller than
\(\alpha\) in order to protect the type 1 error rate at level \(\alpha\).
Such an adjustment of \(\alpha\) is implemented in blindrecalc for the
binary and the continuous case.
Package structure
When a clinical trial with an internal pilot study is planned, it is essential
to know the characteristics of the applied design.
To this end, the performance in terms of achieved power levels, type I error rates,
and sample size distribution has to be known for different values of the
nuisance parameter and the first-stage sample size \(n_1\).
The package blindrecalc provides all necessary tools that are needed to plan
a trial with a blinded sample size recalculation with only a small number of
functions, which makes using the package very accessible.
blindrecalc makes use of R’s S4 class system.
This allows the application of the same methods for different design classes
and facilitates the usage of the package.
Furthermore, this approach makes the package easily extendable without
any changes in the current source code.
The usage of blindrecalc is intended to be as intuitive as possible.
To obtain characteristics of a blinded sample size recalculation procedure,
two steps have to be made.
At first, the user has to define a design object to indicate which test
and which characteristics such as the desired type 1 and type 2 error rates
are to be applied.
To this end, the three functions setupChiSquare,
setupFarringtonManning, and setupStudent exist to define
a design object of the class corresponding to the respective test.
Secondly, the trial characteristic of interest can be calculated.
Currently, the following methods are implemented:
The method toer allows the computation of the actual type 1 error rate
for different values of the nuisance parameter and the sample size of
the internal pilot study.
By means of adjusted_alpha, the adjusted significance level
can be calculated that can be applied as nominal significance level
when strict type 1 error rate control is desired.
The method pow computes the achieved power of the design under a given
set of nuisance parameters or internal pilot sample sizes.
With n_fix, the sample size of the corresponding fixed design
can be computed.
Finally, the method n_dist provides plots and summaries of the
distribution of the sample size.
For all these methods (except for n_dist), the logical parameter
recalculation allows to define whether a fixed design or a design
with blinded sample size recalculation is analyzed.
Usage and example
For each test, there is a setup function (e.g., setupChiSquare
for the chi-squared test) that creates an object of the class of the test.
Each setup function takes the same arguments:
alpha: The one-sided type 1 error rate.beta: The type 2 error rate.r: The allocation ratio between experimental and control group,
with a default of 1.delta: The difference in effect size between alternative and null hypothesis.alternative: Whether the alternative hypothesis contains greater
(default) or smaller values than the null.n_max: The maximal total sample size, with a default value of Inf.
In this example, the nuisance parameter is the overall response rate \(p\).
A difference in response rates between the two treatment groups of
\(\Delta = p_{E} - p_{C} = 0.2\) is to be detected.
Using blindrecalc, a chi-squared test that achieves a power of
\(1 - \beta = 0.8\) to detect this effect of \(\Delta = 0.2\) and that uses a nominal
type 1 error rate of \(\alpha = 0.025\) can be set up by
design <- setupChiSquare(alpha = 0.025, beta = 0.2, delta = 0.2)
The sample size for a fixed design given one or multiple values of the
nuisance parameter (argument nuisance) can then be calculated with
the function n_fix:
n_fix(design, nuisance = c(0.2, 0.3, 0.4, 0.5))
#> [1] 124 164 186 194
The function toer calculates the actual level of a design with blinded
sample size recalculation or of a fixed design (logical argument recalculation)
given either one or more values of the total sample size in a fixed or the
sample size for the first stage in a recalculation design (argument n1)
or one or more values of the nuisance parameter.
Note that all functions are only vectorized in one of the two arguments
n1 and nuisance.
In this example, it is assumed that the internal pilot study contains half of
the fixed sample size that would be needed if the overall response rate \(p\)
equals \(0.2\).
In this setting, blindrecalc can be used to compare the actual levels
of a fixed design and a recalculation design with the same parameters.
n <- n_fix(design, nuisance = 0.2)
p <- seq(0.1, 0.9, by = 0.01)
toer_fix <- toer(design, n1 = n, nuisance = p, recalculation = FALSE)
toer_ips <- toer(design, n1 = n/2, nuisance = p, recalculation = TRUE)
In the figure below, the type 1 error rate in dependence of the
nuisance parameter is depicted for the designs with and without sample size
recalculation.
Note that, as mentioned in Section Binary outcomes, the level of significance
exceeds the pre-defined level of \(\alpha = 0.025\) in both cases.
If strict control of the type 1 error rate is desired, the function
adjusted_alpha can be used to to calculate an adjusted significance level,
such that the nominal significance level is preserved.
adj_sig <- adjusted_alpha(design, n1 = n/2, nuisance = p, precision = 0.0001,
recalculation = TRUE)
design@alpha <- adj_sig
toer_adj <- toer(design, n1 = n/2, nuisance = p, recalculation = TRUE)
In this example, the adjusted significance level equals 0.0232
for the trial with internal pilot study, i.e., using this value as nominal level
ensures that the actual significance level does not exceed \(\alpha = 0.025\).
The figure demonstrates that the type 1 error rate
is protected at level \(\alpha = 0.025\) if the adjusted significance level
is applied.
In the setting of binary outcomes, adjusting the level such that the
nominal type 1 error rate is protected for any realization of the nuisance parameter
in its domain \([0, 1]\) is feasible.
However, when the nuisance parameter has an infinite domain,
such as the variance in the case of continuous outcomes, this is not possible.
The solution in these cases is to compute a \((1 - \gamma)\)
confidence interval of the nuisance parameter in the blinded interim analysis and
adjust the significance level such that the actual level is below \(\alpha - \gamma\)
for all values in this confidence interval (Friede and Kieser 2011).
If this approach is applied, the user can set the parameter gamma.
Actual significance level for different nuisance parameters with (right panel) and without (left panel) adjustment of the nominal significance level. In the adjusted case, the actual level is below the desired level (red line) for all nuisance parameters in contrast to the un-adjusted case.
To calculate the power of either the internal pilot study design or the fixed design,
the function pow can be used.
Again, the function is vectorized in either n1 or nuisance.
This function can be used to compare the power values of the two designs under
different actual values of the nuisance parameter.
pow_fix <- pow(design, n1 = n, nuisance = p, recalculation = FALSE)
pow_ips <- pow(design, n1 = n/2, nuisance = p, recalculation = TRUE)
Power values for different nuisance parameters for a fixed design and a design with blinded sample size recalculation. The design with recalculation meets the target power (red line) for a wider range of nuisance parameters.
As we can see in the figure, the power achieved by the internal pilot
study design is very close to the target power of 0.8 in most cases.
Only when the overall response rate is very close to 0 or 1, the power is exceeded.
On the other hand, the fixed design is much more sensitive to the actual value
of the nuisance parameter and the actual power can either be way too large
or way too small if the sample size was calculated under wrong assumptions.
Finally, the distribution of the total sample size can be computed under
different assumptions on the nuisance parameter with the function n_dist.
This is particularly useful for the planning of internal pilot study designs
since it allows the investigation of what could happen in a certain clinical
trial and helps the applicant to prepare for different scenarios.
p <- seq(0.2, 0.8, by = 0.1)
n_dist(design, n1 = n/2, nuisance = p, plot = TRUE)
Distribution of the sample size for a design with blinded sample size recalculation in dependence of the nuisance parameter. For more extreme values of the nuisance parameter the variance of the sample size distribution becomes larger.
#> p = 0.2 p = 0.3 p = 0.4 p = 0.5 p = 0.6 p = 0.7 p = 0.8
#> Min. 62.000 78.0000 124.0000 158.0000 124.0000 78.0000 62.000
#> 1st Qu. 108.000 152.0000 182.0000 196.0000 182.0000 152.0000 108.000
#> Median 124.000 170.0000 192.0000 198.0000 192.0000 170.0000 124.000
#> Mean 125.021 164.6057 188.4801 196.5075 188.4801 164.6057 125.021
#> 3rd Qu. 138.000 178.0000 196.0000 200.0000 196.0000 178.0000 138.000
#> Max. 190.000 200.0000 200.0000 200.0000 200.0000 200.0000 190.000
By default, n_dist prints a summary of the sample size distribution
for each nuisance parameter.
With plot = TRUE, a series of boxplots is drawn (cf. figure above).
Since the maximum sample size is obtained if the overall response rate is
estimated to be 0.5 in the sample size recalculation, this maximum
can occur under any true value of the nuisance parameter (except for 0 and 1),
albeit with very small probability.
For this reason, sample sizes that occur with a probability of less than
0.01% are ignored. This is not the case for continuous outcomes since there, the
sample size distributions are determined by simulation.
For continuous outcomes, i.e., the (shifted) \(t\)-test, the functional
content of blindrecalc is the same as in the binary case that
was presented in this example.
The only difference is that in the continuous case, the numbers are computed
by simulation.
Thus, the user can set the parameters iters, defining the number
of simulation iterations, and seed, the random seed for the simulation.
Development principles
The utilization of R’s object-oriented programming capabilities implies
that the example that was presented for the chi-squared test could
very similarly be applied to the Farrington-Manning test or the \(t\)-test.
Besides using S4 classes, the following development principles of
blindrecalc should be briefly described.
All calculations for binary outcomes are exact and require nested for-loops.
Since for-loops are known to be very slow in R, all computation-intensive functions
for the chi-squared test and the Farrington-Manning test are implemented in
C++ via the package (Eddelbuettel and François 2011) to speed up the calculations significantly.
blindrecalc is developed open-source on GitHub.com.
The entire source code can be found at
.
This allows anyone to contribute to blindrecalc and, furthermore,
provides maximal transparency.
To ensure a certain quality of the provided code, blindrecalc is checked
by unit tests using the package (Hester 2020).
The unit tests compare numbers for the sample size,
type 1 error rate, and power calculated with blindrecalc with numbers
from peer-reviewed publications and, furthermore, check the technical functionality
of the package such as vectorization and display of error messages.
Thus, the unit tests do not only monitor the technical accuracy
of the package’s results but also their content-related correctness.
The current version blindrecalc 0.1.3 achieves a code coverage of 100%,
i.e., each line of the source code is checked by at least one unit test.