Usage
Sample size calculation with extraSuperpower is done in three steps:
Define the study groups and expected outcomes from your experiment, possibly including interaction effects. If your design has repeated measurements, generate a covariance matrix as well.
Simulate the outcomes with a user specified sample size or a set of sample sizes. Outcomes may have a normal, skewed normal or truncated normal distribution.
Calculate the power for main effects and interaction based on the simulated data.
Creating a cell mean model for independent or repeated measures
The minimal required input is the number of levels \(a\) of factor \(A\), number of levels \(b\) of factor \(B\), the mean outcome value of cell \(a_1, b_1\) (generally a baseline value in control participants) and the ratio of change from one level to the next in each of the factors.
A list with the names of the factors and each of the levels of each factor is highly recommended.
If interaction effects are expected, the combination of levels and the ratio of change of this effect is required.
If the design includes repeated measurements, the correlation within subjects and the factor for which measures are repeated must be specified.
We will start with an example of an independent measures design in which interaction is not expected. Factor A is treatment group simply labelled as “Group” and factor B is timepoint or simply “Time”.
## outcome mean in reference group at baseline is 10
## a control group and an intervention group will be compared over 3 timepoints
## all measurements are independent
refmean <- 10
Alevs <- 2
Blevs <- 3
fAeff <- 1.5
fBeff <- 0.8
## if you do not provide a list with names of factors and levels, factor names are to "fA" and "fB" and level names are set to 'letters[1:nlfA]' and 'letters[1:nlfB]'.
Alevelnames <- c("control", "intervention")
Blevelnames <- 1:Blevs
nameslist <- list("Group" = Alevelnames, "Time" = Blevelnames)
simple_twoway <- calculate_mean_matrix(refmean = refmean, nlfA = Alevs, nlfB = Blevs,
fAeffect = fAeff, fBeffect = fBeff,
label_list = nameslist)In this case the output is two matrices and a plot. The matrices are the cell means and standard deviations for each combination of levels of factors \(A\) and \(B\).
##labelling factors and their levels is convenient
simple_twoway
#> $matrices_obj
#> $matrices_obj$mean.mat
#> Time
#> Group 1 2 3
#> control 10 9 8
#> intervention 15 14 13
#>
#> $matrices_obj$sd.mat
#> Time
#> Group 1 2 3
#> control 2 1.8 1.6
#> intervention 3 2.8 2.6
#>
#>
#> $meansplot
As a default, the standard deviation is one fifth of the mean value of each cell. This proportion can be changed with the ‘sdproportional’ and the ‘sdratio’ options. By setting ‘sdproportional = FALSE’ standard deviation will be the same in all cells and estimated as a proportion of the mean of all cells in the mean matrix. By setting ‘sdratio = 0.1’ this proportion is 10%.
For all other parameters we will use the values from the previous example.
simple_twoway_sdadjusted <- calculate_mean_matrix(refmean = refmean, nlfA = Alevs, nlfB = Blevs,
fAeffect = fAeff, fBeffect = fBeff,
sdproportional = FALSE, sdratio = 0.1,
label_list = nameslist)
simple_twoway_sdadjusted
#> $matrices_obj
#> $matrices_obj$mean.mat
#> Time
#> Group 1 2 3
#> control 10 9 8
#> intervention 15 14 13
#>
#> $matrices_obj$sd.mat
#> [1] 1.15
#>
#>
#> $meansplot
Including an interaction effect
Now, instead of not expecting an interaction effect we are expecting it in the intervention group at times 2 and 3. The other parameters will be as in the first example.
#intervention group is the second row in the means matrix, times 2 and 3 the 2nd and 3rd columns.
cellsinteraction <- c(2, 2, 2, 3)
cellsinteraction <- matrix(cellsinteraction, 2, 2)
interaction_twoway <- calculate_mean_matrix(refmean = refmean, nlfA = Alevs, nlfB = Blevs,
fAeffect = fAeff, fBeffect = fBeff,
groupswinteraction = cellsinteraction, interact = 0.7,
label_list = nameslist)
interaction_twoway
#> $matrices_obj
#> $matrices_obj$mean.mat
#> Time
#> Group 1 2 3
#> control 10 9.0 8.0
#> intervention 15 9.8 9.1
#>
#> $matrices_obj$sd.mat
#> Time
#> Group 1 2 3
#> control 2 1.80 1.60
#> intervention 3 1.96 1.82
#>
#>
#> $meansplot
Modeling a repeated measures experiment
To end this section, we will switch from an independent measurements design to a design were treatment is a between factor, while “Time” is a within factor. In other words, the same study participant has been measured over the different levels of factor “Time”. All other parameters will stay as in the independent measures interaction example.
#Let's suppose within subject correlation is 0.7
rho <- 0.7
interaction_twoway_timewithin <- calculate_mean_matrix(refmean = refmean, nlfA = Alevs, nlfB = Blevs,
fAeffect = fAeff, fBeffect = fBeff,
groupswinteraction = cellsinteraction, interact = 0.7,
rho = rho, withinf = "fB",
label_list = nameslist)
interaction_twoway_timewithin
#> $matrices_obj
#> $matrices_obj$within.factor
#> [1] "fB"
#>
#> $matrices_obj$mean.mat
#> Time
#> Group 1 2 3
#> control 10 9.0 8.0
#> intervention 15 9.8 9.1
#>
#> $matrices_obj$sd.mat
#> Time
#> Group 1 2 3
#> control 2 1.80 1.60
#> intervention 3 1.96 1.82
#>
#> $matrices_obj$cormat
#> control_1 control_2 control_3 intervention_1 intervention_2
#> control_1 1.0 0.7 0.7 0.0 0.0
#> control_2 0.7 1.0 0.7 0.0 0.0
#> control_3 0.7 0.7 1.0 0.0 0.0
#> intervention_1 0.0 0.0 0.0 1.0 0.7
#> intervention_2 0.0 0.0 0.0 0.7 1.0
#> intervention_3 0.0 0.0 0.0 0.7 0.7
#> intervention_3
#> control_1 0.0
#> control_2 0.0
#> control_3 0.0
#> intervention_1 0.7
#> intervention_2 0.7
#> intervention_3 1.0
#>
#> $matrices_obj$sigmat
#> control_1 control_2 control_3 intervention_1 intervention_2
#> control_1 4.00 2.520 2.240 0.000 0.00000
#> control_2 2.52 3.240 2.016 0.000 0.00000
#> control_3 2.24 2.016 2.560 0.000 0.00000
#> intervention_1 0.00 0.000 0.000 9.000 4.11600
#> intervention_2 0.00 0.000 0.000 4.116 3.84160
#> intervention_3 0.00 0.000 0.000 3.822 2.49704
#> intervention_3
#> control_1 0.00000
#> control_2 0.00000
#> control_3 0.00000
#> intervention_1 3.82200
#> intervention_2 2.49704
#> intervention_3 3.31240
#>
#>
#> $meansplot
The cell matrices and the plot are the same, the difference lies in that we obtain correlation and covariance matrices as additional output. The ‘withinf’ option can take the values of ‘NULL’ (in which case the ‘calculate_mean_matrix’ will assume independent measurements), ‘fA’, ‘fB’ or ‘both’, depending on which factor is the “within” factor in the design.
For the within factor correlation we can provide a constant or, if more than 1 value is provided, a correlation gradient can be generated in the matrix.
Note Functions to generate correlation and covariance matrices may also be run separately. To generate te correlation matrix the required input is the mean matrix,
Simulating a two-way factorial experiment
Once we have set values for each combination of levels of factors \(A\) and \(B\) we will sample outcome values under these assumptions with different sample sizes. Sampling can be done from normal, skewed normal or truncated normal distributions. In the repeated measurements case these are multivariate distributions.
Lets start with the independent measurement experiment in which we expect there will be interaction that we stored in the ‘interaction_twoway’ object. We will sample from the normal distribution, which is the default.
We will use a low number of iterations for these examples. The
coverage of the power calculations will be wide. The default value for
nsims is 150. With this value power estimation may take a
couple of minutes in the context of the examples that follow.
iterations <- 50
set.seed(170824)
n <- seq(6, 12, 3)
indepmeasures_normal_sim <- simulate_twoway_nrange(matrices_obj = interaction_twoway,
nset = n, distribution = "normal", nsims = iterations)
#> Simulating independent observations experiment
length(indepmeasures_normal_sim)
#> [1] 3
length(n)
#> [1] 3The output is a list in which each element is a simulation of sample size specified by the ‘nset’ vector.
To sample from a skewed normal distribution.
indepmeasures_skewed_sim <- simulate_twoway_nrange(matrices_obj = interaction_twoway,
nset = n, distribution = "skewed", skewness = 2,
nsims = iterations)
#> Simulating independent observations experimentThe same function is used to sample from a multivariate distribution but you must set the ‘repeated_measurements’ option to ‘TRUE’.
repmeasures_normal_sim <- simulate_twoway_nrange(matrices_obj = interaction_twoway_timewithin,
nset = n, repeated_measurements = TRUE,
nsims = iterations)
#> Simulating repeated observations experimentTo end this section, we sample from a multivariate skewed normal distribution. The skewness parameter must be length 1, \(a\), \(b\) or \(ab\).
Estimating power from simulated data
The final step is estimating the power under the different simulated sample sizes. We will start with the independent measurement example. Default test is ANOVA.
test_power_overkn(indepmeasures_normal_sim)
#> Testing power on an independent observations design experiment.
#> Sample size per group = 6
#> Testing power on an independent observations design experiment.
#> Sample size per group = 9
#> Testing power on an independent observations design experiment.
#> Sample size per group = 12
#> $power_table
#> n power lower.bound.ci upper.bound.ci effect
#> Group 6 0.86 0.7638 0.9562 Group
#> Time 6 1.00 1.0000 1.0000 Time
#> Group:Time 6 0.54 0.4019 0.6781 Group:Time
#> Group1 9 0.98 0.9412 1.0188 Group
#> Time1 9 1.00 1.0000 1.0000 Time
#> Group:Time1 9 0.86 0.7638 0.9562 Group:Time
#> Group2 12 0.98 0.9412 1.0188 Group
#> Time2 12 1.00 1.0000 1.0000 Time
#> Group:Time2 12 0.94 0.8742 1.0058 Group:Time
#>
#> $power_curve
We go on with the simulation in which sampling was from a skewed distribution.
test_power_overkn(indepmeasures_skewed_sim)
#> Testing power on an independent observations design experiment.
#> Sample size per group = 6
#> Testing power on an independent observations design experiment.
#> Sample size per group = 9
#> Testing power on an independent observations design experiment.
#> Sample size per group = 12
#> $power_table
#> n power lower.bound.ci upper.bound.ci effect
#> Group 6 0.92 0.8448 0.9952 Group
#> Time 6 1.00 1.0000 1.0000 Time
#> Group:Time 6 0.68 0.5507 0.8093 Group:Time
#> Group1 9 0.96 0.9057 1.0143 Group
#> Time1 9 1.00 1.0000 1.0000 Time
#> Group:Time1 9 0.88 0.7899 0.9701 Group:Time
#> Group2 12 1.00 1.0000 1.0000 Group
#> Time2 12 1.00 1.0000 1.0000 Time
#> Group:Time2 12 0.88 0.7899 0.9701 Group:Time
#>
#> $power_curve
Now we use a estimate the power of the rank test on the simulation sampled from the skewed normal distribution.
test_power_overkn(indepmeasures_skewed_sim, test = "rank")
#> Testing power on an independent observations design experiment.
#> Sample size per group = 6
#> Testing power on an independent observations design experiment.
#> Sample size per group = 9
#> Testing power on an independent observations design experiment.
#> Sample size per group = 12
#> $power_table
#> n power lower.bound.ci upper.bound.ci effect
#> Group 6 0.94 0.8742 1.0058 Group
#> Time 6 1.00 1.0000 1.0000 Time
#> Group:Time 6 0.78 0.6652 0.8948 Group:Time
#> Group1 9 0.94 0.8742 1.0058 Group
#> Time1 9 1.00 1.0000 1.0000 Time
#> Group:Time1 9 0.92 0.8448 0.9952 Group:Time
#> Group2 12 1.00 1.0000 1.0000 Group
#> Time2 12 1.00 1.0000 1.0000 Time
#> Group:Time2 12 0.86 0.7638 0.9562 Group:Time
#>
#> $power_curve
Likewise, we estimate the power for the simulations sampled from multivariate normal and multivariate skewed normal distributions.
test_power_overkn(repmeasures_normal_sim)
#> Testing power on a repeated observations design experiment.
#> Sample size = 6
#> Testing power on a repeated observations design experiment.
#> Sample size = 9
#> Testing power on a repeated observations design experiment.
#> Sample size = 12
#> $power_table
#> n power lower.bound.ci upper.bound.ci effect
#> Group 6 0.50 0.3614 0.6386 Group
#> Time 6 1.00 1.0000 1.0000 Time
#> Group:Time 6 0.94 0.8742 1.0058 Group:Time
#> Group1 9 0.68 0.5507 0.8093 Group
#> Time1 9 1.00 1.0000 1.0000 Time
#> Group:Time1 9 1.00 1.0000 1.0000 Group:Time
#> Group2 12 0.82 0.7135 0.9265 Group
#> Time2 12 1.00 1.0000 1.0000 Time
#> Group:Time2 12 1.00 1.0000 1.0000 Group:Time
#>
#> $power_curve
In this case the power for the group effect is apparently lower than the independent measures design, while the power for interaction is higher. Lets see what happens with the rank test.
test_power_overkn(repmeasures_normal_sim, test = "rank")
#> $power_table
#> n power lower.bound.ci upper.bound.ci effect
#> Group 6 0.46 0.3219 0.5981 Group
#> Time 6 1.00 1.0000 1.0000 Time
#> Group:Time 6 0.36 0.2270 0.4930 Group:Time
#> Group1 9 0.54 0.4019 0.6781 Group
#> Time1 9 1.00 1.0000 1.0000 Time
#> Group:Time1 9 0.78 0.6652 0.8948 Group:Time
#> Group2 12 0.78 0.6652 0.8948 Group
#> Time2 12 1.00 1.0000 1.0000 Time
#> Group:Time2 12 0.88 0.7899 0.9701 Group:Time
#>
#> $power_curve
The power is lower using a rank test.
What about the simulation in which we sampled from the skewed distribution?
test_power_overkn(repmeasures_skewed_sim)
#> Testing power on a repeated observations design experiment.
#> Sample size = 6
#> Testing power on a repeated observations design experiment.
#> Sample size = 9
#> Testing power on a repeated observations design experiment.
#> Sample size = 12
#> $power_table
#> n power lower.bound.ci upper.bound.ci effect
#> Group 6 0.42 0.2832 0.5568 Group
#> Time 6 1.00 1.0000 1.0000 Time
#> Group:Time 6 0.98 0.9412 1.0188 Group:Time
#> Group1 9 0.70 0.5730 0.8270 Group
#> Time1 9 1.00 1.0000 1.0000 Time
#> Group:Time1 9 1.00 1.0000 1.0000 Group:Time
#> Group2 12 0.90 0.8168 0.9832 Group
#> Time2 12 1.00 1.0000 1.0000 Time
#> Group:Time2 12 1.00 1.0000 1.0000 Group:Time
#>
#> $power_curve
Finally, the rank test for a skewed outcome, which is probably more appropriate.
test_power_overkn(repmeasures_skewed_sim, test = "rank")
#> Testing power on a repeated observations design experiment.
#> Sample size = 6
#> Testing power on a repeated observations design experiment.
#> Sample size = 9
#> Testing power on a repeated observations design experiment.
#> Sample size = 12
#> $power_table
#> n power lower.bound.ci upper.bound.ci effect
#> Group 6 0.30 0.1730 0.4270 Group
#> Time 6 1.00 1.0000 1.0000 Time
#> Group:Time 6 0.50 0.3614 0.6386 Group:Time
#> Group1 9 0.50 0.3614 0.6386 Group
#> Time1 9 1.00 1.0000 1.0000 Time
#> Group:Time1 9 0.62 0.4855 0.7545 Group:Time
#> Group2 12 0.74 0.6184 0.8616 Group
#> Time2 12 1.00 1.0000 1.0000 Time
#> Group:Time2 12 0.82 0.7135 0.9265 Group:Time
#>
#> $power_curve