---
title: "Pilot Trial Sample Size"
author: "Joel Ofoe Dadiboe"
format: 
  html:
    embed-resources: true
    toc: true
vignette: >
  %\VignetteIndexEntry{ssPilot}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---




# Introduction

This document develops and validates the R package `ssPilot` that implements the sample size calculations described by Whitehead et al. (2016).

```{r}
library(ssPilot)
```

The methods considered in this project are:

1. Upper confidence limit (UCL) approach
2. Non-central t-distribution (NCT) approach
3. Optimal pilot trial sample size based on UCL and NCT approaches

The final calculations will be implemented as functions in the
`ssPilot` R package.

# Standard Sample Size

For a two-arm trial with a continuous normally distributed outcome,
the standard sample size per treatment arm is

$$
n =
\frac{(r+1)(z_{1-\beta}+z_{1-\alpha/2})^2\sigma^2}
{rd^2}.
$$

where:

- $r$ is the allocation ratio,
- $\sigma$ is the population standard deviation,
- $d$ is the treatment difference,
- $\alpha$ is the Type I error rate,
- $1-\beta$ is the desired power.

## Example

Suppose that

$$
\sigma = 1,
$$

and

$$
d = 0.5.
$$

We use a two-sided significance level of 5% and 90% power.

```{r}
standard_sample_size(
  sd = 1,
  effect = 0.5,
  power = 0.90,
  alpha = 0.05
)
```


# Upper Confidence Limit Method

Whitehead et al. describe an approach in which the uncertainty in
the estimate of the population variance is incorporated into the
main-trial sample-size calculation.

The method uses an upper confidence limit for the variance estimated
from the pilot trial. This produces a more conservative estimate of
the variance than simply using the pilot estimate itself.

The upper confidence limit for the variance is given by

$$
s^2_{UCL}
=
\frac{k s^2}
{\chi^2_{1-X,k}},
$$

where $s^2$ is the variance estimated from the pilot trial, $k$ is
the degrees of freedom associated with the variance estimate, and
$X$ is the confidence level used for the upper confidence limit.

The quantity $\chi^2_{1-X,k}$ is the lower-tail chi-squared quantile
with $k$ degrees of freedom.

For a two-arm pilot trial with equal allocation and $m$ participants
per treatment arm, the degrees of freedom for the pooled variance
estimate are

$$
k = 2m - 2.
$$

Therefore, the degrees of freedom depend directly on the number of
participants included in the pilot trial.

For example, if the pilot trial contains 12 participants per
treatment arm, then

$$
k = 2(12)-2 = 22.
$$

We can calculate this directly in R.
```{r}
pilot_n <- 12
k <- 2 * pilot_n - 2
k
```



The next step is to calculate the upper confidence limit for the
variance.

Suppose that the standard deviation estimated from the pilot trial
is 1. The corresponding variance is therefore

$$
s^2 = 1^2 = 1.
$$

Whitehead et al. consider an 80% confidence level for the upper
confidence limit. Thus, we set

$$
X = 0.80.
$$

Using the UCL equation,

$$
s^2_{UCL}
=
\frac{k s^2}
{\chi^2_{1-X,k}},
$$

we can calculate the upper confidence limit for the variance in R.
```{r}
sd <- 1
variance <- sd^2
conf_level <- 0.80

variance_ucl <-
  k * variance /
  qchisq(1 - conf_level, df = k)

variance_ucl
```



The UCL calculation above produces an upper confidence limit for the
variance.

Because the main-trial sample-size equation is expressed in terms
of the standard deviation, we take the square root of the upper
confidence limit for the variance.

Thus,

$$
s_{UCL}
=
\sqrt{s^2_{UCL}}.
$$

We calculate this quantity in R as follows.
```{r}
sd_ucl <- sqrt(variance_ucl)

sd_ucl
```


The purpose of the UCL method is to incorporate the uncertainty in
the pilot estimate of the standard deviation into the sample-size
calculation for the main trial.

The upper confidence limit for the standard deviation is therefore
used in place of the original pilot estimate of the standard
deviation.

The standard sample-size equation is

$$
n_M =
\frac{(r+1)(z_{1-\beta}+z_{1-\alpha/2})^2s^2_{UCL}}
{rd^2}.
$$

Under the UCL approach, the variance estimate is replaced by the
upper confidence limit for the variance. Equivalently, we can use
$s_{UCL}$ as the standard deviation in our implementation of the
standard sample-size function.

For illustration, suppose that the standardized treatment
difference is

$$
\delta = \frac{d}{\sigma} = 0.50.
$$

We can represent this using

$$
s = 1
$$

and

$$
d = 0.50.
$$

We use 90% power, a two-sided 5% significance level, and equal
allocation between the two treatment arms.
```{r}
effect <- 0.50
power <- 0.90
alpha <- 0.05
allocation <- 1

main_n <- standard_sample_size(
  sd = sd_ucl,
  effect = effect,
  power = power,
  alpha = alpha,
  allocation = allocation
)

main_n
```


The previous calculations can be combined into a single sequence.

Starting with the pilot sample size and the pilot standard deviation,
we first calculate the degrees of freedom. We then calculate the
upper confidence limit for the variance, obtain the corresponding
upper confidence limit for the standard deviation, and finally use
that value in the main-trial sample-size calculation.

This sequence represents the computational structure that we will
later implement as the `ucl_sample_size()` function in the
`ssPilot` package.
```{r}
pilot_n <- 12
sd <- 1
conf_level <- 0.80
effect <- 0.50
power <- 0.90
alpha <- 0.05
allocation <- 1

# Degrees of freedom
k <- 2 * pilot_n - 2

# Upper confidence limit for the variance
variance_ucl <-
  k * sd^2 /
  qchisq(1 - conf_level, df = k)

# Upper confidence limit for the standard deviation
sd_ucl <- sqrt(variance_ucl)

# Main-trial sample size
main_n <- standard_sample_size(
  sd = sd_ucl,
  effect = effect,
  power = power,
  alpha = alpha,
  allocation = allocation
)

main_n
```

The calculations above demonstrate the complete UCL procedure.

The pilot sample size determines the degrees of freedom. The pilot
standard deviation is then used to obtain an upper confidence limit
for the variance. The square root of this quantity gives the upper
confidence limit for the standard deviation, which is subsequently
used to calculate the required main-trial sample size.

The next step is to translate this sequence into the
`ucl_sample_size()` function in the R package.

```{r}
ucl_sample_size(
  pilot_n = 12,
  sd = 1,
  effect = 0.50,
  power = 0.90,
  alpha = 0.05,
  conf_level = 0.80
)
```


<!-- # Non-Central t-Distribution Method -->

<!-- The second method considered by Whitehead et al. is based on the -->
<!-- non-central t-distribution. -->

<!-- The conventional sample-size calculation uses a normal approximation -->
<!-- to determine the required sample size for the main trial. The -->
<!-- non-central t-distribution approach instead calculates the power -->
<!-- using the exact distribution of the test statistic under the -->
<!-- alternative hypothesis. -->

<!-- This approach is particularly useful when the sample size is small, -->
<!-- because the normal approximation can be less accurate in small -->
<!-- samples. -->

<!-- The non-central t-distribution is characterized by two important -->
<!-- quantities: -->

<!-- - the degrees of freedom, and -->
<!-- - the non-centrality parameter. -->

<!-- For a two-arm trial with allocation ratio $r$, treatment difference -->
<!-- $d$, standard deviation $\sigma$, and sample size $n$ in the control -->
<!-- group, the non-centrality parameter is -->

<!-- $$ -->
<!-- \lambda = -->
<!-- \frac{d} -->
<!-- {\sigma\sqrt{1/n + 1/(rn)}}. -->
<!-- $$ -->

<!-- For equal allocation, where $r=1$, this simplifies to -->

<!-- $$ -->
<!-- \lambda = -->
<!-- \frac{d\sqrt{n}}{\sigma\sqrt{2}}. -->
<!-- $$ -->

<!-- For a two-arm trial with equal allocation, the degrees of freedom -->
<!-- for the two-sample t-test are -->

<!-- $$ -->
<!-- df = 2n - 2. -->
<!-- $$ -->

<!-- Thus, the sample size determines both the degrees of freedom and -->
<!-- the non-centrality parameter of the non-central t-distribution. -->

<!-- For example, if there are 50 participants per treatment arm, the -->
<!-- degrees of freedom are -->

<!-- $$ -->
<!-- df = 2(50)-2 = 98. -->
<!-- $$ -->

<!-- We can calculate this directly in R. -->
<!-- ```{r} -->
<!-- n <- 50 -->

<!-- df <- 2 * n - 2 -->

<!-- df -->
<!-- ``` -->


<!-- The next quantity we need is the non-centrality parameter. -->

<!-- For equal allocation, the non-centrality parameter is -->

<!-- $$ -->
<!-- \lambda = -->
<!-- \frac{d\sqrt{n}} -->
<!-- {\sigma\sqrt{2}}. -->
<!-- $$ -->

<!-- Suppose that the standardized treatment difference is 0.50. We can -->
<!-- represent this using -->

<!-- $$ -->
<!-- d = 0.50 -->
<!-- $$ -->

<!-- and -->

<!-- $$ -->
<!-- \sigma = 1. -->
<!-- $$ -->

<!-- If there are 50 participants per treatment arm, the non-centrality -->
<!-- parameter is therefore calculated using the treatment difference, -->
<!-- standard deviation, and sample size. -->
<!-- ```{r} -->
<!-- effect <- 0.50 -->

<!-- sd <- 1 -->

<!-- n <- 50 -->

<!-- ncp <- effect * sqrt(n) / (sd * sqrt(2)) -->

<!-- ncp -->
<!-- ``` -->


<!-- R provides the function `pt()` for calculating probabilities from -->
<!-- the t-distribution. -->

<!-- When the argument `ncp` is supplied, `pt()` calculates probabilities -->
<!-- from the non-central t-distribution. -->

<!-- For a two-sided hypothesis test, the rejection region contains both -->
<!-- the lower and upper tails of the t-distribution. Therefore, the -->
<!-- power is obtained by calculating the probability of falling in the -->
<!-- rejection region under the alternative hypothesis. -->

<!-- For a significance level of $\alpha$, the critical values are -->

<!-- $$ -->
<!-- -t_{1-\alpha/2,df} -->
<!-- $$ -->

<!-- and -->

<!-- $$ -->
<!-- t_{1-\alpha/2,df}. -->
<!-- $$ -->

<!-- For a two-sided test with $\alpha=0.05$, the upper critical value is -->

<!-- $$ -->
<!-- t_{1-\alpha/2,df} -->
<!-- = -->
<!-- t_{0.975,df}. -->
<!-- $$ -->

<!-- We can obtain this value in R using the `qt()` function. -->

<!-- The `qt()` function returns quantiles from the t-distribution. -->
<!-- ```{r} -->
<!-- alpha <- 0.05 -->

<!-- critical_value <- qt( -->
<!--   1 - alpha / 2, -->
<!--   df = df -->
<!-- ) -->

<!-- critical_value -->
<!-- ``` -->


<!-- The power of the two-sided test is the probability that the test -->
<!-- statistic falls in the rejection region when the alternative -->
<!-- hypothesis is true. -->

<!-- Using the non-central t-distribution, the power can therefore be -->
<!-- calculated as -->

<!-- $$ -->
<!-- P(T < -t_{1-\alpha/2,df}) -->
<!-- + -->
<!-- P(T > t_{1-\alpha/2,df}), -->
<!-- $$ -->

<!-- where $T$ follows a non-central t-distribution with the specified -->
<!-- degrees of freedom and non-centrality parameter. -->

<!-- In R, the two probabilities can be obtained using `pt()`. -->

<!-- The lower-tail probability is -->

<!-- $$ -->
<!-- P(T < -t_{1-\alpha/2,df}), -->
<!-- $$ -->

<!-- and the upper-tail probability is -->

<!-- $$ -->
<!-- P(T > t_{1-\alpha/2,df}). -->
<!-- $$ -->

<!-- ```{r} -->
<!-- lower_tail <- pt( -->
<!--   -critical_value, -->
<!--   df = df, -->
<!--   ncp = ncp -->
<!-- ) -->

<!-- upper_tail <- 1 - pt( -->
<!--   critical_value, -->
<!--   df = df, -->
<!--   ncp = ncp -->
<!-- ) -->

<!-- power_nct <- lower_tail + upper_tail -->

<!-- power_nct -->
<!-- ``` -->

<!-- The important difference between the standard and non-central -->
<!-- t-distribution approaches is that the standard approach directly -->
<!-- solves a normal-approximation sample-size formula, whereas the NCT -->
<!-- approach calculates the actual power associated with a proposed -->
<!-- sample size. -->

<!-- Therefore, to determine the required sample size using the NCT -->
<!-- method, we need to search over possible sample sizes until the -->
<!-- calculated power reaches the desired target power. -->

<!-- For example, if the desired power is 90%, we seek the smallest -->
<!-- sample size for which -->

<!-- $$ -->
<!-- Power(n) \geq 0.90. -->
<!-- $$ -->



<!-- We can evaluate the NCT power over a sequence of possible sample -->
<!-- sizes. -->

<!-- For each candidate sample size, we calculate: -->

<!-- 1. the degrees of freedom, -->
<!-- 2. the non-centrality parameter, -->
<!-- 3. the critical t-value, -->
<!-- 4. the lower-tail probability, -->
<!-- 5. the upper-tail probability, and -->
<!-- 6. the resulting power. -->

<!-- We then identify the smallest sample size for which the calculated -->
<!-- power is at least the desired power. -->
<!-- ```{r} -->
<!-- effect <- 0.50 -->
<!-- sd <- 1 -->
<!-- alpha <- 0.05 -->
<!-- target_power <- 0.90 -->

<!-- candidate_n <- 2:200 -->

<!-- nct_power <- sapply(candidate_n, function(n) { -->

<!--   df <- 2 * n - 2 -->

<!--   ncp <- effect * sqrt(n) / (sd * sqrt(2)) -->

<!--   critical_value <- qt( -->
<!--     1 - alpha / 2, -->
<!--     df = df -->
<!--   ) -->

<!--   lower_tail <- pt( -->
<!--     -critical_value, -->
<!--     df = df, -->
<!--     ncp = ncp -->
<!--   ) -->

<!--   upper_tail <- 1 - pt( -->
<!--     critical_value, -->
<!--     df = df, -->
<!--     ncp = ncp -->
<!--   ) -->

<!--   lower_tail + upper_tail -->
<!-- }) -->

<!-- required_n <- candidate_n[ -->
<!--   which(nct_power >= target_power)[1] -->
<!-- ] -->

<!-- required_n -->
<!-- ``` -->

<!-- The value returned above is the smallest number of participants per -->
<!-- treatment arm for which the calculated NCT power reaches the target -->
<!-- power of 90%. -->

<!-- We can also inspect the calculated power at this sample size. -->
<!-- ```{r} -->
<!-- nct_power[ -->
<!--   candidate_n == required_n -->
<!-- ] -->
<!-- ``` -->

<!-- We can also inspect the previous sample size to confirm that it does -->
<!-- not achieve the desired power. -->

<!-- ```{r} -->
<!-- nct_power[ -->
<!--   candidate_n == required_n - 1 -->
<!-- ] -->
<!-- ``` -->


<!-- ### Verification of the Required Sample Size -->

<!-- The NCT function identifies the smallest sample size for which the -->
<!-- target power is achieved. -->

<!-- For the example above, the required sample size is 86 participants -->
<!-- per treatment arm. -->

<!-- We can verify that this is the minimum by examining the power at -->
<!-- the preceding sample size, 85 participants per treatment arm. -->

<!-- If the implementation is correct, the power at 85 should be below -->
<!-- 90%. -->

<!-- ```{r} -->
<!-- n <- 85 -->

<!-- df <- 2 * n - 2 -->

<!-- ncp <- effect * sqrt(n) / (sd * sqrt(2)) -->

<!-- critical_value <- qt( -->
<!--   1 - alpha / 2, -->
<!--   df = df -->
<!-- ) -->

<!-- lower_tail <- pt( -->
<!--   -critical_value, -->
<!--   df = df, -->
<!--   ncp = ncp -->
<!-- ) -->

<!-- upper_tail <- 1 - pt( -->
<!--   critical_value, -->
<!--   df = df, -->
<!--   ncp = ncp -->
<!-- ) -->

<!-- power_85 <- lower_tail + upper_tail -->

<!-- power_85 -->
<!-- ``` -->
<!-- The power is then calculated again at 86 participants per treatment -->
<!-- arm. -->

<!-- Because 86 is the sample size returned by the NCT function, the -->
<!-- result should be at least 90%. -->
<!-- ```{r} -->
<!-- n <- 86 -->

<!-- df <- 2 * n - 2 -->

<!-- ncp <- effect * sqrt(n) / (sd * sqrt(2)) -->

<!-- critical_value <- qt( -->
<!--   1 - alpha / 2, -->
<!--   df = df -->
<!-- ) -->

<!-- lower_tail <- pt( -->
<!--   -critical_value, -->
<!--   df = df, -->
<!--   ncp = ncp -->
<!-- ) -->

<!-- upper_tail <- 1 - pt( -->
<!--   critical_value, -->
<!--   df = df, -->
<!--   ncp = ncp -->
<!-- ) -->

<!-- power_86 <- lower_tail + upper_tail -->

<!-- power_86 -->
<!-- ``` -->

<!-- The results confirm that 85 participants per arm do not achieve the -->
<!-- target power, whereas 86 participants per arm do. -->

<!-- Therefore, 86 is the smallest sample size per treatment arm that -->
<!-- achieves at least 90% power under this NCT calculation. -->

<!-- ```{r} -->
<!-- nct_sample_size( -->
<!--   pilot_n = 12, -->
<!--   sd = 1, -->
<!--   effect = 0.50, -->
<!--   power = 0.90, -->
<!--   alpha = 0.05 -->
<!-- ) -->
<!-- ``` -->


# Optimized UCL Pilot Sample Size

The UCL method requires a pilot study to estimate the population
standard deviation. Increasing the pilot sample size improves the
precision of the variance estimate, but it also increases the total
number of participants required.

The optimized approach considers this trade-off.

For each possible pilot sample size, we:

1. calculate the degrees of freedom,
2. calculate the upper confidence limit for the variance,
3. obtain the corresponding upper confidence limit for the standard
   deviation,
4. calculate the main-trial sample size using the UCL and NCT methodS, and
5. calculate the combined pilot and main-trial sample size.

The optimal pilot sample size is the value that minimizes the total
sample size.

The optimization is subject to a minimum pilot sample size of
10 participants per treatment group.

We begin by defining a sequence of possible pilot sample sizes.

For illustration, we consider pilot sample sizes from 10 through
50 participants per treatment arm.
```{r}
pilot_sizes <- 10:50

pilot_sizes
```


For each candidate pilot sample size, the degrees of freedom are

$$
k = 2m-2.
$$

The upper confidence limit for the variance is then

$$
s^2_{UCL}
=
\frac{k s^2}
{\chi^2_{1-X,k}}.
$$

We can therefore calculate the UCL standard deviation for every
candidate pilot sample size.
```{r}
sd <- 1
conf_level <- 0.80

variance_ucl <- sapply(pilot_sizes, function(m) {
  k <- 2 * m - 2
  k * sd^2 /
    qchisq(1 - conf_level, df = k)
})

sd_ucl <- sqrt(variance_ucl)
sd_ucl
```


Each upper confidence limit for the standard deviation is then used
to calculate the corresponding main-trial sample size.

Thus, every candidate pilot sample size produces a corresponding
main-trial sample size.
```{r}
effect <- 0.50
power <- 0.90
alpha <- 0.05

main_sizes <- sapply(sd_ucl, function(s) {
  standard_sample_size(
    sd = s,
    effect = effect,
    power = power,
    alpha = alpha,
    allocation = 1
  )
})

main_sizes
```


The objective of the optimization is to minimize the combined
number of participants in the pilot and main trials.

Therefore, for each candidate pilot size $m$, we calculate

$$
N_{\text{total}} = m+n_M,
$$

where $m$ is the pilot sample size per treatment arm and $n_M$ is
the resulting main-trial sample size per treatment arm.
```{r}
total_sizes <- pilot_sizes + main_sizes

total_sizes
```


The optimal pilot sample size is the candidate that produces the
smallest total sample size.
```{r}
optimal_index <- which.min(total_sizes)
optimal_pilot_n <- pilot_sizes[optimal_index]
optimal_main_n <- main_sizes[optimal_index]
optimal_total_n <- total_sizes[optimal_index]

optimal_pilot_n
optimal_main_n
optimal_total_n
```


The candidate pilot sample sizes and their corresponding main-trial
and total sample sizes can be displayed together in a table.
```{r}
optimization_results <- data.frame(
  pilot_n_per_arm = pilot_sizes,
  main_n_per_arm = main_sizes,
  total_n_per_arm = total_sizes
)

optimization_results
```



# Comparison of Methods

We can now compare the sample-size calculations implemented in the
package.

For the following example, we use a standard deviation of 1, a
treatment effect of 0.50, 90% power, and a two-sided significance
level of 0.05.

The UCL and UCL optimized methods use an 80% upper confidence limit.
```{r}
sd <- 1
effect <- 0.50
power <- 0.90
alpha <- 0.05
conf_level <- 0.80

standard_result <- standard_sample_size(
  sd = sd,
  effect = effect,
  power = power,
  alpha = alpha
)

ucl_result <- ucl_sample_size(
  pilot_n = 12, # as suggested by Julious (2005)
  sd = sd,
  effect = effect,
  power = power,
  alpha = alpha,
  conf_level = conf_level
)

optimized_ucl_result <- optimized_ucl_sample_size(
  sd = sd,
  effect = effect,
  power = power,
  alpha = alpha,
  conf_level = conf_level
)


```

<!-- nct_result <- nct_sample_size( -->
<!--   pilot_n = 16, -->
<!--   sd = sd, -->
<!--   effect = effect, -->
<!--   power = power, -->
<!--   alpha = alpha -->
<!-- ) -->

<!-- optimized_nct_result <- optimized_nct_sample_size( -->
<!--   sd = sd, -->
<!--   effect = effect, -->
<!--   power = power, -->
<!--   alpha = alpha -->
<!-- ) -->


The results from the four calculations can be summarized in a single
table.
```{r}
comparison <- data.frame(
  method = c(
    "Standard",
    "UCL",
    "Optimized UCL"
  ),
  main_n_per_arm = c(
    standard_result,
    ucl_result$main_n_per_arm,
    optimized_ucl_result$main_n_per_arm
  )
)

comparison
```

# References

Whitehead AL, Julious SA, Cooper CL, Campbell MJ. Estimating the sample size for a pilot randomised trial to minimise the overall trial sample size for the external pilot and main trial for a continuous outcome variable. Stat Methods Med Res. 2016 Jun;25(3):1057–73. doi:10.1177/0962280215588241

Browne RH. On the use of a pilot sample for sample size determination. Statistics in Medicine. 1995 Sep 15;14(17):1933–40. doi:10.1002/sim.4780141709

Julious SA, Owen RJ. Sample size calculations for clinical studies allowing for uncertainty about the variance. Pharmaceutical Statistics. 2006 Jan;5(1):29–37. doi:10.1002/pst.197

Sim J, Lewis M. The size of a pilot study for a clinical trial should be calculated in relation to considerations of precision and efficiency. Journal of Clinical Epidemiology. 2012 Mar;65(3):301–8. doi:10.1016/j.jclinepi.2011.07.011

Julious SA. Sample size of 12 per group rule of thumb for a pilot study. Pharmaceutical Statistics. 2005 Oct;4(4):287–91. doi:10.1002/pst.185