1. Table of Contents


This project explores different sample size and power calculations for mean comparison tests in clinical research using various helpful packages in R. The important factors to be assessed prior to determining the appropriate sample sizes were evaluated for the One-Sample, Two-Sample and Multiple-Sample One-Way ANOVA Pairwise designs. Power analyses were conducted to address the requirements across different study hypotheses including Tests of Equality, Non-Inferiority, Superiority and Equivalence. All results were consolidated in a Summary presented at the end of the document.

Mean comparison tests applied during clinical research refer to trials evaluated in terms of mean responses of certain primary study endpoints. The objectives of the intended clinical trials usually include the evaluation of intervention effects, demonstration of therapeutic equivalence or non-inferiority and establishment of superiority. The equations applied in this study (mostly contained from the book Sample Size Calculations in Clinical Research and implemented using the base package) attempt to provide calculations for the sample size and statistical power based on the study objectives and hypotheses of interest.

1.1 One-Sample Design (OSD)


One-Sample Design involves individual numeric response values from sampled subjects. In clinical research, responses could be the difference between matched pairs such as the pre-treatment and post-treatment responses or changes from baseline to endpoint within a treatment group. Responses are assumed to be independent and identically distributed normal random variables with zero mean and a given variance. Parameters needed for the computations are as follows:

[A] Mean Response : True mean response value obtained from sampled subjects

[B] Reference : Baseline value for comparison

[C] Standard Deviation : Data variability based from literature, pilot studies or preliminary experiments

[D] Alpha : Maximum allowed Type I error

[E] Beta : Maximum allowed Type II error

[F] Delta :
     [F.1] Minimum testing margin to establish that the mean response is much better than reference
     [F.2] Minimum testing margin to establish that the mean response is not much worse as the reference
     [F.3] Minimum testing margin to establish that the mean response is no better and no worse than the reference

A case scenario was provided as follows to demonstrate sample size computation and power analysis as applied on different types of hypotheses:

[Case Scenario] A clinical research study is being pursued concerning osteoporosis in post-menopausal women. Osteoporosis and osteopenia (or decreased bone mass) most commonly develop in post-menopausal women. The consequences of osteoporosis are vertebral crush fractures and hip fractures. The diagnosis of osteoporosis is made when vertebral bone density is more than 10% below what is expected for sex, age, height, weight, and race. Usually, bone density is reported in terms of standard deviation (SD) from mean values. The World Health Organization (WHO) defines osteopenia as bone density value greater than one SD below peak bone mass levels in young women and osteoporosis as a value of greater than 2.5 SD below the same measurement scale. In medical practice, most clinicians suggest therapeutic intervention should be begun in patients with osteopenia to prevent progression to osteoporosis.

1.1.1 Test for Equality (OSD_EQUALITY)


[Research Question] Suppose that the mean bone density before the treatment is 1.5 SD. After treatment, the mean value is expected to range from 2.0 to 4.0 SD with increments of 0.2. For prevention of progression from osteopenia to osteoporosis, the study aims to demonstrate that the mean bone density post-treatment and pre-treatment are different. Given a level of significance of 0.05, evaluate the sample size range for a test of equality to achieve statistical test powers equal to 70%, 80% and 90%.

[Null Hypothesis] Post-treatment mean bone density - Pre-treatment mean bone density equals 0

[Alternative Hypothesis] Post-treatment mean bone density - Pre-treatment mean bone density does not equal 0

[Hypothesis Test] Rejection of the null hypothesis implies difference between the mean bone density values post-treatment and pre-treatment

[A] Mean Response : 2.0 to 4.0 SD at 0.20 intervals

[B] Reference : 1.5 SD

[C] Standard Deviation : Not explicitly given but already factored into the mean values and will cancel out during computation

[D] Alpha : 0.05

[E] Beta : 0.30, 0.20 and 0.10 corresponding to 70%, 80% and 90% statistical powers

[F] Delta : 0 testing margin for the Test of Equality hypothesis

Code Chunk | Output
##################################
# Loading R libraries
##################################
library(moments)
library(car)
library(multcomp)
library(effects)
library(psych)
library(ggplot2)
library(dplyr)
library(ggpubr)
library(rstatix)
library(ggfortify)
library(trend)

##################################
# Defining a range of possible hypothesized values for mu_Response
##################################
mu_Response=seq(2.0,4.0,0.2)

##################################
# Defining a fixed hypothesized value for mu_Reference
##################################
mu_Reference=1.5

##################################
# Defining a fixed hypothesized value for the standard deviation
##################################
sd=1

##################################
# Defining a fixed hypothesized value for the Type I error
##################################
alpha=0.05

##################################
# Defining a range of possible hypothesized values for the Type II error
##################################
beta=c(0.30,0.20,0.10)

##################################
# Computing the range of samples sizes based on different levels of Type II errors
##################################
beta_1=beta[1]
(n_Beta30=(sd*(qnorm(1-alpha/2)+qnorm(1-beta_1))/(mu_Response-mu_Reference))^2)
##  [1] 24.6882678 12.5960550  7.6198357  5.1008818  3.6521106  2.7431409
##  [7]  2.1356633  1.7097138  1.3995617  1.1667423  0.9875307
ceiling(n_Beta30)
##  [1] 25 13  8  6  4  3  3  2  2  2  1
z_Beta30=(mu_Response-mu_Reference)/(sd/sqrt(n_Beta30))
(Power_Beta30=pnorm(z_Beta30-qnorm(1-alpha/2))+pnorm(-z_Beta30-qnorm(1-alpha/2)))
##  [1] 0.7000044 0.7000044 0.7000044 0.7000044 0.7000044 0.7000044 0.7000044
##  [8] 0.7000044 0.7000044 0.7000044 0.7000044
beta_2=beta[2]
(n_Beta20=(sd*(qnorm(1-alpha/2)+qnorm(1-beta_2))/(mu_Response-mu_Reference))^2)
##  [1] 31.395519 16.018122  9.689975  6.486677  4.644308  3.488391  2.715875
##  [8]  2.174205  1.779791  1.483720  1.255821
ceiling(n_Beta20)
##  [1] 32 17 10  7  5  4  3  3  2  2  2
z_Beta20=(mu_Response-mu_Reference)/(sd/sqrt(n_Beta20))
(Power_Beta20=pnorm(z_Beta20-qnorm(1-alpha/2))+pnorm(-z_Beta20-qnorm(1-alpha/2)))
##  [1] 0.800001 0.800001 0.800001 0.800001 0.800001 0.800001 0.800001 0.800001
##  [9] 0.800001 0.800001 0.800001
beta_3=beta[3]
(n_Beta10=(sd*(qnorm(1-alpha/2)+qnorm(1-beta_3))/(mu_Response-mu_Reference))^2)
##  [1] 42.029692 21.443721 12.972127  8.683821  6.217410  4.669966  3.635787
##  [8]  2.910644  2.382636  1.986280  1.681188
ceiling(n_Beta10)
##  [1] 43 22 13  9  7  5  4  3  3  2  2
z_Beta10=(mu_Response-mu_Reference)/(sd/sqrt(n_Beta10))
(Power_Beta10=pnorm(z_Beta10-qnorm(1-alpha/2))+pnorm(-z_Beta10-qnorm(1-alpha/2)))
##  [1] 0.9000001 0.9000001 0.9000001 0.9000001 0.9000001 0.9000001 0.9000001
##  [8] 0.9000001 0.9000001 0.9000001 0.9000001
OneSampleDesign.Equality <- as.data.frame(cbind(mu_Response,n_Beta30,n_Beta20,n_Beta10))
names(OneSampleDesign.Equality) <- c("Posttreatment.Mean.Bone.Density",
                                 "Power=70%",
                                 "Power=80%",
                                 "Power=90%")

##################################
# Restructuring the data
##################################
OneSampleDesign.Equality.Reshaped <- gather(OneSampleDesign.Equality,
                                            "Power=70%","Power=80%","Power=90%",
                                            key="Power",
                                            value="Sample.Size")
OneSampleDesign.Equality.Reshaped$Label <- rep("OSD_EQUALITY",nrow(OneSampleDesign.Equality.Reshaped))

OneSampleDesign.Equality.Reshaped
##    Posttreatment.Mean.Bone.Density     Power Sample.Size        Label
## 1                              2.0 Power=70%  24.6882678 OSD_EQUALITY
## 2                              2.2 Power=70%  12.5960550 OSD_EQUALITY
## 3                              2.4 Power=70%   7.6198357 OSD_EQUALITY
## 4                              2.6 Power=70%   5.1008818 OSD_EQUALITY
## 5                              2.8 Power=70%   3.6521106 OSD_EQUALITY
## 6                              3.0 Power=70%   2.7431409 OSD_EQUALITY
## 7                              3.2 Power=70%   2.1356633 OSD_EQUALITY
## 8                              3.4 Power=70%   1.7097138 OSD_EQUALITY
## 9                              3.6 Power=70%   1.3995617 OSD_EQUALITY
## 10                             3.8 Power=70%   1.1667423 OSD_EQUALITY
## 11                             4.0 Power=70%   0.9875307 OSD_EQUALITY
## 12                             2.0 Power=80%  31.3955189 OSD_EQUALITY
## 13                             2.2 Power=80%  16.0181219 OSD_EQUALITY
## 14                             2.4 Power=80%   9.6899750 OSD_EQUALITY
## 15                             2.6 Power=80%   6.4866775 OSD_EQUALITY
## 16                             2.8 Power=80%   4.6443075 OSD_EQUALITY
## 17                             3.0 Power=80%   3.4883910 OSD_EQUALITY
## 18                             3.2 Power=80%   2.7158753 OSD_EQUALITY
## 19                             3.4 Power=80%   2.1742049 OSD_EQUALITY
## 20                             3.6 Power=80%   1.7797913 OSD_EQUALITY
## 21                             3.8 Power=80%   1.4837202 OSD_EQUALITY
## 22                             4.0 Power=80%   1.2558208 OSD_EQUALITY
## 23                             2.0 Power=90%  42.0296922 OSD_EQUALITY
## 24                             2.2 Power=90%  21.4437205 OSD_EQUALITY
## 25                             2.4 Power=90%  12.9721272 OSD_EQUALITY
## 26                             2.6 Power=90%   8.6838207 OSD_EQUALITY
## 27                             2.8 Power=90%   6.2174101 OSD_EQUALITY
## 28                             3.0 Power=90%   4.6699658 OSD_EQUALITY
## 29                             3.2 Power=90%   3.6357865 OSD_EQUALITY
## 30                             3.4 Power=90%   2.9106435 OSD_EQUALITY
## 31                             3.6 Power=90%   2.3826356 OSD_EQUALITY
## 32                             3.8 Power=90%   1.9862804 OSD_EQUALITY
## 33                             4.0 Power=90%   1.6811877 OSD_EQUALITY
##################################
# Plotting the sample sizes and power curves
##################################
(OneSampleDesign.Equality.PowerAnalysis <- ggplot(OneSampleDesign.Equality.Reshaped,
                                                      aes(x=Posttreatment.Mean.Bone.Density,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Post-Treatment Mean Bone Density",
                      limits=c(2,4),
                      breaks=seq(2,4,
                                 by=0.2)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50)) +
  ggtitle("One-Sample Design : Test for Equality"))

1.1.2 Test for Non-Inferiority (OSD_NONINFERIORITY)


[Research Question] Suppose that the mean bone density before the treatment is 1.5 SD. After treatment, the mean value is expected to range from 2.0 to 4.0 SD with increments of 0.2. For prevention of progression from osteopenia to osteoporosis, the study aims to demonstrate that the mean bone density post-treatment is not much worse than pre-treatment by a clinically meaningful difference equal to −0.5 SD. Given a level of significance of 0.05, evaluate the sample size range for a test of non-inferiority to achieve statistical test powers equal to 70%, 80% and 90%.

[Null Hypothesis] Post-treatment mean bone density - Pre-treatment mean bone density <= Clinically meaningful difference

[Alternative Hypothesis] Post-treatment mean bone density - Pre-treatment mean bone density > Clinically meaningful difference

[Hypothesis Test] Rejection of the null hypothesis implies non-inferiority of the mean bone density treatment as compared to non-treatment, given the delta

[A] Mean Response : 2.0 to 4.0 SD at 0.20 intervals

[B] Reference : 1.5 SD

[C] Standard Deviation : Not explicitly given but already factored into the mean values and will cancel out during computation

[D] Alpha : 0.05

[E] Beta : 0.30, 0.20 and 0.10 corresponding to 70%, 80% and 90% statistical powers

[F] Delta : -0.50 SD testing margin for the Test of Non-Inferiority hypothesis

Code Chunk | Output
##################################
# Defining a range of possible hypothesized values for mu_Response
##################################
mu_Response=seq(2.0,4.0,0.2)

##################################
# Defining a fixed hypothesized value for mu_Reference
##################################
mu_Reference=1.5

##################################
# Defining a fixed hypothesized value for the standard deviation
##################################
sd=1

##################################
# Defining a fixed hypothesized value for the Type I error
##################################
alpha=0.05

##################################
# Defining a range of possible hypothesized values for the Type II error
##################################
beta=c(0.30,0.20,0.10)

##################################
# Defining the non-inferiority delta
##################################
delta=-0.50

##################################
# Computing the range of samples sizes based on different levels of Type II errors
##################################
beta_1=beta[1]
(n_Beta30=(sd*(qnorm(1-alpha)+qnorm(1-beta_1))/(mu_Response-mu_Reference-delta))^2)
##  [1] 4.7056635 3.2678219 2.4008487 1.8381498 1.4523653 1.1764159 0.9722445
##  [8] 0.8169555 0.6961041 0.6002122 0.5228515
ceiling(n_Beta30)
##  [1] 5 4 3 2 2 2 1 1 1 1 1
z_Beta30=(mu_Response-mu_Reference-delta)/(sd/sqrt(n_Beta30))
(Power_Beta30=pnorm(z_Beta30-qnorm(1-alpha))+pnorm(-z_Beta30-qnorm(1-alpha)))
##  [1] 0.7000683 0.7000683 0.7000683 0.7000683 0.7000683 0.7000683 0.7000683
##  [8] 0.7000683 0.7000683 0.7000683 0.7000683
beta_2=beta[2]
(n_Beta20=(sd*(qnorm(1-alpha)+qnorm(1-beta_2))/(mu_Response-mu_Reference-delta))^2)
##  [1] 6.1825572 4.2934425 3.1543659 2.4150614 1.9081967 1.5456393 1.2773879
##  [8] 1.0733606 0.9145795 0.7885915 0.6869508
ceiling(n_Beta20)
##  [1] 7 5 4 3 2 2 2 2 1 1 1
z_Beta20=(mu_Response-mu_Reference-delta)/(sd/sqrt(n_Beta20))
(Power_Beta20=pnorm(z_Beta20-qnorm(1-alpha))+pnorm(-z_Beta20-qnorm(1-alpha)))
##  [1] 0.800018 0.800018 0.800018 0.800018 0.800018 0.800018 0.800018 0.800018
##  [9] 0.800018 0.800018 0.800018
beta_3=beta[3]
(n_Beta10=(sd*(qnorm(1-alpha)+qnorm(1-beta_3))/(mu_Response-mu_Reference-delta))^2)
##  [1] 8.5638474 5.9471162 4.3693099 3.3452529 2.6431628 2.1409618 1.7693899
##  [8] 1.4867791 1.2668413 1.0923275 0.9515386
ceiling(n_Beta10)
##  [1] 9 6 5 4 3 3 2 2 2 2 1
z_Beta10=(mu_Response-mu_Reference-delta)/(sd/sqrt(n_Beta10))
(Power_Beta10=pnorm(z_Beta10-qnorm(1-alpha))+pnorm(-z_Beta10-qnorm(1-alpha)))
##  [1] 0.9000024 0.9000024 0.9000024 0.9000024 0.9000024 0.9000024 0.9000024
##  [8] 0.9000024 0.9000024 0.9000024 0.9000024
OneSampleDesign.NonInferiority <- as.data.frame(cbind(mu_Response,n_Beta30,n_Beta20,n_Beta10))
names(OneSampleDesign.NonInferiority) <- c("Posttreatment.Mean.Bone.Density",
                                 "Power=70%",
                                 "Power=80%",
                                 "Power=90%")

##################################
# Restructuring the data
##################################
OneSampleDesign.NonInferiority.Reshaped <- gather(OneSampleDesign.NonInferiority,
                                            "Power=70%","Power=80%","Power=90%",
                                            key="Power",
                                            value="Sample.Size")
OneSampleDesign.NonInferiority.Reshaped$Label <- rep("OSD_NONINFERIORITY",nrow(OneSampleDesign.NonInferiority.Reshaped))

OneSampleDesign.NonInferiority.Reshaped
##    Posttreatment.Mean.Bone.Density     Power Sample.Size              Label
## 1                              2.0 Power=70%   4.7056635 OSD_NONINFERIORITY
## 2                              2.2 Power=70%   3.2678219 OSD_NONINFERIORITY
## 3                              2.4 Power=70%   2.4008487 OSD_NONINFERIORITY
## 4                              2.6 Power=70%   1.8381498 OSD_NONINFERIORITY
## 5                              2.8 Power=70%   1.4523653 OSD_NONINFERIORITY
## 6                              3.0 Power=70%   1.1764159 OSD_NONINFERIORITY
## 7                              3.2 Power=70%   0.9722445 OSD_NONINFERIORITY
## 8                              3.4 Power=70%   0.8169555 OSD_NONINFERIORITY
## 9                              3.6 Power=70%   0.6961041 OSD_NONINFERIORITY
## 10                             3.8 Power=70%   0.6002122 OSD_NONINFERIORITY
## 11                             4.0 Power=70%   0.5228515 OSD_NONINFERIORITY
## 12                             2.0 Power=80%   6.1825572 OSD_NONINFERIORITY
## 13                             2.2 Power=80%   4.2934425 OSD_NONINFERIORITY
## 14                             2.4 Power=80%   3.1543659 OSD_NONINFERIORITY
## 15                             2.6 Power=80%   2.4150614 OSD_NONINFERIORITY
## 16                             2.8 Power=80%   1.9081967 OSD_NONINFERIORITY
## 17                             3.0 Power=80%   1.5456393 OSD_NONINFERIORITY
## 18                             3.2 Power=80%   1.2773879 OSD_NONINFERIORITY
## 19                             3.4 Power=80%   1.0733606 OSD_NONINFERIORITY
## 20                             3.6 Power=80%   0.9145795 OSD_NONINFERIORITY
## 21                             3.8 Power=80%   0.7885915 OSD_NONINFERIORITY
## 22                             4.0 Power=80%   0.6869508 OSD_NONINFERIORITY
## 23                             2.0 Power=90%   8.5638474 OSD_NONINFERIORITY
## 24                             2.2 Power=90%   5.9471162 OSD_NONINFERIORITY
## 25                             2.4 Power=90%   4.3693099 OSD_NONINFERIORITY
## 26                             2.6 Power=90%   3.3452529 OSD_NONINFERIORITY
## 27                             2.8 Power=90%   2.6431628 OSD_NONINFERIORITY
## 28                             3.0 Power=90%   2.1409618 OSD_NONINFERIORITY
## 29                             3.2 Power=90%   1.7693899 OSD_NONINFERIORITY
## 30                             3.4 Power=90%   1.4867791 OSD_NONINFERIORITY
## 31                             3.6 Power=90%   1.2668413 OSD_NONINFERIORITY
## 32                             3.8 Power=90%   1.0923275 OSD_NONINFERIORITY
## 33                             4.0 Power=90%   0.9515386 OSD_NONINFERIORITY
##################################
# Plotting the sample sizes and power curves
##################################
(OneSampleDesign.NonInferiority.PowerAnalysis <-ggplot(OneSampleDesign.NonInferiority.Reshaped,
                                                      aes(x=Posttreatment.Mean.Bone.Density,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Post-Treatment Mean Bone Density",
                      limits=c(2,4),
                      breaks=seq(2,4,
                                 by=0.2)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50)) +
  ggtitle("One-Sample Design : Test for Non-Inferiority"))

1.1.3 Test for Superiority (OSD_SUPERIORITY)


[Research Question] Suppose that the mean bone density before the treatment is 1.5 SD. After treatment, the mean value is expected to range from 2.0 to 4.0 SD with increments of 0.2. For prevention of progression from osteopenia to osteoporosis, the study aims to demonstrate that the mean bone density post-treatment is much better than pre-treatment by a clinically meaningful difference equal to +0.10 SD. Given a level of significance of 0.05, evaluate the sample size range for a test of superiority to achieve statistical test powers equal to 70%, 80% and 90%.

[Null Hypothesis] Post-treatment mean bone density - Pre-treatment mean bone density <= Clinically meaningful difference

[Alternative Hypothesis] Post-treatment mean bone density - Pre-treatment mean bone density > Clinically meaningful difference

[Hypothesis Test] Rejection of the null hypothesis implies non-inferiority of the mean bone density treatment as compared to non-treatment, given the delta

[A] Mean Response : 2.0 to 4.0 SD at 0.20 intervals

[B] Reference : 1.5 SD

[C] Standard Deviation : Not explicitly given but already factored into the mean values and will cancel out during computation

[D] Alpha : 0.05

[E] Beta : 0.30, 0.20 and 0.10 corresponding to 70%, 80% and 90% statistical powers

[F] Delta : +0.10 SD testing margin for the Test of Superiority hypothesis

Code Chunk | Output
##################################
# Defining a range of possible hypothesized values for mu_Response
##################################
mu_Response=seq(2.0,4.0,0.2)

##################################
# Defining a fixed hypothesized value for mu_Reference
##################################
mu_Reference=1.5

##################################
# Defining a fixed hypothesized value for the standard deviation
##################################
sd=1

##################################
# Defining a fixed hypothesized value for the Type I error
##################################
alpha=0.05

##################################
# Defining a range of possible hypothesized values for the Type II error
##################################
beta=c(0.30,0.20,0.10)

##################################
# Defining the superiority delta
##################################
delta= 0.10

##################################
# Computing the range of samples sizes based on different levels of Type II errors
##################################
beta_1=beta[1]
(n_Beta30=(sd*(qnorm(1-alpha)+qnorm(1-beta_1))/(mu_Response-mu_Reference-delta))^2)
##  [1] 29.4103970 13.0712876  7.3525993  4.7056635  3.2678219  2.4008487
##  [7]  1.8381498  1.4523653  1.1764159  0.9722445  0.8169555
ceiling(n_Beta30)
##  [1] 30 14  8  5  4  3  2  2  2  1  1
z_Beta30=(mu_Response-mu_Reference-delta)/(sd/sqrt(n_Beta30))
(Power_Beta30=pnorm(z_Beta30-qnorm(1-alpha))+pnorm(-z_Beta30-qnorm(1-alpha)))
##  [1] 0.7000683 0.7000683 0.7000683 0.7000683 0.7000683 0.7000683 0.7000683
##  [8] 0.7000683 0.7000683 0.7000683 0.7000683
beta_2=beta[2]
(n_Beta20=(sd*(qnorm(1-alpha)+qnorm(1-beta_2))/(mu_Response-mu_Reference-delta))^2)
##  [1] 38.640983 17.173770  9.660246  6.182557  4.293443  3.154366  2.415061
##  [8]  1.908197  1.545639  1.277388  1.073361
ceiling(n_Beta20)
##  [1] 39 18 10  7  5  4  3  2  2  2  2
z_Beta20=(mu_Response-mu_Reference-delta)/(sd/sqrt(n_Beta20))
(Power_Beta20=pnorm(z_Beta20-qnorm(1-alpha))+pnorm(-z_Beta20-qnorm(1-alpha)))
##  [1] 0.800018 0.800018 0.800018 0.800018 0.800018 0.800018 0.800018 0.800018
##  [9] 0.800018 0.800018 0.800018
beta_3=beta[3]
(n_Beta10=(sd*(qnorm(1-alpha)+qnorm(1-beta_3))/(mu_Response-mu_Reference-delta))^2)
##  [1] 53.524046 23.788465 13.381011  8.563847  5.947116  4.369310  3.345253
##  [8]  2.643163  2.140962  1.769390  1.486779
ceiling(n_Beta10)
##  [1] 54 24 14  9  6  5  4  3  3  2  2
z_Beta10=(mu_Response-mu_Reference-delta)/(sd/sqrt(n_Beta10))
(Power_Beta10=pnorm(z_Beta10-qnorm(1-alpha))+pnorm(-z_Beta10-qnorm(1-alpha)))
##  [1] 0.9000024 0.9000024 0.9000024 0.9000024 0.9000024 0.9000024 0.9000024
##  [8] 0.9000024 0.9000024 0.9000024 0.9000024
OneSampleDesign.Superiority <- as.data.frame(cbind(mu_Response,n_Beta30,n_Beta20,n_Beta10))
names(OneSampleDesign.Superiority) <- c("Posttreatment.Mean.Bone.Density",
                                 "Power=70%",
                                 "Power=80%",
                                 "Power=90%")

##################################
# Restructuring the data
##################################
OneSampleDesign.Superiority.Reshaped <- gather(OneSampleDesign.Superiority,
                                            "Power=70%","Power=80%","Power=90%",
                                            key="Power",
                                            value="Sample.Size")
OneSampleDesign.Superiority.Reshaped$Label <- rep("OSD_SUPERIORITY",nrow(OneSampleDesign.Superiority.Reshaped))

OneSampleDesign.Superiority.Reshaped
##    Posttreatment.Mean.Bone.Density     Power Sample.Size           Label
## 1                              2.0 Power=70%  29.4103970 OSD_SUPERIORITY
## 2                              2.2 Power=70%  13.0712876 OSD_SUPERIORITY
## 3                              2.4 Power=70%   7.3525993 OSD_SUPERIORITY
## 4                              2.6 Power=70%   4.7056635 OSD_SUPERIORITY
## 5                              2.8 Power=70%   3.2678219 OSD_SUPERIORITY
## 6                              3.0 Power=70%   2.4008487 OSD_SUPERIORITY
## 7                              3.2 Power=70%   1.8381498 OSD_SUPERIORITY
## 8                              3.4 Power=70%   1.4523653 OSD_SUPERIORITY
## 9                              3.6 Power=70%   1.1764159 OSD_SUPERIORITY
## 10                             3.8 Power=70%   0.9722445 OSD_SUPERIORITY
## 11                             4.0 Power=70%   0.8169555 OSD_SUPERIORITY
## 12                             2.0 Power=80%  38.6409827 OSD_SUPERIORITY
## 13                             2.2 Power=80%  17.1737701 OSD_SUPERIORITY
## 14                             2.4 Power=80%   9.6602457 OSD_SUPERIORITY
## 15                             2.6 Power=80%   6.1825572 OSD_SUPERIORITY
## 16                             2.8 Power=80%   4.2934425 OSD_SUPERIORITY
## 17                             3.0 Power=80%   3.1543659 OSD_SUPERIORITY
## 18                             3.2 Power=80%   2.4150614 OSD_SUPERIORITY
## 19                             3.4 Power=80%   1.9081967 OSD_SUPERIORITY
## 20                             3.6 Power=80%   1.5456393 OSD_SUPERIORITY
## 21                             3.8 Power=80%   1.2773879 OSD_SUPERIORITY
## 22                             4.0 Power=80%   1.0733606 OSD_SUPERIORITY
## 23                             2.0 Power=90%  53.5240459 OSD_SUPERIORITY
## 24                             2.2 Power=90%  23.7884649 OSD_SUPERIORITY
## 25                             2.4 Power=90%  13.3810115 OSD_SUPERIORITY
## 26                             2.6 Power=90%   8.5638474 OSD_SUPERIORITY
## 27                             2.8 Power=90%   5.9471162 OSD_SUPERIORITY
## 28                             3.0 Power=90%   4.3693099 OSD_SUPERIORITY
## 29                             3.2 Power=90%   3.3452529 OSD_SUPERIORITY
## 30                             3.4 Power=90%   2.6431628 OSD_SUPERIORITY
## 31                             3.6 Power=90%   2.1409618 OSD_SUPERIORITY
## 32                             3.8 Power=90%   1.7693899 OSD_SUPERIORITY
## 33                             4.0 Power=90%   1.4867791 OSD_SUPERIORITY
##################################
# Plotting the sample sizes and power curves
##################################
(OneSampleDesign.Superiority.PowerAnalysis <- ggplot(OneSampleDesign.Superiority.Reshaped,
                                                      aes(x=Posttreatment.Mean.Bone.Density,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Post-Treatment Mean Bone Density",
                      limits=c(2,4),
                      breaks=seq(2,4,
                                 by=0.2)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50)) +
  ggtitle("One-Sample Design : Test for Superiority"))

1.1.4 Test for Equivalence (OSD_EQUIVALENCE)


[Research Question] Suppose that the mean bone density before the treatment is 1.5 SD. After treatment, the mean value is expected to range from 2.0 to 4.0 SD with increments of 0.2. For a special investigation, the study aims to demonstrate that the mean bone density post-treatment is not much worse or better than pre-treatment by a clinically meaningful difference equal to +0.10 SD. Given a level of significance of 0.05, evaluate the sample size range for a test of equivalence to achieve statistical test powers equal to 70%, 80% and 90%.

[Null Hypothesis] Absolute value of (Post-treatment mean bone density - Pre-treatment mean bone density) >= Clinically meaningful difference

[Alternative Hypothesis] Absolute value of (Post-treatment mean bone density - Pre-treatment) mean bone density < Clinically meaningful difference

[Hypothesis Test] Rejection of the null hypothesis implies equivalence of the mean bone density treatment as compared to non-treatment, given the delta

[A] Mean Response : 2.0 to 4.0 SD at 0.20 intervals

[B] Reference : 1.5 SD

[C] Standard Deviation : Not explicitly given but already factored into the mean values and will cancel out during computation

[D] Alpha : 0.05

[E] Beta : 0.30, 0.20 and 0.10 corresponding to 70%, 80% and 90% statistical powers

[F] Delta : +0.10 SD testing margin for the Test of Equivalence hypothesis

Code Chunk | Output
##################################
# Defining a range of possible hypothesized values for mu_Response
##################################
mu_Response=seq(2.0,4.0,0.2)

##################################
# Defining a fixed hypothesized value for mu_Reference
##################################
mu_Reference=1.5

##################################
# Defining a fixed hypothesized value for the standard deviation
##################################
sd=1

##################################
# Defining a fixed hypothesized value for the Type I error
##################################
alpha=0.05

##################################
# Defining a range of possible hypothesized values for the Type II error
##################################
beta=c(0.30,0.20,0.10)

##################################
# Defining the equivalence delta
##################################
delta= 0.10

##################################
# Computing the range of samples sizes based on different levels of Type II errors
##################################
beta_1=beta[1]
(n_Beta30=(sd*(qnorm(1-alpha)+qnorm(1-beta_1/2))/(delta-abs(mu_Response-mu_Reference)))^2)
##  [1] 44.933125 19.970278 11.233281  7.189300  4.992569  3.668010  2.808320
##  [8]  2.218920  1.797325  1.485393  1.248142
ceiling(n_Beta30)
##  [1] 45 20 12  8  5  4  3  3  2  2  2
z_Beta30=(abs(mu_Response-mu_Reference)-delta)/(sd/sqrt(n_Beta30))
(Power_Beta30=2*(pnorm(z_Beta30-qnorm(1-alpha))+pnorm(-z_Beta30-qnorm(1-alpha)))-1)
##  [1] 0.7000152 0.7000152 0.7000152 0.7000152 0.7000152 0.7000152 0.7000152
##  [8] 0.7000152 0.7000152 0.7000152 0.7000152
beta_2=beta[2]
(n_Beta20=(sd*(qnorm(1-alpha)+qnorm(1-beta_2/2))/(delta-abs(mu_Response-mu_Reference)))^2)
##  [1] 53.524046 23.788465 13.381011  8.563847  5.947116  4.369310  3.345253
##  [8]  2.643163  2.140962  1.769390  1.486779
ceiling(n_Beta20)
##  [1] 54 24 14  9  6  5  4  3  3  2  2
z_Beta20=(abs(mu_Response-mu_Reference)-delta)/(sd/sqrt(n_Beta20))
(Power_Beta20=2*(pnorm(z_Beta20-qnorm(1-alpha))+pnorm(-z_Beta20-qnorm(1-alpha)))-1)
##  [1] 0.8000048 0.8000048 0.8000048 0.8000048 0.8000048 0.8000048 0.8000048
##  [8] 0.8000048 0.8000048 0.8000048 0.8000048
(Power_Beta20=2*(pnorm(z_Beta20-qnorm(alpha)))-1)
##  [1] 0.9999952 0.9999952 0.9999952 0.9999952 0.9999952 0.9999952 0.9999952
##  [8] 0.9999952 0.9999952 0.9999952 0.9999952
beta_3=beta[3]
(n_Beta10=(sd*(qnorm(1-alpha)+qnorm(1-beta_3/2))/(delta-abs(mu_Response-mu_Reference)))^2)
##  [1] 67.638586 30.061594 16.909647 10.822174  7.515398  5.521517  4.227412
##  [8]  3.340177  2.705543  2.235986  1.878850
ceiling(n_Beta10)
##  [1] 68 31 17 11  8  6  5  4  3  3  2
z_Beta10=(abs(mu_Response-mu_Reference)-delta)/(sd/sqrt(n_Beta10))
(Power_Beta10=2*(pnorm(z_Beta10-qnorm(1-alpha))+pnorm(-z_Beta10-qnorm(1-alpha)))-1)
##  [1] 0.9000008 0.9000008 0.9000008 0.9000008 0.9000008 0.9000008 0.9000008
##  [8] 0.9000008 0.9000008 0.9000008 0.9000008
OneSampleDesign.Equivalence <- as.data.frame(cbind(mu_Response,n_Beta30,n_Beta20,n_Beta10))
names(OneSampleDesign.Equivalence) <- c("Posttreatment.Mean.Bone.Density",
                                 "Power=70%",
                                 "Power=80%",
                                 "Power=90%")

##################################
# Restructuring the data
##################################
OneSampleDesign.Equivalence.Reshaped <- gather(OneSampleDesign.Equivalence,
                                            "Power=70%","Power=80%","Power=90%",
                                            key="Power",
                                            value="Sample.Size")
OneSampleDesign.Equivalence.Reshaped$Label <- rep("OSD_EQUIVALENCE",nrow(OneSampleDesign.Equivalence.Reshaped))

OneSampleDesign.Equivalence.Reshaped
##    Posttreatment.Mean.Bone.Density     Power Sample.Size           Label
## 1                              2.0 Power=70%   44.933125 OSD_EQUIVALENCE
## 2                              2.2 Power=70%   19.970278 OSD_EQUIVALENCE
## 3                              2.4 Power=70%   11.233281 OSD_EQUIVALENCE
## 4                              2.6 Power=70%    7.189300 OSD_EQUIVALENCE
## 5                              2.8 Power=70%    4.992569 OSD_EQUIVALENCE
## 6                              3.0 Power=70%    3.668010 OSD_EQUIVALENCE
## 7                              3.2 Power=70%    2.808320 OSD_EQUIVALENCE
## 8                              3.4 Power=70%    2.218920 OSD_EQUIVALENCE
## 9                              3.6 Power=70%    1.797325 OSD_EQUIVALENCE
## 10                             3.8 Power=70%    1.485393 OSD_EQUIVALENCE
## 11                             4.0 Power=70%    1.248142 OSD_EQUIVALENCE
## 12                             2.0 Power=80%   53.524046 OSD_EQUIVALENCE
## 13                             2.2 Power=80%   23.788465 OSD_EQUIVALENCE
## 14                             2.4 Power=80%   13.381011 OSD_EQUIVALENCE
## 15                             2.6 Power=80%    8.563847 OSD_EQUIVALENCE
## 16                             2.8 Power=80%    5.947116 OSD_EQUIVALENCE
## 17                             3.0 Power=80%    4.369310 OSD_EQUIVALENCE
## 18                             3.2 Power=80%    3.345253 OSD_EQUIVALENCE
## 19                             3.4 Power=80%    2.643163 OSD_EQUIVALENCE
## 20                             3.6 Power=80%    2.140962 OSD_EQUIVALENCE
## 21                             3.8 Power=80%    1.769390 OSD_EQUIVALENCE
## 22                             4.0 Power=80%    1.486779 OSD_EQUIVALENCE
## 23                             2.0 Power=90%   67.638586 OSD_EQUIVALENCE
## 24                             2.2 Power=90%   30.061594 OSD_EQUIVALENCE
## 25                             2.4 Power=90%   16.909647 OSD_EQUIVALENCE
## 26                             2.6 Power=90%   10.822174 OSD_EQUIVALENCE
## 27                             2.8 Power=90%    7.515398 OSD_EQUIVALENCE
## 28                             3.0 Power=90%    5.521517 OSD_EQUIVALENCE
## 29                             3.2 Power=90%    4.227412 OSD_EQUIVALENCE
## 30                             3.4 Power=90%    3.340177 OSD_EQUIVALENCE
## 31                             3.6 Power=90%    2.705543 OSD_EQUIVALENCE
## 32                             3.8 Power=90%    2.235986 OSD_EQUIVALENCE
## 33                             4.0 Power=90%    1.878850 OSD_EQUIVALENCE
##################################
# Plotting the sample sizes and power curves
##################################
(OneSampleDesign.Equivalence.PowerAnalysis <- ggplot(OneSampleDesign.Equivalence.Reshaped,
                                                      aes(x=Posttreatment.Mean.Bone.Density,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Post-Treatment Mean Bone Density",
                      limits=c(2,4),
                      breaks=seq(2,4,
                                 by=0.2)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50)) +
  ggtitle("One-Sample Design : Test for Equivalence"))

1.1.5 Consolidated Evaluation


[A] The minimum sample size required for the treatment group in a One-Sample Design study is dependent upon the following factors defined and hypothesized prior to the clinical research:
     [A.1] Study Hypothesis (Equality, Non-Inferiority, Superiority, Equivalence)
     [A.2] Level of significance (Alpha) : Lower value increases sample size
     [A.3] Power of the test (Beta) : Higher value increases sample size
     [A.4] Data variability (Standard Deviation) : Higher value increases sample size
     [A.5] Epsilon (Difference between Mean Response and Reference) : Lower value increases sample size
     [A.6] Delta (Testing margin, as applicable) : Higher value increases sample size
[B] Given the level of significance = 5%, standard deviation = 1 and mean bone density values before treatment = 1.5 SD and after treatment = 2.0 to 4.0 SD representing epsilon values = 0.5 to 2.5, the range of sample sizes is given on the respective power analyses provided for each study hypotheses -
     [B.1] OSD_EQUALITY :
            [B.1.1] For 70% power, sample size = 1 to 25
            [B.1.2] For 80% power, sample size = 2 to 32
            [B.1.3] For 90% power, sample size = 2 to 43
     [B.2] OSD_NONINFERIORITY with Delta = -0.50 SD :
            [B.2.1] For 70% power, sample size = 1 to 5
            [B.2.2] For 80% power, sample size = 1 to 7
            [B.2.3] For 90% power, sample size = 1 to 9
     [B.3] OSD_SUPERIORITY with Delta = +0.10 SD :
            [B.3.1] For 70% power, sample size = 1 to 30
            [B.3.2] For 80% power, sample size = 2 to 39
            [B.3.3] For 90% power, sample size = 2 to 54
     [B.4] OSD_EQUIVALENCE with Delta = +0.10 SD:
            [B.4.1] For 70% power, sample size = 2 to 45
            [B.4.2] For 80% power, sample size = 2 to 54
            [B.4.3] For 90% power, sample size = 2 to 68

Code Chunk | Output
##################################
# Consolidating the power analyses charts
##################################

OSD_EQUALITY.PowerAnalysis <- ggplot(OneSampleDesign.Equality.Reshaped,
                                     aes(x=Posttreatment.Mean.Bone.Density,
                                         y=Sample.Size, 
                                         color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   facet_grid(. ~ Label) +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Post-Treatment Mean Bone Density",
                      limits=c(2,4),
                      breaks=seq(2,4,by=0.2)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50))

OSD_NONINFERIORITY.PowerAnalysis <- ggplot(OneSampleDesign.NonInferiority.Reshaped,
                                     aes(x=Posttreatment.Mean.Bone.Density,
                                         y=Sample.Size, 
                                         color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   facet_grid(. ~ Label) +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Post-Treatment Mean Bone Density",
                      limits=c(2,4),
                      breaks=seq(2,4,by=0.2)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50))

OSD_SUPERIORITY.PowerAnalysis <- ggplot(OneSampleDesign.Superiority.Reshaped,
                                     aes(x=Posttreatment.Mean.Bone.Density,
                                         y=Sample.Size, 
                                         color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   facet_grid(. ~ Label) +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Post-Treatment Mean Bone Density",
                      limits=c(2,4),
                      breaks=seq(2,4,by=0.2)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50))

OSD_EQUIVALENCE.PowerAnalysis <- ggplot(OneSampleDesign.Equivalence.Reshaped,
                                     aes(x=Posttreatment.Mean.Bone.Density,
                                         y=Sample.Size, 
                                         color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   facet_grid(. ~ Label) +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Post-Treatment Mean Bone Density",
                      limits=c(2,4),
                      breaks=seq(2,4,by=0.2)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50))

OSD_PowerAnalysis <- ggarrange(OSD_EQUALITY.PowerAnalysis,
                               OSD_NONINFERIORITY.PowerAnalysis,
                               OSD_SUPERIORITY.PowerAnalysis,
                               OSD_EQUIVALENCE.PowerAnalysis,
                               ncol=2, nrow=2)

annotate_figure(OSD_PowerAnalysis,
                top = text_grob("One-Sample Design Power Analysis by Statistical Hypothesis",
                                color = "black",
                                face = "bold",
                                size = 14))

1.2 Two-Sample Design (TSD)


Two-Sample Design involves individual numeric response values from sampled subjects between treatment and control groups, or between two competing groups. In clinical research, responses could be between a test drug and a placebo control or an active control agent; or between two competing agents. Responses are assumed to be independent normal random variables with zero mean and a given pooled variance. Parameters needed for the computations are as follows:

[A] Mean Treatment Response : True mean response value obtained from sampled subjects in the treatment group

[B] Mean Control Response : True mean response value obtained from sampled subjects in the control group

[C] Epsilon : Mean response differences between the treatment and control groups

[D] Standard Deviation : Data variability based from the from sampled subjects from both the Treatment and Control groups

[E] Alpha : Maximum allowed Type I error

[F] Beta : Maximum allowed Type II error

[G] Kappa : Matching ratio between the Treatment and Control groups

[H] Delta :
     [H.1] Minimum testing margin to establish that the mean treatment response is much better than control
     [H.2] Minimum testing margin to establish that the mean treatment response is much better than control
     [H.3] Minimum testing margin to establish that the mean treatment response is much better than control

A case scenario was provided as follows to demonstrate sample size computation and power analysis as applied on different types of hypotheses:

[Case Scenario] A clinical research study is being pursued concerning the evaluation of the effect of a test drug on cholesterol in patients with coronary heart disease (CHD). Cholesterol is the main lipid associated with arteriosclerotic vascular disease. The purpose of cholesterol testing is to identify patients at risk for arteriosclerotic heart disease. The liver metabolizes the cholesterol to its free form which is transported in the bloodtream by lipoproteins. As indicated by Pagana and Pagana (1998), nearly 75% of the cholesterol is bound to low density lipoproteins (LDLs) and 25% is bound to high density lipoproteins (HDLs). Therefore, cholesterol is the main component of LDLs and only a minimal component of HDLs and very low density lipoproteins. LDL is the most directly associated with increased risk of CHD. A pharmaceutical company is interested in conducting a clinical trial to compare two cholesterol lowering agents for treatment of patients with CHD through a parallel design. The primary efficacy parameter is the LDL.

1.2.1 Test for Equality (TSD_EQUALITY)


[Research Question] Suppose that a difference of 6 to 12% in percent change of LDL is hypothesized between two cholesterol lowering agents. The study aims to demonstrate that the mean LDL percent change between the two cholesterol lowering agents are different. Given a level of significance of 0.05, standard deviation of 10% and matching ratio between both groups equal to 1, evaluate the sample size range for a test of equality to achieve statistical test powers equal to 70%, 80% and 90%.

[Null Hypothesis] Lowering agent 1 mean LDL percent change - Lowering agent 2 mean LDL percent change equals 0

[Alternative Hypothesis] Lowering agent 1 mean LDL percent change - Lowering agent 2 mean LDL percent does not equal 0

[Hypothesis Test] Rejection of the null hypothesis implies difference between the mean LDL percent change between both lowering agents

[A] Epsilon : 6 to 12% at 1% intervals

[B] Standard Deviation : 10%

[C] Alpha : 0.05

[D] Beta : 0.30, 0.20 and 0.10 corresponding to 70%, 80% and 90% statistical powers

[E] Kappa : 1.0

[F] Delta : 0 testing margin for the Test of Equality hypothesis

Code Chunk | Output
##################################
# Defining a range of possible hypothesized values for epsilon
##################################
epsilon=seq(6,12,1)

##################################
# Defining a fixed hypothesized value for the standard deviation
##################################
sd=10

##################################
# Defining a fixed hypothesized value for the Type I error
##################################
alpha=0.05

##################################
# Defining a range of possible hypothesized values for the Type II error
##################################
beta=c(0.30,0.20,0.10)

##################################
# Defining a fixed hypothesized value for kappa
##################################
kappa=1.0

##################################
# Computing the range of samples sizes based on different levels of Type II errors
##################################
beta_1=beta[1]
(n_Beta30=(1+1/kappa)*(sd*(qnorm(1-alpha/2)+qnorm(1-beta_1))/(epsilon))^2)
## [1] 34.289261 25.192110 19.287709 15.239671 12.344134 10.201764  8.572315
ceiling(n_Beta30)
## [1] 35 26 20 16 13 11  9
z_Beta30=(epsilon)/(sd*sqrt((1+1/kappa)/n_Beta30))
(Power_Beta30=pnorm(z_Beta30-qnorm(1-alpha/2))+pnorm(-z_Beta30-qnorm(1-alpha/2)))
## [1] 0.7000044 0.7000044 0.7000044 0.7000044 0.7000044 0.7000044 0.7000044
beta_2=beta[2]
(n_Beta20=(1+1/kappa)*(sd*(qnorm(1-alpha/2)+qnorm(1-beta_2))/(epsilon))^2)
## [1] 43.60489 32.03624 24.52775 19.37995 15.69776 12.97335 10.90122
ceiling(n_Beta20)
## [1] 44 33 25 20 16 13 11
z_Beta20=(epsilon)/(sd*sqrt((1+1/kappa)/n_Beta20))
(Power_Beta20=pnorm(z_Beta20-qnorm(1-alpha/2))+pnorm(-z_Beta20-qnorm(1-alpha/2)))
## [1] 0.800001 0.800001 0.800001 0.800001 0.800001 0.800001 0.800001
beta_3=beta[3]
(n_Beta10=(1+1/kappa)*(sd*(qnorm(1-alpha/2)+qnorm(1-beta_3))/(epsilon))^2)
## [1] 58.37457 42.88744 32.83570 25.94425 21.01485 17.36764 14.59364
ceiling(n_Beta10)
## [1] 59 43 33 26 22 18 15
z_Beta10=(epsilon)/(sd*sqrt((1+1/kappa)/n_Beta10))
(Power_Beta10=pnorm(z_Beta10-qnorm(1-alpha/2))+pnorm(-z_Beta10-qnorm(1-alpha/2)))
## [1] 0.9000001 0.9000001 0.9000001 0.9000001 0.9000001 0.9000001 0.9000001
TwoSampleDesign.Equality <- as.data.frame(cbind(epsilon,n_Beta30,n_Beta20,n_Beta10))
names(TwoSampleDesign.Equality) <- c("Mean.LDL.Percent.Change",
                                 "Power=70%",
                                 "Power=80%",
                                 "Power=90%")

##################################
# Restructuring the data
##################################
TwoSampleDesign.Equality.Reshaped <- gather(TwoSampleDesign.Equality,
                                            "Power=70%","Power=80%","Power=90%",
                                            key="Power",
                                            value="Sample.Size")
TwoSampleDesign.Equality.Reshaped$Label <- rep("TSD_EQUALITY",nrow(TwoSampleDesign.Equality.Reshaped))

TwoSampleDesign.Equality.Reshaped
##    Mean.LDL.Percent.Change     Power Sample.Size        Label
## 1                        6 Power=70%   34.289261 TSD_EQUALITY
## 2                        7 Power=70%   25.192110 TSD_EQUALITY
## 3                        8 Power=70%   19.287709 TSD_EQUALITY
## 4                        9 Power=70%   15.239671 TSD_EQUALITY
## 5                       10 Power=70%   12.344134 TSD_EQUALITY
## 6                       11 Power=70%   10.201764 TSD_EQUALITY
## 7                       12 Power=70%    8.572315 TSD_EQUALITY
## 8                        6 Power=80%   43.604887 TSD_EQUALITY
## 9                        7 Power=80%   32.036244 TSD_EQUALITY
## 10                       8 Power=80%   24.527749 TSD_EQUALITY
## 11                       9 Power=80%   19.379950 TSD_EQUALITY
## 12                      10 Power=80%   15.697759 TSD_EQUALITY
## 13                      11 Power=80%   12.973355 TSD_EQUALITY
## 14                      12 Power=80%   10.901222 TSD_EQUALITY
## 15                       6 Power=90%   58.374573 TSD_EQUALITY
## 16                       7 Power=90%   42.887441 TSD_EQUALITY
## 17                       8 Power=90%   32.835697 TSD_EQUALITY
## 18                       9 Power=90%   25.944254 TSD_EQUALITY
## 19                      10 Power=90%   21.014846 TSD_EQUALITY
## 20                      11 Power=90%   17.367641 TSD_EQUALITY
## 21                      12 Power=90%   14.593643 TSD_EQUALITY
##################################
# Plotting the sample sizes and power curves
##################################
(TwoSampleDesign.Equality.PowerAnalysis <- ggplot(TwoSampleDesign.Equality.Reshaped,
                                                      aes(x=Mean.LDL.Percent.Change,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Mean LDL Percent Change",
                      limits=c(6,12),
                      breaks=seq(6,12,by=1)) +
  scale_y_continuous(name="Sample Size Per Group",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50)) +
  ggtitle("Two-Sample Design : Test for Equality"))

1.2.2 Test for Non-Inferiority (TSD_NONINFERIORITY)


[Research Question] Suppose that a difference of 6 to 12% in percent change of LDL is hypothesized between two cholesterol lowering agents. The study aims to demonstrate that the mean LDL percent change for the first cholesterol lowering agent is not much worse than the second cholesterol lowering agent by a clinically meaningful difference equal to −5. Given a level of significance of 0.05, standard deviation of 10% and matching ratio between both groups equal to 1, evaluate the sample size range for a test of non-inferiority to achieve statistical test powers equal to 70%, 80% and 90%.

[Null Hypothesis] Lowering agent 1 mean LDL percent change - Lowering agent 2 mean LDL percent change <= Clinically meaningful difference

[Alternative Hypothesis] Lowering agent 1 mean LDL percent change - Lowering agent 2 mean LDL percent > Clinically meaningful difference

[Hypothesis Test] Rejection of the null hypothesis implies non-inferiority of the Lowering agent 1 mean LDL percent change as compared to Lowering agent 1, given the testing margin

[A] Epsilon : 6 to 12% at 1% intervals

[B] Standard Deviation : 10%

[C] Alpha : 0.05

[D] Beta : 0.30, 0.20 and 0.10 corresponding to 70%, 80% and 90% statistical powers

[E] Kappa : 1.0

[F] Delta : -5 testing margin for the Test of Non-Inferiority hypothesis

Code Chunk | Output
##################################
# Defining a range of possible hypothesized values for epsilon
##################################
epsilon=seq(6,12,1)

##################################
# Defining a fixed hypothesized value for the standard deviation
##################################
sd=10

##################################
# Defining a fixed hypothesized value for the Type I error
##################################
alpha=0.05

##################################
# Defining a range of possible hypothesized values for the Type II error
##################################
beta=c(0.30,0.20,0.10)

##################################
# Defining a fixed hypothesized value for kappa
##################################
kappa=1.0

##################################
# Defining the non-inferiority delta
##################################
delta=-5

##################################
# Computing the range of samples sizes based on different levels of Type II errors
##################################
beta_1=beta[1]
(n_Beta30=(1+1/kappa)*(sd*(qnorm(1-alpha)+qnorm(1-beta_1))/(epsilon-delta))^2)
## [1] 7.777956 6.535644 5.568833 4.801697 4.182812 3.676300 3.256515
ceiling(n_Beta30)
## [1] 8 7 6 5 5 4 4
z_Beta30=(epsilon-delta)/(sd*sqrt((1+1/kappa)/n_Beta30))
(Power_Beta30=pnorm(z_Beta30-qnorm(1-alpha))+pnorm(-z_Beta30-qnorm(1-alpha)))
## [1] 0.7000683 0.7000683 0.7000683 0.7000683 0.7000683 0.7000683 0.7000683
beta_2=beta[2]
(n_Beta20=(1+1/kappa)*(sd*(qnorm(1-alpha)+qnorm(1-beta_2))/(epsilon-delta))^2)
## [1] 10.219103  8.586885  7.316636  6.308732  5.495606  4.830123  4.278586
ceiling(n_Beta20)
## [1] 11  9  8  7  6  5  5
z_Beta20=(epsilon-delta)/(sd*sqrt((1+1/kappa)/n_Beta20))
(Power_Beta20=pnorm(z_Beta20-qnorm(1-alpha))+pnorm(-z_Beta20-qnorm(1-alpha)))
## [1] 0.800018 0.800018 0.800018 0.800018 0.800018 0.800018 0.800018
beta_3=beta[3]
(n_Beta10=(1+1/kappa)*(sd*(qnorm(1-alpha)+qnorm(1-beta_3))/(epsilon-delta))^2)
## [1] 14.155120 11.894232 10.134731  8.738620  7.612309  6.690506  5.926538
ceiling(n_Beta10)
## [1] 15 12 11  9  8  7  6
z_Beta10=(epsilon-delta)/(sd*sqrt((1+1/kappa)/n_Beta10))
(Power_Beta10=pnorm(z_Beta10-qnorm(1-alpha))+pnorm(-z_Beta10-qnorm(1-alpha)))
## [1] 0.9000024 0.9000024 0.9000024 0.9000024 0.9000024 0.9000024 0.9000024
TwoSampleDesign.NonInferiority <- as.data.frame(cbind(epsilon,n_Beta30,n_Beta20,n_Beta10))
names(TwoSampleDesign.NonInferiority) <- c("Mean.LDL.Percent.Change",
                                 "Power=70%",
                                 "Power=80%",
                                 "Power=90%")

##################################
# Restructuring the data
##################################
TwoSampleDesign.NonInferiority.Reshaped <- gather(TwoSampleDesign.NonInferiority,
                                            "Power=70%","Power=80%","Power=90%",
                                            key="Power",
                                            value="Sample.Size")
TwoSampleDesign.NonInferiority.Reshaped$Label <- rep("TSD_NONINFERIORITY",nrow(TwoSampleDesign.NonInferiority.Reshaped))

TwoSampleDesign.NonInferiority.Reshaped
##    Mean.LDL.Percent.Change     Power Sample.Size              Label
## 1                        6 Power=70%    7.777956 TSD_NONINFERIORITY
## 2                        7 Power=70%    6.535644 TSD_NONINFERIORITY
## 3                        8 Power=70%    5.568833 TSD_NONINFERIORITY
## 4                        9 Power=70%    4.801697 TSD_NONINFERIORITY
## 5                       10 Power=70%    4.182812 TSD_NONINFERIORITY
## 6                       11 Power=70%    3.676300 TSD_NONINFERIORITY
## 7                       12 Power=70%    3.256515 TSD_NONINFERIORITY
## 8                        6 Power=80%   10.219103 TSD_NONINFERIORITY
## 9                        7 Power=80%    8.586885 TSD_NONINFERIORITY
## 10                       8 Power=80%    7.316636 TSD_NONINFERIORITY
## 11                       9 Power=80%    6.308732 TSD_NONINFERIORITY
## 12                      10 Power=80%    5.495606 TSD_NONINFERIORITY
## 13                      11 Power=80%    4.830123 TSD_NONINFERIORITY
## 14                      12 Power=80%    4.278586 TSD_NONINFERIORITY
## 15                       6 Power=90%   14.155120 TSD_NONINFERIORITY
## 16                       7 Power=90%   11.894232 TSD_NONINFERIORITY
## 17                       8 Power=90%   10.134731 TSD_NONINFERIORITY
## 18                       9 Power=90%    8.738620 TSD_NONINFERIORITY
## 19                      10 Power=90%    7.612309 TSD_NONINFERIORITY
## 20                      11 Power=90%    6.690506 TSD_NONINFERIORITY
## 21                      12 Power=90%    5.926538 TSD_NONINFERIORITY
##################################
# Plotting the sample sizes and power curves
##################################
(TwoSampleDesign.NonInferiority.PowerAnalysis <-ggplot(TwoSampleDesign.NonInferiority.Reshaped,
                                                      aes(x=Mean.LDL.Percent.Change,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Mean LDL Percent Change",
                      limits=c(6,12),
                      breaks=seq(6,12,by=1)) +
  scale_y_continuous(name="Sample Size Per Group",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50)) +
  ggtitle("Two-Sample Design : Test for Non-Inferiority"))

1.2.3 Test for Superiority (TSD_SUPERIORITY)


[Research Question] Suppose that a difference of 6 to 12% in percent change of LDL is hypothesized between two cholesterol lowering agents. The study aims to demonstrate that the mean LDL percent change for the first cholesterol lowering agent is much better than the second cholesterol lowering agent by a clinically meaningful difference equal to +5. Given a level of significance of 0.05, standard deviation of 10% and matching ratio between both groups equal to 1, evaluate the sample size range for a test of superiority to achieve statistical test powers equal to 70%, 80% and 90%.

[Null Hypothesis] Lowering agent 1 mean LDL percent change - Lowering agent 2 mean LDL percent change <= Clinically meaningful difference

[Alternative Hypothesis] Lowering agent 1 mean LDL percent change - Lowering agent 2 mean LDL percent > Clinically meaningful difference

[Hypothesis Test] Rejection of the null hypothesis implies superiority of the Lowering agent 1 mean LDL percent change as compared to Lowering agent 1, given the testing margin

[A] Epsilon : 6 to 12% at 1% intervals

[B] Standard Deviation : 10%

[C] Alpha : 0.05

[D] Beta : 0.30, 0.20 and 0.10 corresponding to 70%, 80% and 90% statistical powers

[E] Kappa : 1.0

[F] Delta : +1 testing margin for the Test of Superiority hypothesis

Code Chunk | Output
##################################
# Defining a range of possible hypothesized values for epsilon
##################################
epsilon=seq(6,12,1)

##################################
# Defining a fixed hypothesized value for the standard deviation
##################################
sd=10

##################################
# Defining a fixed hypothesized value for the Type I error
##################################
alpha=0.05

##################################
# Defining a range of possible hypothesized values for the Type II error
##################################
beta=c(0.30,0.20,0.10)

##################################
# Defining a fixed hypothesized value for kappa
##################################
kappa=1.0

##################################
# Defining the superiority delta
##################################
delta=1

##################################
# Computing the range of samples sizes based on different levels of Type II errors
##################################
beta_1=beta[1]
(n_Beta30=(1+1/kappa)*(sd*(qnorm(1-alpha)+qnorm(1-beta_1))/(epsilon-delta))^2)
## [1] 37.645308 26.142575 19.206790 14.705199 11.618922  9.411327  7.777956
ceiling(n_Beta30)
## [1] 38 27 20 15 12 10  8
z_Beta30=(epsilon-delta)/(sd*sqrt((1+1/kappa)/n_Beta30))
(Power_Beta30=pnorm(z_Beta30-qnorm(1-alpha))+pnorm(-z_Beta30-qnorm(1-alpha)))
## [1] 0.7000683 0.7000683 0.7000683 0.7000683 0.7000683 0.7000683 0.7000683
beta_2=beta[2]
(n_Beta20=(1+1/kappa)*(sd*(qnorm(1-alpha)+qnorm(1-beta_2))/(epsilon-delta))^2)
## [1] 49.46046 34.34754 25.23493 19.32049 15.26557 12.36511 10.21910
ceiling(n_Beta20)
## [1] 50 35 26 20 16 13 11
z_Beta20=(epsilon-delta)/(sd*sqrt((1+1/kappa)/n_Beta20))
(Power_Beta20=pnorm(z_Beta20-qnorm(1-alpha))+pnorm(-z_Beta20-qnorm(1-alpha)))
## [1] 0.800018 0.800018 0.800018 0.800018 0.800018 0.800018 0.800018
beta_3=beta[3]
(n_Beta10=(1+1/kappa)*(sd*(qnorm(1-alpha)+qnorm(1-beta_3))/(epsilon-delta))^2)
## [1] 68.51078 47.57693 34.95448 26.76202 21.14530 17.12769 14.15512
ceiling(n_Beta10)
## [1] 69 48 35 27 22 18 15
z_Beta10=(epsilon-delta)/(sd*sqrt((1+1/kappa)/n_Beta10))
(Power_Beta10=pnorm(z_Beta10-qnorm(1-alpha))+pnorm(-z_Beta10-qnorm(1-alpha)))
## [1] 0.9000024 0.9000024 0.9000024 0.9000024 0.9000024 0.9000024 0.9000024
TwoSampleDesign.Superiority <- as.data.frame(cbind(epsilon,n_Beta30,n_Beta20,n_Beta10))
names(TwoSampleDesign.Superiority) <- c("Mean.LDL.Percent.Change",
                                 "Power=70%",
                                 "Power=80%",
                                 "Power=90%")

##################################
# Restructuring the data
##################################
TwoSampleDesign.Superiority.Reshaped <- gather(TwoSampleDesign.Superiority,
                                            "Power=70%","Power=80%","Power=90%",
                                            key="Power",
                                            value="Sample.Size")
TwoSampleDesign.Superiority.Reshaped$Label <- rep("TSD_SUPERIORITY",nrow(TwoSampleDesign.Superiority.Reshaped))

TwoSampleDesign.Superiority.Reshaped
##    Mean.LDL.Percent.Change     Power Sample.Size           Label
## 1                        6 Power=70%   37.645308 TSD_SUPERIORITY
## 2                        7 Power=70%   26.142575 TSD_SUPERIORITY
## 3                        8 Power=70%   19.206790 TSD_SUPERIORITY
## 4                        9 Power=70%   14.705199 TSD_SUPERIORITY
## 5                       10 Power=70%   11.618922 TSD_SUPERIORITY
## 6                       11 Power=70%    9.411327 TSD_SUPERIORITY
## 7                       12 Power=70%    7.777956 TSD_SUPERIORITY
## 8                        6 Power=80%   49.460458 TSD_SUPERIORITY
## 9                        7 Power=80%   34.347540 TSD_SUPERIORITY
## 10                       8 Power=80%   25.234927 TSD_SUPERIORITY
## 11                       9 Power=80%   19.320491 TSD_SUPERIORITY
## 12                      10 Power=80%   15.265573 TSD_SUPERIORITY
## 13                      11 Power=80%   12.365114 TSD_SUPERIORITY
## 14                      12 Power=80%   10.219103 TSD_SUPERIORITY
## 15                       6 Power=90%   68.510779 TSD_SUPERIORITY
## 16                       7 Power=90%   47.576930 TSD_SUPERIORITY
## 17                       8 Power=90%   34.954479 TSD_SUPERIORITY
## 18                       9 Power=90%   26.762023 TSD_SUPERIORITY
## 19                      10 Power=90%   21.145302 TSD_SUPERIORITY
## 20                      11 Power=90%   17.127695 TSD_SUPERIORITY
## 21                      12 Power=90%   14.155120 TSD_SUPERIORITY
##################################
# Plotting the sample sizes and power curves
##################################
(TwoSampleDesign.Superiority.PowerAnalysis <-ggplot(TwoSampleDesign.Superiority.Reshaped,
                                                      aes(x=Mean.LDL.Percent.Change,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Mean LDL Percent Change",
                      limits=c(6,12),
                      breaks=seq(6,12,by=1)) +
  scale_y_continuous(name="Sample Size Per Group",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50)) +
  ggtitle("Two-Sample Design : Test for Superiority"))

1.2.4 Test for Equivalence (TSD_EQUIVALENCE)


[Research Question] Suppose that a difference of 6 to 12% in percent change of LDL is hypothesized between two cholesterol lowering agents. The study aims to demonstrate that the mean LDL percent change for the first cholesterol lowering agent is not much worse or better than the second cholesterol lowering agent by a clinically meaningful difference equal to +5. Given a level of significance of 0.05, standard deviation of 10% and matching ratio between both groups equal to 1, evaluate the sample size range for a test of equivalence to achieve statistical test powers equal to 70%, 80% and 90%.

[Null Hypothesis] Absolute value of (Lowering agent 1 mean LDL percent change - Lowering agent 2 mean LDL percent change) >= Clinically meaningful difference

[Alternative Hypothesis] Absolute value of (Lowering agent 1 mean LDL percent change - Lowering agent 2 mean LDL percent) < Clinically meaningful difference

[Hypothesis Test] Rejection of the null hypothesis implies equivalence of the Lowering agent 1 mean LDL percent change as compared to Lowering agent 1, given the testing margin

[A] Epsilon : 6 to 12% at 1% intervals

[B] Standard Deviation : 10%

[C] Alpha : 0.05

[D] Beta : 0.30, 0.20 and 0.10 corresponding to 70%, 80% and 90% statistical powers

[E] Kappa : 1.0

[F] Delta : +1 testing margin for the Test of Equivalence hypothesis

Code Chunk | Output
##################################
# Defining a range of possible hypothesized values for epsilon
##################################
epsilon=seq(6,12,1)

##################################
# Defining a fixed hypothesized value for the standard deviation
##################################
sd=10

##################################
# Defining a fixed hypothesized value for the Type I error
##################################
alpha=0.05

##################################
# Defining a range of possible hypothesized values for the Type II error
##################################
beta=c(0.30,0.20,0.10)

##################################
# Defining a fixed hypothesized value for kappa
##################################
kappa=1.0

##################################
# Defining the equivalence delta
##################################
delta=1

##################################
# Computing the range of samples sizes based on different levels of Type II errors
##################################
beta_1=beta[1]
(n_Beta30=(1+1/kappa)*(sd*(qnorm(1-alpha)+qnorm(1-beta_1/2))/(abs(epsilon)-delta))^2)
## [1] 57.51440 39.94056 29.34408 22.46656 17.75136 14.37860 11.88314
ceiling(n_Beta30)
## [1] 58 40 30 23 18 15 12
z_Beta30=(abs(epsilon)-delta)/(sd*sqrt((1+1/kappa)/n_Beta30))
(Power_Beta30=2*(pnorm(z_Beta30-qnorm(1-alpha))+pnorm(-z_Beta30-qnorm(1-alpha)))-1)
## [1] 0.7000152 0.7000152 0.7000152 0.7000152 0.7000152 0.7000152 0.7000152
beta_2=beta[2]
(n_Beta20=(1+1/kappa)*(sd*(qnorm(1-alpha)+qnorm(1-beta_2/2))/(abs(epsilon)-delta))^2)
## [1] 68.51078 47.57693 34.95448 26.76202 21.14530 17.12769 14.15512
ceiling(n_Beta20)
## [1] 69 48 35 27 22 18 15
z_Beta20=(abs(epsilon)-delta)/(sd*sqrt((1+1/kappa)/n_Beta20))
(Power_Beta20=2*(pnorm(z_Beta20-qnorm(1-alpha))+pnorm(-z_Beta20-qnorm(1-alpha)))-1)
## [1] 0.8000048 0.8000048 0.8000048 0.8000048 0.8000048 0.8000048 0.8000048
beta_3=beta[3]
(n_Beta10=(1+1/kappa)*(sd*(qnorm(1-alpha)+qnorm(1-beta_3/2))/(abs(epsilon)-delta))^2)
## [1] 86.57739 60.12319 44.17214 33.81929 26.72142 21.64435 17.88789
ceiling(n_Beta10)
## [1] 87 61 45 34 27 22 18
z_Beta10=(abs(epsilon)-delta)/(sd*sqrt((1+1/kappa)/n_Beta10))
(Power_Beta10=2*(pnorm(z_Beta10-qnorm(1-alpha))+pnorm(-z_Beta10-qnorm(1-alpha)))-1)
## [1] 0.9000008 0.9000008 0.9000008 0.9000008 0.9000008 0.9000008 0.9000008
TwoSampleDesign.Equivalence <- as.data.frame(cbind(epsilon,n_Beta30,n_Beta20,n_Beta10))
names(TwoSampleDesign.Equivalence) <- c("Mean.LDL.Percent.Change",
                                 "Power=70%",
                                 "Power=80%",
                                 "Power=90%")

##################################
# Restructuring the data
##################################
TwoSampleDesign.Equivalence.Reshaped <- gather(TwoSampleDesign.Equivalence,
                                            "Power=70%","Power=80%","Power=90%",
                                            key="Power",
                                            value="Sample.Size")
TwoSampleDesign.Equivalence.Reshaped$Label <- rep("TSD_EQUIVALENCE",nrow(TwoSampleDesign.Equivalence.Reshaped))

TwoSampleDesign.Equivalence.Reshaped
##    Mean.LDL.Percent.Change     Power Sample.Size           Label
## 1                        6 Power=70%    57.51440 TSD_EQUIVALENCE
## 2                        7 Power=70%    39.94056 TSD_EQUIVALENCE
## 3                        8 Power=70%    29.34408 TSD_EQUIVALENCE
## 4                        9 Power=70%    22.46656 TSD_EQUIVALENCE
## 5                       10 Power=70%    17.75136 TSD_EQUIVALENCE
## 6                       11 Power=70%    14.37860 TSD_EQUIVALENCE
## 7                       12 Power=70%    11.88314 TSD_EQUIVALENCE
## 8                        6 Power=80%    68.51078 TSD_EQUIVALENCE
## 9                        7 Power=80%    47.57693 TSD_EQUIVALENCE
## 10                       8 Power=80%    34.95448 TSD_EQUIVALENCE
## 11                       9 Power=80%    26.76202 TSD_EQUIVALENCE
## 12                      10 Power=80%    21.14530 TSD_EQUIVALENCE
## 13                      11 Power=80%    17.12769 TSD_EQUIVALENCE
## 14                      12 Power=80%    14.15512 TSD_EQUIVALENCE
## 15                       6 Power=90%    86.57739 TSD_EQUIVALENCE
## 16                       7 Power=90%    60.12319 TSD_EQUIVALENCE
## 17                       8 Power=90%    44.17214 TSD_EQUIVALENCE
## 18                       9 Power=90%    33.81929 TSD_EQUIVALENCE
## 19                      10 Power=90%    26.72142 TSD_EQUIVALENCE
## 20                      11 Power=90%    21.64435 TSD_EQUIVALENCE
## 21                      12 Power=90%    17.88789 TSD_EQUIVALENCE
##################################
# Plotting the sample sizes and power curves
##################################
(TwoSampleDesign.Equivalence.PowerAnalysis <-ggplot(TwoSampleDesign.Equivalence.Reshaped,
                                                      aes(x=Mean.LDL.Percent.Change,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Mean LDL Percent Change",
                      limits=c(6,12),
                      breaks=seq(6,12,by=1)) +
  scale_y_continuous(name="Sample Size Per Group",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50)) +
  ggtitle("Two-Sample Design : Test for Equivalence"))

1.2.5 Consolidated Evaluation


[A] The minimum sample size required for the treatment group in a Two-Sample Design study is dependent upon the following factors defined and hypothesized prior to the clinical research:
     [A.1] Study Hypothesis (Equality, Non-Inferiority, Superiority, Equivalence)
     [A.2] Level of significance (Alpha) : Lower value increases sample size
     [A.3] Power of the test (Beta) : Higher value increases sample size
     [A.4] Data variability (Standard Deviation) : Higher value increases sample size
     [A.5] Epsilon (Difference between Mean Response and Reference) : Lower value increases sample size
     [A.6] Delta (Testing margin, as applicable) : Higher value increases sample size
[B] Given the level of significance = 5%, standard deviation = 10, matching sample size ratio between groups = 1 and LDL percent change between two cholesterol lowering agents representing epsilon values = 6 to 12%, the range of sample sizes is given on the respective power analyses provided for each study hypotheses -
     [B.1] OSD_EQUALITY :
            [B.1.1] For 70% power, sample size = 9 to 35
            [B.1.2] For 80% power, sample size = 11 to 44
            [B.1.3] For 90% power, sample size = 15 to 59
     [B.2] OSD_NONINFERIORITY with Delta = -0.50 SD :
            [B.2.1] For 70% power, sample size = 4 to 8
            [B.2.2] For 80% power, sample size = 5 to 11
            [B.2.3] For 90% power, sample size = 6 to 15
     [B.3] OSD_SUPERIORITY with Delta = +0.10 SD :
            [B.3.1] For 70% power, sample size = 8 to 38
            [B.3.2] For 80% power, sample size = 11 to 50
            [B.3.3] For 90% power, sample size = 15 to 69
     [B.4] OSD_EQUIVALENCE with Delta = +0.10 SD:
            [B.4.1] For 70% power, sample size = 12 to 58
            [B.4.2] For 80% power, sample size = 15 to 69
            [B.4.3] For 90% power, sample size = 18 to 87

Code Chunk | Output
##################################
# Consolidating the power analyses charts
##################################
TSD_EQUALITY.PowerAnalysis <- ggplot(TwoSampleDesign.Equality.Reshaped,
                                     aes(x=Mean.LDL.Percent.Change,
                                         y=Sample.Size, 
                                         color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   facet_grid(. ~ Label) +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Mean LDL Percent Change",
                      limits=c(6,12),
                      breaks=seq(6,12,by=1)) +
  scale_y_continuous(name="Sample Size Per Group",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50))

TSD_NONINFERIORITY.PowerAnalysis <- ggplot(TwoSampleDesign.NonInferiority.Reshaped,
                                     aes(x=Mean.LDL.Percent.Change,
                                         y=Sample.Size, 
                                         color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   facet_grid(. ~ Label) +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Mean LDL Percent Change",
                      limits=c(6,12),
                      breaks=seq(6,12,by=1)) +
  scale_y_continuous(name="Sample Size Per Group",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50))

TSD_SUPERIORITY.PowerAnalysis <- ggplot(TwoSampleDesign.Superiority.Reshaped,
                                     aes(x=Mean.LDL.Percent.Change,
                                         y=Sample.Size, 
                                         color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   facet_grid(. ~ Label) +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Mean LDL Percent Change",
                      limits=c(6,12),
                      breaks=seq(6,12,by=1)) +
  scale_y_continuous(name="Sample Size Per Group",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50))

TSD_EQUIVALENCE.PowerAnalysis <- ggplot(TwoSampleDesign.Equivalence.Reshaped,
                                     aes(x=Mean.LDL.Percent.Change,
                                         y=Sample.Size, 
                                         color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   facet_grid(. ~ Label) +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Mean LDL Percent Change",
                      limits=c(6,12),
                      breaks=seq(6,12,by=1)) +
  scale_y_continuous(name="Sample Size Per Group",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50))

TSD_PowerAnalysis <- ggarrange(TSD_EQUALITY.PowerAnalysis,
                               TSD_NONINFERIORITY.PowerAnalysis,
                               TSD_SUPERIORITY.PowerAnalysis,
                               TSD_EQUIVALENCE.PowerAnalysis,
                               ncol=2, nrow=2)

annotate_figure(TSD_PowerAnalysis,
                top = text_grob("Two-Sample Design Power Analysis by Statistical Hypothesis",
                                color = "black",
                                face = "bold",
                                size = 14))

1.3 Multiple-Sample Design One-Way ANOVA Pairwise (MSDOWAP)


Multiple-Sample Design One-Way ANOVA Pairwise involves individual numeric response values from sampled subjects across three or more competing treatment groups. In clinical research, responses could be between multiple test drugs and a placebo control or an active control agent; or between multiple competing agents. The response and treatment variables are modeled using a One-Way ANOVA model involving the fixed effect of the treatment and random errors which are assumed to be independent and identically distributed normal random variables with zero mean and a given variance. For a pairwise comparison between treatments, all possible combinations must be analyzed. Due to the multiple comparisons inflating the Type I error, adjustments arre made for controlling the overall Type I error rate at the desired significance level. Parameters needed for the computations are as follows:

[A] Mean Treatment 1 Response : True mean response value obtained from sampled subjects in the first treatment group

[B] Mean Treatment 2 Response : True mean response value obtained from sampled subjects in the second treatment group

[C] Epsilon : Mean response differences between the two treatment groups

[D] Standard Deviation : Data variability based from the from sampled subjects from both treatment groups

[E] Alpha : Maximum allowed Type I error

[F] Beta : Maximum allowed Type II error

[G] Tau : Number of treatment groups

A case scenario was provided as follows to demonstrate sample size computation and power analysis as applied on different treatment pairs:

[Case Scenario] A clinical research study was conducted to investigate pharmaceutical compounds for treatment of patients with cancer. The investigator is not only interested in showing efficacy of the test drug but also in establishing the dose response curve. To achieve this study objective, a three-arm parallel group, double-blind, randomized clinical trial was designed to compare three treatments: A (Placebo), B (10 mg) and C (30 mg). Based on the results from the pilot study, it is estimated that the standard deviation is 3.5.

1.3.1 Test for Equality (MSDOWAP_EQUALITY)


[Research Question] Suppose that the hypothesized mean clinical response for the placebo treatment is 9.0. However, clinical responses for Treatments B and C were expected to range from 11.0 to 13.0 and 15.0 to 17.0, respectively, with increments of 0.5. The study aims to demonstrate the efficacy of each treatment as pairwise compared against each other. Given a level of significance of 0.05, evaluate the sample size range for a test of equality to achieve statistical test powers equal to 70%, 80% and 90%.

[Null Hypothesis] Mean clinical response of the first evaluated treatment equals the Mean clinical response of the second evaluated treatment

[Alternative Hypothesis] Mean clinical response of the first evaluated treatment is not equal to the Mean clinical response of the second evaluated treatment

[Hypothesis Test] Rejection of the null hypothesis implies difference between the mean clinical response values between both evaluated treatments

[A] Mean Response : Treatment A : 9.0

[B] Mean Response : Treatment B : 11.0 to 13.0 at 0.5 intervals

[C] Mean Response : Treatment C : 15.0 to 17.0 at 0.5 intervals

[D] Standard Deviation : 3.50

[E] Alpha : 0.05

[F] Beta : 0.30, 0.20 and 0.10 corresponding to 70%, 80% and 90% statistical powers

[G] Tau : 3 representing the number of treatments

Code Chunk | Output
##################################
# Defining a fixed hypothesized value for mu_A
##################################
mu_AF=9.0

##################################
# Defining a range of hypothesized values for mu_B
##################################
mu_BR=seq(11.0,13.0,0.5)
##################################
# Defining a fixed hypothesized value for mu_B
# using the highest value from the range
##################################
mu_BF=13.0

##################################
# Defining a range of hypothesized value for mu_C
##################################
mu_CR=seq(15.0,17.0,0.5)

##################################
# Defining a fixed hypothesized value for the standard deviation
##################################
sd=3.5

##################################
# Defining a fixed hypothesized value for the Type I error
##################################
alpha=0.05

##################################
# Defining a range of possible hypothesized values for the Type II error
##################################
beta=c(0.30,0.20,0.10)

##################################
# Defining the number of treatments
##################################
tau=3

##################################
# Conducting a pairwise comparison
# between Treatments A and B
##################################

##################################
# Computing the range of samples sizes based on different levels of Type II errors
##################################
beta_1=beta[1]
(n_Beta30=2*(sd*(qnorm(1-alpha/(2/tau))+qnorm(1-beta_1))/(mu_BR-mu_AF))^2)
## [1] 23.624302 15.119553 10.499690  7.714058  5.906075
ceiling(n_Beta30)
## [1] 24 16 11  8  6
z_Beta30=(mu_BR-mu_AF)/(sd*sqrt(2/n_Beta30))
(Power_Beta30=pnorm(z_Beta30-qnorm(1-alpha/(2/tau)))+pnorm(-z_Beta30-qnorm(1-alpha/(2/tau))))
## [1] 0.7003327 0.7003327 0.7003327 0.7003327 0.7003327
beta_2=beta[2]
(n_Beta20=2*(sd*(qnorm(1-alpha/(2/tau))+qnorm(1-beta_2))/(mu_BR-mu_AF))^2)
## [1] 31.872403 20.398338 14.165513 10.407315  7.968101
ceiling(n_Beta20)
## [1] 32 21 15 11  8
z_Beta20=(mu_BR-mu_AF)/(sd*sqrt(2/n_Beta20))
(Power_Beta20=pnorm(z_Beta20-qnorm(1-alpha/(2/tau)))+pnorm(-z_Beta20-qnorm(1-alpha/(2/tau))))
## [1] 0.8000993 0.8000993 0.8000993 0.8000993 0.8000993
beta_3=beta[3]
(n_Beta10=2*(sd*(qnorm(1-alpha/(2/tau))+qnorm(1-beta_3))/(mu_BR-mu_AF))^2)
## [1] 45.35129 29.02483 20.15613 14.80859 11.33782
ceiling(n_Beta10)
## [1] 46 30 21 15 12
z_Beta10=(mu_BR-mu_AF)/(sd*sqrt(2/n_Beta10))
(Power_Beta10=pnorm(z_Beta10-qnorm(1-alpha/(2/tau)))+pnorm(-z_Beta10-qnorm(1-alpha/(2/tau))))
## [1] 0.9000159 0.9000159 0.9000159 0.9000159 0.9000159
MultipleSampleDesign.AB.Equality <- as.data.frame(cbind(mu_BR,n_Beta30,n_Beta20,n_Beta10))
names(MultipleSampleDesign.AB.Equality) <- c("Clinical.Response",
                                 "Power=70%",
                                 "Power=80%",
                                 "Power=90%")

##################################
# Restructuring the data
##################################
MultipleSampleDesign.AB.Equality.Reshaped <- gather(MultipleSampleDesign.AB.Equality,
                                            "Power=70%","Power=80%","Power=90%",
                                            key="Power",
                                            value="Sample.Size")
MultipleSampleDesign.AB.Equality.Reshaped$Label <- rep("MSDOWAP_AB_EQUALITY",nrow(MultipleSampleDesign.AB.Equality.Reshaped))

MultipleSampleDesign.AB.Equality.Reshaped
##    Clinical.Response     Power Sample.Size               Label
## 1               11.0 Power=70%   23.624302 MSDOWAP_AB_EQUALITY
## 2               11.5 Power=70%   15.119553 MSDOWAP_AB_EQUALITY
## 3               12.0 Power=70%   10.499690 MSDOWAP_AB_EQUALITY
## 4               12.5 Power=70%    7.714058 MSDOWAP_AB_EQUALITY
## 5               13.0 Power=70%    5.906075 MSDOWAP_AB_EQUALITY
## 6               11.0 Power=80%   31.872403 MSDOWAP_AB_EQUALITY
## 7               11.5 Power=80%   20.398338 MSDOWAP_AB_EQUALITY
## 8               12.0 Power=80%   14.165513 MSDOWAP_AB_EQUALITY
## 9               12.5 Power=80%   10.407315 MSDOWAP_AB_EQUALITY
## 10              13.0 Power=80%    7.968101 MSDOWAP_AB_EQUALITY
## 11              11.0 Power=90%   45.351294 MSDOWAP_AB_EQUALITY
## 12              11.5 Power=90%   29.024828 MSDOWAP_AB_EQUALITY
## 13              12.0 Power=90%   20.156131 MSDOWAP_AB_EQUALITY
## 14              12.5 Power=90%   14.808586 MSDOWAP_AB_EQUALITY
## 15              13.0 Power=90%   11.337823 MSDOWAP_AB_EQUALITY
##################################
# Plotting the sample sizes and power curves
##################################
(MultipleSampleDesign.AB.Equality.PowerAnalysis <-    ggplot(MultipleSampleDesign.AB.Equality.Reshaped,
                                                      aes(x=Clinical.Response,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Clinical Response",
                      limits=c(11,13),
                      breaks=seq(11,13,by=0.5)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50)) +
  ggtitle("Multiple-Sample Design : Test for Equality Between A Versus B"))

##################################
# Conducting a pairwise comparison
# between Treatments A and C
##################################

##################################
# Computing the range of samples sizes based on different levels of Type II errors
##################################
beta_1=beta[1]
(n_Beta30=2*(sd*(qnorm(1-alpha/(2/tau))+qnorm(1-beta_1))/(mu_CR-mu_AF))^2)
## [1] 2.624922 2.236620 1.928514 1.679950 1.476519
ceiling(n_Beta30)
## [1] 3 3 2 2 2
z_Beta30=(mu_CR-mu_AF)/(sd*sqrt(2/n_Beta30))
(Power_Beta30=pnorm(z_Beta30-qnorm(1-alpha/(2/tau)))+pnorm(-z_Beta30-qnorm(1-alpha/(2/tau))))
## [1] 0.7003327 0.7003327 0.7003327 0.7003327 0.7003327
beta_2=beta[2]
(n_Beta20=2*(sd*(qnorm(1-alpha/(2/tau))+qnorm(1-beta_2))/(mu_CR-mu_AF))^2)
## [1] 3.541378 3.017506 2.601829 2.266482 1.992025
ceiling(n_Beta20)
## [1] 4 4 3 3 2
z_Beta20=(mu_CR-mu_AF)/(sd*sqrt(2/n_Beta20))
(Power_Beta20=pnorm(z_Beta20-qnorm(1-alpha/(2/tau)))+pnorm(-z_Beta20-qnorm(1-alpha/(2/tau))))
## [1] 0.8000993 0.8000993 0.8000993 0.8000993 0.8000993
beta_3=beta[3]
(n_Beta10=2*(sd*(qnorm(1-alpha/(2/tau))+qnorm(1-beta_3))/(mu_CR-mu_AF))^2)
## [1] 5.039033 4.293614 3.702146 3.224981 2.834456
ceiling(n_Beta10)
## [1] 6 5 4 4 3
z_Beta10=(mu_CR-mu_AF)/(sd*sqrt(2/n_Beta10))
(Power_Beta10=pnorm(z_Beta10-qnorm(1-alpha/(2/tau)))+pnorm(-z_Beta10-qnorm(1-alpha/(2/tau))))
## [1] 0.9000159 0.9000159 0.9000159 0.9000159 0.9000159
MultipleSampleDesign.AC.Equality <- as.data.frame(cbind(mu_CR,n_Beta30,n_Beta20,n_Beta10))
names(MultipleSampleDesign.AC.Equality) <- c("Clinical.Response",
                                 "Power=70%",
                                 "Power=80%",
                                 "Power=90%")

##################################
# Restructuring the data
##################################
MultipleSampleDesign.AC.Equality.Reshaped <- gather(MultipleSampleDesign.AC.Equality,
                                            "Power=70%","Power=80%","Power=90%",
                                            key="Power",
                                            value="Sample.Size")
MultipleSampleDesign.AC.Equality.Reshaped$Label <- rep("MSDOWAP_AC_EQUALITY",nrow(MultipleSampleDesign.AC.Equality.Reshaped))

MultipleSampleDesign.AC.Equality.Reshaped
##    Clinical.Response     Power Sample.Size               Label
## 1               15.0 Power=70%    2.624922 MSDOWAP_AC_EQUALITY
## 2               15.5 Power=70%    2.236620 MSDOWAP_AC_EQUALITY
## 3               16.0 Power=70%    1.928514 MSDOWAP_AC_EQUALITY
## 4               16.5 Power=70%    1.679950 MSDOWAP_AC_EQUALITY
## 5               17.0 Power=70%    1.476519 MSDOWAP_AC_EQUALITY
## 6               15.0 Power=80%    3.541378 MSDOWAP_AC_EQUALITY
## 7               15.5 Power=80%    3.017506 MSDOWAP_AC_EQUALITY
## 8               16.0 Power=80%    2.601829 MSDOWAP_AC_EQUALITY
## 9               16.5 Power=80%    2.266482 MSDOWAP_AC_EQUALITY
## 10              17.0 Power=80%    1.992025 MSDOWAP_AC_EQUALITY
## 11              15.0 Power=90%    5.039033 MSDOWAP_AC_EQUALITY
## 12              15.5 Power=90%    4.293614 MSDOWAP_AC_EQUALITY
## 13              16.0 Power=90%    3.702146 MSDOWAP_AC_EQUALITY
## 14              16.5 Power=90%    3.224981 MSDOWAP_AC_EQUALITY
## 15              17.0 Power=90%    2.834456 MSDOWAP_AC_EQUALITY
##################################
# Plotting the sample sizes and power curves
##################################
(MultipleSampleDesign.AC.Equality.PowerAnalysis <-    ggplot(MultipleSampleDesign.AC.Equality.Reshaped,
                                                      aes(x=Clinical.Response,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Clinical Response",
                      limits=c(15,17),
                      breaks=seq(15,17,by=0.5)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50)) +
  ggtitle("Multiple-Sample Design : Test for Equality Between A Versus C"))

##################################
# Conducting a pairwise comparison
# between Treatments B and C
##################################

##################################
# Computing the range of samples sizes based on different levels of Type II errors
##################################
beta_1=beta[1]
(n_Beta30=2*(sd*(qnorm(1-alpha/(2/tau))+qnorm(1-beta_1))/(mu_CR-mu_BF))^2)
## [1] 23.624302 15.119553 10.499690  7.714058  5.906075
ceiling(n_Beta30)
## [1] 24 16 11  8  6
z_Beta30=(mu_CR-mu_BF)/(sd*sqrt(2/n_Beta30))
(Power_Beta30=pnorm(z_Beta30-qnorm(1-alpha/(2/tau)))+pnorm(-z_Beta30-qnorm(1-alpha/(2/tau))))
## [1] 0.7003327 0.7003327 0.7003327 0.7003327 0.7003327
beta_2=beta[2]
(n_Beta20=2*(sd*(qnorm(1-alpha/(2/tau))+qnorm(1-beta_2))/(mu_CR-mu_BF))^2)
## [1] 31.872403 20.398338 14.165513 10.407315  7.968101
ceiling(n_Beta20)
## [1] 32 21 15 11  8
z_Beta20=(mu_CR-mu_BF)/(sd*sqrt(2/n_Beta20))
(Power_Beta20=pnorm(z_Beta20-qnorm(1-alpha/(2/tau)))+pnorm(-z_Beta20-qnorm(1-alpha/(2/tau))))
## [1] 0.8000993 0.8000993 0.8000993 0.8000993 0.8000993
beta_3=beta[3]
(n_Beta10=2*(sd*(qnorm(1-alpha/(2/tau))+qnorm(1-beta_3))/(mu_CR-mu_BF))^2)
## [1] 45.35129 29.02483 20.15613 14.80859 11.33782
ceiling(n_Beta10)
## [1] 46 30 21 15 12
z_Beta10=(mu_CR-mu_BF)/(sd*sqrt(2/n_Beta10))
(Power_Beta10=pnorm(z_Beta10-qnorm(1-alpha/(2/tau)))+pnorm(-z_Beta10-qnorm(1-alpha/(2/tau))))
## [1] 0.9000159 0.9000159 0.9000159 0.9000159 0.9000159
MultipleSampleDesign.BC.Equality <- as.data.frame(cbind(mu_CR,n_Beta30,n_Beta20,n_Beta10))
names(MultipleSampleDesign.BC.Equality) <- c("Clinical.Response",
                                 "Power=70%",
                                 "Power=80%",
                                 "Power=90%")

##################################
# Restructuring the data
##################################
MultipleSampleDesign.BC.Equality.Reshaped <- gather(MultipleSampleDesign.BC.Equality,
                                            "Power=70%","Power=80%","Power=90%",
                                            key="Power",
                                            value="Sample.Size")
MultipleSampleDesign.BC.Equality.Reshaped$Label <- rep("MSDOWAP_BC_EQUALITY",nrow(MultipleSampleDesign.BC.Equality.Reshaped))

MultipleSampleDesign.BC.Equality.Reshaped
##    Clinical.Response     Power Sample.Size               Label
## 1               15.0 Power=70%   23.624302 MSDOWAP_BC_EQUALITY
## 2               15.5 Power=70%   15.119553 MSDOWAP_BC_EQUALITY
## 3               16.0 Power=70%   10.499690 MSDOWAP_BC_EQUALITY
## 4               16.5 Power=70%    7.714058 MSDOWAP_BC_EQUALITY
## 5               17.0 Power=70%    5.906075 MSDOWAP_BC_EQUALITY
## 6               15.0 Power=80%   31.872403 MSDOWAP_BC_EQUALITY
## 7               15.5 Power=80%   20.398338 MSDOWAP_BC_EQUALITY
## 8               16.0 Power=80%   14.165513 MSDOWAP_BC_EQUALITY
## 9               16.5 Power=80%   10.407315 MSDOWAP_BC_EQUALITY
## 10              17.0 Power=80%    7.968101 MSDOWAP_BC_EQUALITY
## 11              15.0 Power=90%   45.351294 MSDOWAP_BC_EQUALITY
## 12              15.5 Power=90%   29.024828 MSDOWAP_BC_EQUALITY
## 13              16.0 Power=90%   20.156131 MSDOWAP_BC_EQUALITY
## 14              16.5 Power=90%   14.808586 MSDOWAP_BC_EQUALITY
## 15              17.0 Power=90%   11.337823 MSDOWAP_BC_EQUALITY
##################################
# Plotting the sample sizes and power curves
##################################
(MultipleSampleDesign.BC.Equality.PowerAnalysis <-    ggplot(MultipleSampleDesign.BC.Equality.Reshaped,
                                                      aes(x=Clinical.Response,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Clinical Response",
                      limits=c(15,17),
                      breaks=seq(15,17,by=0.5)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50)) +
  ggtitle("Multiple-Sample Design : Test for Equality Between B Versus C"))

1.3.2 Consolidated Evaluation


[A] The minimum sample size required for each treatment group in a Multiple-Sample Design One-Way ANOVA Pairwise study is dependent upon the following factors defined and hypothesized prior to the clinical research:
     [A.1] Level of significance (Alpha) : Lower value increases sample size
     [A.2] Power of the test (Beta) : Higher value increases sample size
     [A.3] Data variability (Standard Deviation) : Higher value increases sample size
     [A.4] Epsilon (Difference between pairwise Mean Responses) : Lower value increases sample size
     [A.5] Tau (Number of treatment groups) : Higher value increases sample size
[B] Given the level of significance = 5%, standard deviation = 3.5, number of treatments = 3 and a varying range of clinical response values across pairwise combinations of treatments, the sample size is determined as the maximum number obtained for all computations under each statistical power -
     [B.1] 70% power, sample size = 24
     [B.2] 80% power, sample size = 32
     [B.3] 90% power, sample size = 45

Code Chunk | Output
##################################
# Consolidating the power analyses charts
##################################
MSDOWAP.AB.Equality.PowerAnalysis <- ggplot(MultipleSampleDesign.AB.Equality.Reshaped,
                                                      aes(x=Clinical.Response,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   facet_grid(. ~ Label) +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Clinical Response",
                      limits=c(11,13),
                      breaks=seq(11,13,by=0.5)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50))

MSDOWAP.AC.Equality.PowerAnalysis <- ggplot(MultipleSampleDesign.AC.Equality.Reshaped,
                                                      aes(x=Clinical.Response,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   facet_grid(. ~ Label) +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Clinical Response",
                      limits=c(15,17),
                      breaks=seq(15,17,by=0.5)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50))

MSDOWAP.BC.Equality.PowerAnalysis <- ggplot(MultipleSampleDesign.BC.Equality.Reshaped,
                                                      aes(x=Clinical.Response,
                                                          y=Sample.Size, 
                                                          color=Power)) +
   geom_line(size=1)+
   geom_point(size=4)+
   theme_bw() +
   facet_grid(. ~ Label) +
   scale_color_brewer(palette="Paired") +
   scale_x_continuous(name="Hypothesized Clinical Response",
                      limits=c(15,17),
                      breaks=seq(15,17,by=0.5)) +
  scale_y_continuous(name="Treatment Group Sample Size",
                     limits=c(0,100),
                     breaks=seq(0,100,by=10)) +
  theme(axis.title.x=element_text(color="black",face="bold",size=12),
        legend.position="top",
        axis.title.y=element_text(color="black",face="bold",size=12),
        plot.title=element_text(color="black",size=14,face="bold",hjust=0.50))


MSDOWAP_PowerAnalysis <- ggarrange(MSDOWAP.AB.Equality.PowerAnalysis,
                                   MSDOWAP.AC.Equality.PowerAnalysis,
                                   MSDOWAP.BC.Equality.PowerAnalysis,
                                   ncol=2, nrow=2)

annotate_figure(MSDOWAP_PowerAnalysis,
                top = text_grob("Multiple-Sample Design One-Way ANOVA Power Analysis by Pairwise Comparison",
                                color = "black",
                                face = "bold",
                                size = 14))

2. Summary



3. References


[Book] Sample Size Calculations in Clinical Research by Shein-Chung Chow, Jun Shao, Hansheng Wang and Yuliya Lokhnygina
[Book] Clinical Trial Data Analysis Using R and SAS by Ding-Geng Chen, Karl Peace and Pinggao Zhang
[Book] Design and Analysis of Experiments with R by John Lawson
[Book] Design and Analysis of Experiments by Angela Dean , Daniel Voss and Danel Draguljić
[Book] Design and Analysis of Experiments by Douglas Montgomery
[R Package] rstatix by Alboukadel Kassambara
[R Package] car by John Fox and Sanford Weisberg
[R Package] effects by John Fox and Sanford Weisberg
[R Package] multcomp by Torsten Hothorn, Frank Bretz and Peter Westfall
[R Package] psych by William Revelle
[R Package] bootstrap by Robert Tibshirani
[R Package] moments by Lukasz Komsta and Frederick Novomestky
[R Package] dplyr by Hadley Wickham
[R Package] ggplot2 by Hadley Wickham
[R Package] ggpubr by Alboukadel Kassambara
[R Package] ggfortify by Yuan Tang
[R Package] trend by Thorsten Pohlert
[Article] Power and Sample Size by HyLown Consulting Team
[Article] An Introduction to Different Types of Study Design by Hadi Abbas
[Article] Case-Control and Cohort Studies: A Brief Overview by Saul Crandon
[Article] Cohort Studies: Prospective and Retrospective Designs by Izabel de Oliveira
[Article] Prevalence vs. Incidence: What is the Difference by Georgina Ford
[Article] Randomized Clinical Trial (RCT): Simple Definition, Phases, and Types by Statistics-How-To Team