Create input for simulation based two-way factorial experiments

Description

This function will generate a matrix of expected mean values for ab factor level combinations of a two-way factorial design by assuming linear effects with possible departure from linearity by interaction. It will also provide a standard deviation matrix for these ab combinations of factor levels. If the design has repeated measures, it will additionally provide correlation and covariance matrices calculated depending on which factor has repeated measurements or is the 'within' factor.

Usage

calculate_mean_matrix(
  refmean,
  nlfA,
  nlfB,
  fAeffect,
  fBeffect,
  groupswinteraction = NULL,
  interact = 1,
  label_list = NULL,
  sdproportional = TRUE,
  sdratio = 0.2,
  endincrement = TRUE,
  rho = 0,
  withinf = NULL,
  plot = TRUE
)

Arguments

Details

The user must provide a reference mean (usually mean in control or untreated group), the expected change for each factor from first to last level (or from one level to the next) and the number of levels in each factor.

The user can also specify factor level combinations in which interaction is assumed and its magnitude with respect to the reference mean. The cell mean matrix will be modified accordingly and this can also have an effect of the standard deviation matrix.

We were motivated by sample size calculation for two-way factorial designs with 1,2,...,a levels of factor A and 1,2,...,b levels of factor B in which the mean outcome value for replicates of cell A=1, B=1 are known. Furthermore, there is an expected change in level mean for each of the factors. Finally, interaction can be explicitly introduced to level combinations in which it is expected to occur.

If a repeated measures experiment is intended withinf must be set to "fA", "fB" or "both", depending on which is the 'within' factor. If rho is a vector length 1, the within subject correlation will be constant for the factor defined in withinf. If rho is a vector length 2 and withinf is either "fA" or "fB" a correlation gradient will be created from the first to second value of rho. If rho is a vector length 2 and withinf="both", the first element of rho will be the correlation within factor A, while the second element will be the correlation within factor B. If rho is a 2*2 matrix, only possible if withinf="both", a correlation gradient will be created across rows of rho for each of the factors.

Value

If rho and withinf are left at their default values of 0 and NULL, respectively, a cell mean matrix, a cell standard deviation matrix and optionally a graph that represents both.

If rho is between -1 and 1 but different to 0 and withinf is either "fA", "fB" or "both", correlation and covariance matrices are generated along with the aforementioned output.

Examples

refmean <- 1
treatgroups <- 4
timepoints <- 5
treateff <- 1.5
timeeff <- 0.85
factors_levels_names <- list(treatment=letters[1:treatgroups], time=1:timepoints)
## Independent design
effects_treat_time <- calculate_mean_matrix(refmean = refmean,
                                            fAeffect = treateff, fBeffect = timeeff,
                                            nlfA = treatgroups, nlfB = timepoints,
                                            label_list = factors_levels_names)
## Inspect plot to check if matrices correspond to design
effects_treat_time$meansplot
n <- 20
independent_experiment <- twoway_simulation_independent(group_size = n,
                                      matrices_obj = effects_treat_time)

head(independent_experiment, 10)

## Repeated measures design, suppose subjects from 4 independent treatment groups measured
## at 5 different timepoints.
## We use the same parameters as the independent design example, except we add within factor level
## correlation and we specify that factor B is the within factor.


refmean <- 1
treatgroups <- 4
timepoints <- 5
treateff <- 1.5
timeeff <- 0.85
rho <- 0.8
withinf <- "fB"

factors_levels_names <- list(treatment=letters[1:treatgroups], time=1:timepoints)

effects_treat_time <- calculate_mean_matrix(refmean = refmean, fAeffect = treateff,
                      fBeffect = timeeff, nlfA = treatgroups,  nlfB = timepoints,
                      rho = rho, withinf = withinf, label_list = factors_levels_names)

## Plot should look the same, structure within data can be checked once simulated
effects_treat_time$meansplot

n <- 20
repeatedmeasures_experiment <- twoway_simulation_correlated(group_size = n,
                                          matrices_obj = effects_treat_time)
head(repeatedmeasures_experiment, 10)

Effect size calculation

Description

Calculate effect sizes for two-way factorial designs from matrices of expected mean and standard deviation values at each combination of factor levels. The output given is Cohen's f. Calculations are done as exemplified in the G*Power 3.1 manual.

Usage

effsize(matrices_obj)

Arguments

Value

Vector of length 3. The first two elements are the effect sizes for the main effects factor A and factor B, respectively. The third element is the interaction effect size.

Examples

# no interaction effect expected
refmean <- 1
treatgroups <- 4
timepoints <- 5
treateff <- 1.5
timeeff <- 0.85

factors_levels_names <- list(treatment=letters[1:treatgroups], time=1:timepoints)

effects_treat_time <- calculate_mean_matrix(refmean = refmean,
                                            nlfA = treatgroups, nlfB = timepoints,
                                            fAeffect = treateff, fBeffect = timeeff,
                                            label_list = factors_levels_names)
effsize(effects_treat_time)

# we add cell specific interaction effect keeping design and main effect coefficients
cellswithinteraction <- matrix(c(rep(2,3), 3:5), 3,2)
#second level of factor A interacts with 3rd, 4th and 5th level of factor B

effects_treat_time_interact <- calculate_mean_matrix(refmean = refmean,
                                                     nlfA = treatgroups, nlfB = timepoints,
                                                     fAeffect = treateff, fBeffect = timeeff,
                                                     label_list = factors_levels_names,
                                                     groupswinteraction = cellswithinteraction,
                                                     interact=1.3)
effsize(effects_treat_time_interact)

Two-way factorial ANOVA exact sample size calculation for independent samples

Description

This functions takes the effect sizes (Cohen's f) for two main effects and their interaction and estimates power a range of sample sizes. The input for this function can be generated by effsize.

Usage

exact_twoway_anova_power(
  a,
  b,
  effect_sizes,
  n,
  alpha = 0.05,
  factor_names = NULL,
  plot = TRUE,
  title = NULL,
  target_power = NULL,
  target_line = TRUE,
  alpha_line = TRUE
)

Arguments

Details

Probably the best way to calculate power for independent balanced designs

Value

A list that contains the number of levels for each factor, the chosen significance level and a data.frame in which the first column is the group sample size and the remaining three columns are the power for the main effect of the first factor, the main effect of the second factor and their interaction, respectively.

Optionally, a graph that displays the power curves.

Examples

refmean <- 1
treatgroups <- 4
timepoints <- 5
treateff <- 1.25
timeeff <- 0.85
cellswithinteraction <- matrix(c(rep(2,3), 3:5), 3,2)
#second level of factor A interacts with 3rd, 4th and 5th level of factor B

factors_levels_names <- list(treatment=letters[1:treatgroups], time=1:timepoints)

effects_treat_time_interact <- calculate_mean_matrix(refmean = refmean,
                                                     nlfA = treatgroups, nlfB = timepoints,
                                                     fAeffect = treateff, fBeffect = timeeff,
                                                     label_list = factors_levels_names,
                                                     groupswinteraction = cellswithinteraction,
                                                     interact=1.3)

fxs <- effsize(effects_treat_time_interact)
exact_twoway_anova_power(a= treatgroups, b=timepoints, effect_sizes=fxs, n=5:20)

Function that generates a correlation matrix taking as input number of factors for each level, factor or factors that present correlation and rho value or values. Additionally, a mean matrix is required to check consistency.

Description

May be run independently or internally as part of calculate_mean_matrix.

Usage

gencorrelationmat(mean_matrix, rho, label_list = NULL, withinf, nlfA, nlfB)

Arguments

Details

For a repeated measures experiment withinf must be set to "fA", "fB" or "both", depending on which is the 'within' factor. If rho is a vector length 1, the within subject correlation will be constant for the factor defined in withinf. If rho is a vector length 2 and withinf is either "fA" or "fB" a correlation gradient will be created from the first to second value of rho. If rho is a vector length 2 and withinf="both", the first element of rho will be the correlation within factor A, while the second element will be the correlation within factor B. If rho is a 2*2 matrix, only possible if withinf="both", a correlation gradient will be created across rows of rho for each of the factors.

Value

Correlation matrix

Examples

meanvals <- c(seq(3,9,2),seq(2,8,2),seq(1,7,2))
mean_matrix <- matrix(meanvals, 3, 4, byrow = TRUE,
                   dimnames = list(A=LETTERS[1:3], B=letters[1:4]))

mean_matrix

gencorrelationmat(mean_matrix = mean_matrix, rho = 0.7, withinf = "fB", nlfA = 3, nlfB = 4)

##correlation gradient over levels of factor B
gencorrelationmat(mean_matrix = mean_matrix, rho = c(0.7, 0.4), withinf = "fB", nlfA = 3, nlfB = 4)

##gradient both factors

rhovals <- matrix(c(0.7, 0.4), 2, 2, byrow = TRUE)
gencorrelationmat(mean_matrix = mean_matrix, rho = rhovals, withinf = "both", nlfA = 3, nlfB = 4)

Function that generates a covariance matrix taking as input a correlation matrix and a standard deviation matrix or value.

Description

May be run independently or internally as part of 'calculate_mean_matrix'.

Usage

gencovariancemat(
  correlation_matrix,
  sd_matrix,
  withinf,
  label_list = NULL,
  nlfA,
  nlfB
)

Arguments

Value

Covariance matrix

Examples

meanvals <- c(seq(3,9,2),seq(2,8,2),seq(1,7,2))
mean_matrix <- matrix(meanvals, 3, 4, byrow = TRUE,
                   dimnames = list(A=LETTERS[1:3], B=letters[1:4]))

mean_matrix
sd_matrix <- mean_matrix*0.2

cor_matrix <- gencorrelationmat(mean_matrix = mean_matrix,
                                rho = 0.7, withinf = "fB", nlfA = 3, nlfB = 4)

gencovariancemat(cor_matrix, sd_matrix, withinf = "fB", nlfA = 3, nlfB = 4)

##correlation gradient over levels of factor B
cor_matrix <- gencorrelationmat(mean_matrix = mean_matrix,
                                rho = c(0.7, 0.4), withinf = "fB", nlfA = 3, nlfB = 4)

gencovariancemat(cor_matrix, sd_matrix, withinf = "fB", nlfA = 3, nlfB = 4)

Graph modeled means and standard deviations of groups in two-way factorial design

Description

Internal function that plots modeled cell means and standard deviations and covariance matrices. Takes input generated by the calculate_mean_matrix function and runs inside of it.

Usage

graph_twoway_assumptions(group_size = 100, matrices_obj)

Arguments

Value

Line plot with expected mean and standard deviation for each combination of factor levels

Plots the output of test_twoway_nrange

Description

Internal function, called by test_twoway_nrange, to plot power against sample size.

Usage

plot_powercurves(
  power_over_nrange,
  target_power = NULL,
  title = NULL,
  target_line = TRUE,
  alpha_line = TRUE,
  alpha = 0.05
)

Arguments

Value

Plot with power curves.

Examples

## 'cornorm_model' is created with the calculate_mean_matrix function
refmean <- 10
treateff <- 1.2
timeeff <- 0.75

treatgroups <- 3
treatgroups_names <- c("wt", "DrugA", "DrugB")

timepoints <- 4
timepoints_names <- paste0("t", 1:timepoints)

nameslist <- list(treatment=treatgroups_names, time=timepoints_names)

rho = 0.7

cornorm_model <- calculate_mean_matrix(refmean = refmean, fAeffect = treateff, fBeffect = timeeff,
nlfA = treatgroups, nlfB = timepoints,
rho = rho, withinf = "fB", label_list = nameslist)

nset <- seq(7, 14, 2)
cornorm_sim <- simulate_twoway_nrange(cornorm_model, nset, repeated_measurements=TRUE, nsims=5)

##used small number of iterations to reduce computation time

power_results <- test_power_overkn(cornorm_sim, test="rank", plot=TRUE)

Simulated independent and repeated measures two-way experiments over a set of sample sizes

Description

Wrapper for both independent and repeated measures two-way simulations. A vector of defined sample sizes is simulated under the model provided.

Usage

simulate_twoway_nrange(
  matrices_obj,
  nset,
  balanced = TRUE,
  group_size = NULL,
  loss = NULL,
  repeated_measurements = FALSE,
  distribution = "normal",
  skewness = 1,
  shape = 0,
  inferior_limit = -Inf,
  superior_limit = Inf,
  nsims = 200
)

Arguments

Details

For unbalanced independent measures designs, this function generates a simulation with 'max(group_size)' for all combinations of factors and then eliminates observations at random in those factor combinations that have less participants or study subjects. This is also the behavior for unbalanced repeated measures designs when loss="random".

For unbalanced repeated measures designs in which loss="sequential" the participants or subjects from the groups with less observations will be a subset of participants or subjects of groups with more observations. The elimination strategy may not sound like the most efficient way to proceed, is quite fast anyhow.

The 'n' column in the output will reflect how many observations each factor combination has. This should match the input matrix.

Value

List with of data.frames of simulated outcome values under different sample sizes. Each data.frame includes factor level labels, iteration number and sample size.

Examples

refmean <- 1
treatgroups <- 4
timepoints <- 5
treateff <- 1.5
timeeff <- 0.85
factors_levels_names <- list(treatment=letters[1:treatgroups], time=1:timepoints)
## Independent design
effects_treat_time <- calculate_mean_matrix(refmean = refmean,
                                            fAeffect = treateff,fBeffect = timeeff,
                                            nlfA = treatgroups, nlfB = timepoints,
                                            label_list = factors_levels_names)
## Inspect plot to check if matrices correspond to design
effects_treat_time$meansplot
n <- seq(from = 16, to = 24, by = 2)

## In this case, the default 'repeated_measurements', 'distribution' and options are used.
indep_simulation <- simulate_twoway_nrange(effects_treat_time, n)

## Simulate from a truncated distribution
indep_simulation_trunc <- simulate_twoway_nrange(matrices_obj = effects_treat_time, nset = n,
                          distribution="truncated.normal", inferior_limit= 0.8)

##randomly select iteration, select a condition
k <- sample(1:max(indep_simulation_trunc[[1]]$iteration), 1)
toviewdist <- indep_simulation_trunc[[1]]
toviewdist <- subset(toviewdist, iteration==k)
toviewdist <- subset(toviewdist, cond=="V6")
hist(toviewdist$y)

Test simulated two-way factorial design experiments over different sample sizes.

Description

Wrapper to test data simulated under independent or repeated measurements and under different outcome distributions with different sample sizes. Takes output from simulate_twoway_nrange as input, along with test and plotting options.

Usage

test_power_overkn(
  data,
  test = "ANOVA",
  plot = TRUE,
  target_power = NULL,
  title = NULL,
  target_line = TRUE,
  alpha_line = TRUE,
  alpha = 0.05
)

Arguments

Value

Data frame with power and confidence intervals for the main effects and interaction for each of the sample sizes. Also presented in graphical form if plot=TRUE.

Examples

## In this example we simulate an independent sample design with skewed outcome
## Model was specified with the 'calculate_mean_matrix function' (see ?calculate_mean_matrix)
refmean <- 1
treatgroups <- 4
timepoints <- 5
treateff <- 1.25
timeeff <- 0.85
factors_levels_names <- list(treatment=letters[1:treatgroups], time=1:timepoints)

indep_matrix <- calculate_mean_matrix(refmean = refmean,
                                      fAeffect = treateff, fBeffect = timeeff,
                                      nlfA = treatgroups, nlfB = timepoints,
                                      label_list = factors_levels_names)

indep_skewsim <- simulate_twoway_nrange(indep_matrix, seq(6, 12, 2),
                            distribution = "skewed", skewness = 1.8, nsims=5)
##used low number of iterations to reduce computation time

test_power_overkn(indep_skewsim, test="rank")

Simulate measurements repeated over either or both factors of a two-way design

Description

Both regular and internal function. As regular function takes input generated by the calculate_mean_matrix function and iteratively simulates repeated measures two-way factorial experiments. Data are sampled from a normal, skewed normal or truncated normal distribution.

Usage

twoway_simulation_correlated(
  group_size,
  matrices_obj,
  distribution = "normal",
  shape = 0,
  inferior_limit = -Inf,
  superior_limit = Inf,
  balanced = TRUE,
  loss = NULL,
  nsims = 200
)

Arguments

Details

As internal function runs with a single iteration inside graph_twoway_assumptions, which in itself is inside 'calculate_mean_matrix' to generate data for the cell mean and standard deviation plot.

For unbalanced repeated measures designs, this function generates a simulation with max(group_size) for all combinations of factors and then eliminates observations. If loss="random" elimination of in those factor combinations that have less participants or study subjects will occur at random. If loss="sequential" the participants or subjects from the groups with less observations will be a subset of participants or subjects of groups with more observations. This may not sound like the most efficient way to proceed, is quite fast anyhow.

The 'n' column in the output will reflect how many observations each factor combination has. This should match the input matrix.

Value

Dataframe with simulated outcome values, factor level labels and iteration number.

Examples

## Repeated measures design, suppose subjects from 4 independent treatment groups
## measured at 5 different timepoints.

refmean <- 1
treatgroups <- 4
timepoints <- 5
treateff <- 1.5
timeeff <- 0.85
rho <- 0.8
withinf <- "fB"
factors_levels_names <- list(treatment=letters[1:treatgroups], time=1:timepoints)

effects_treat_time <- calculate_mean_matrix(refmean = refmean,
                                            fAeffect = treateff, fBeffect = timeeff,
                                            nlfA = treatgroups,  nlfB = timepoints,
                                            rho = rho, withinf = withinf,
                                            label_list = factors_levels_names)

## Inspect plot to check if matrices correspond to design
effects_treat_time$meansplot

n <- 20
repeatedmeasures_experiment <- twoway_simulation_correlated(group_size = n,
                                matrices_obj = effects_treat_time)

head(repeatedmeasures_experiment, 10)

Simulate independent measurements in a two-way factorial design

Description

Both regular and internal function. As regular function takes input generated by the calculate_mean_matrix function and iteratively simulates independent measures two-way factorial experiments. Outcome may be normally distributed, have a skewed normal distribution or a truncated normal distribution.

Usage

twoway_simulation_independent(
  group_size,
  matrices_obj,
  distribution = "normal",
  skewness = 1,
  inferior_limit = -Inf,
  superior_limit = Inf,
  balanced = TRUE,
  nsims = 200
)

Arguments

Details

As internal function runs with a single iteration inside graph_twoway_assumptions, which in itself is inside calculate_mean_matrix to generate data for the cell mean and standard deviation plot.

For unbalanced independent measures designs, this function generates a simulation with max(group_size) for all factors combinations and then eliminates observations at random in those factor combinations that have less participants or study subjects. This may not sound like the most efficient way to proceed, is quite fast anyhow. The 'n' column in the output will reflect how many observations each factor combination has. This should match the input matrix.

Value

data.frame with modeled outcome values, factor level labels, iteration number and sample size.

Examples

refmean <- 1
treatgroups <- 4
timepoints <- 5
treateff <- 1.5
timeeff <- 0.85
factors_levels_names <- list(treatment=letters[1:treatgroups], time=1:timepoints)

## Independent design
effects_treat_time <- calculate_mean_matrix(refmean = refmean,
                                            fAeffect = treateff, fBeffect = timeeff,
                                            nlfA = treatgroups, nlfB = timepoints,
                                            label_list = factors_levels_names)
## Inspect plot to check if matrices correspond to design
n <- 20
independent_experiment <- twoway_simulation_independent(group_size = n,
                                                        matrices_obj = effects_treat_time)
head(independent_experiment, 10)

Calculate power for global main effects and interaction from two-way factorial simulated data

Description

This functions takes the output of either the twoway_simulation_independent or the twoway_simulation_correlated functions and calculates the power of the sample size used in the simulation under parametric analysis of variance, rank based analysis of variance or permutation testing.

Usage

twoway_simulation_testing(data, test = "ANOVA", alpha = 0.05)

Arguments

Value

A data.frame with the power and 95% confidence interval for each of the main effects and their interaction.

Examples

## After creating a 'matrices_obj' with the 'calculate_mean_matrix' function.

refmean <- 1
treatgroups <- 4
timepoints <- 5
treateff <- 1.5
timeeff <- 0.85
rho <- 0.8
withinf <- "fB"
factors_levels_names <- list(treatment=letters[1:treatgroups], time=1:timepoints)

effects_treat_time <- calculate_mean_matrix(refmean = refmean,
                                            fAeffect = treateff, fBeffect = timeeff,
                                            nlfA = treatgroups,  nlfB = timepoints,
                                            rho = rho, withinf = withinf,
                                            label_list = factors_levels_names)

n <- 7
correlated_sim <- twoway_simulation_correlated(group_size=n, matrices_obj=effects_treat_time,
                  nsims=20)
##used smaller number of iterations to reduce computation time

twoway_simulation_testing(correlated_sim)
## defaults to parametric analysis of variance

twoway_simulation_testing(correlated_sim, test="rank")
## rank based analysis of variance

## permutation test is another option