if (FALSE) {
# TASK: Find specific sample size in each group for independent t-test
# corresponding to a power rate of .8
#
# For ease of the setup, assume the groups are the same size, and the mean
# difference corresponds to Cohen's d values of .2, .5, and .8
# This example can be solved numerically using the pwr package (see below),
# though the following simulation setup is far more general and can be
# used for any generate-analyse combination of interest
# SimFunctions(SimSolve=TRUE)
#### Step 1 --- Define your conditions under study and create design data.frame.
#### However, use NA placeholder for sample size as it must be solved,
#### and add desired power rate to object
Design <- createDesign(N = NA,
d = c(.2, .5, .8),
sig.level = .05)
Design # solve for NA's
#~~~~~~~~~~~~~~~~~~~~~~~~
#### Step 2 --- Define generate, analyse, and summarise functions
Generate <- function(condition, fixed_objects = NULL) {
Attach(condition)
group1 <- rnorm(N)
group2 <- rnorm(N, mean=d)
dat <- data.frame(group = gl(2, N, labels=c('G1', 'G2')),
DV = c(group1, group2))
dat
}
Analyse <- function(condition, dat, fixed_objects = NULL) {
p <- t.test(DV ~ group, dat, var.equal=TRUE)$p.value
p
}
Summarise <- function(condition, results, fixed_objects = NULL) {
# Must return a single number corresponding to f(x) in the
# root equation f(x) = b
ret <- EDR(results, alpha = condition$sig.level)
ret
}
#~~~~~~~~~~~~~~~~~~~~~~~~
#### Step 3 --- Optimize N over the rows in design
# Initial search between N = [10,500] for each row using the default
# integer solver (integer = TRUE)
# In this example, b = target power
solved <- SimSolve(design=Design, b=.8, interval=c(10, 500),
generate=Generate, analyse=Analyse,
summarise=Summarise)
solved
summary(solved)
plot(solved, 1)
plot(solved, 2)
plot(solved, 3)
# also can plot median history and estimate precision
plot(solved, 1, type = 'history')
plot(solved, 1, type = 'density')
# verify with true power from pwr package
library(pwr)
pwr.t.test(d=.2, power = .8, sig.level = .05)
pwr.t.test(d=.5, power = .8, sig.level = .05)
pwr.t.test(d=.8, power = .8, sig.level = .05)
# use estimated N results to see how close power was
N <- solved$N
pwr.t.test(d=.2, n=N[1], sig.level = .05)
pwr.t.test(d=.5, n=N[2], sig.level = .05)
pwr.t.test(d=.8, n=N[3], sig.level = .05)
# with rounding
N <- ceiling(solved$N)
pwr.t.test(d=.2, n=N[1], sig.level = .05)
pwr.t.test(d=.5, n=N[2], sig.level = .05)
pwr.t.test(d=.8, n=N[3], sig.level = .05)
# failing analytic formula, confirm results with more precise
# simulation via runSimulation()
csolved <- solved
csolved$N <- ceiling(solved$N)
confirm <- runSimulation(design=csolved, replications=10000, parallel=TRUE,
generate=Generate, analyse=Analyse,
summarise=Summarise)
confirm
############
# Similar setup as above, however goal is now to solve d given sample
# size and power inputs (inputs for root no longer required to be an integer)
Design <- createDesign(N = c(100, 50, 25),
d = NA,
sig.level = .05)
Design # solve for NA's
#~~~~~~~~~~~~~~~~~~~~~~~~
#### Step 2 --- Define generate, analyse, and summarise functions (same as above)
#~~~~~~~~~~~~~~~~~~~~~~~~
#### Step 3 --- Optimize d over the rows in design
# search between d = [.1, 2] for each row
# In this example, b = target power
# note that integer = FALSE to allow smooth updates of d
solved <- SimSolve(design=Design, b = .8, interval=c(.1, 2),
generate=Generate, analyse=Analyse,
summarise=Summarise, integer=FALSE)
solved
summary(solved)
plot(solved, 1)
plot(solved, 2)
plot(solved, 3)
# verify with true power from pwr package
library(pwr)
pwr.t.test(n=100, power = .8, sig.level = .05)
pwr.t.test(n=50, power = .8, sig.level = .05)
pwr.t.test(n=25, power = .8, sig.level = .05)
# use estimated d results to see how close power was
pwr.t.test(n=100, d = solved$d[1], sig.level = .05)
pwr.t.test(n=50, d = solved$d[2], sig.level = .05)
pwr.t.test(n=25, d = solved$d[3], sig.level = .05)
# failing analytic formula, confirm results with more precise
# simulation via runSimulation()
confirm <- runSimulation(design=solved, replications=10000, parallel=TRUE,
generate=Generate, analyse=Analyse,
summarise=Summarise)
confirm
}
Run the code above in your browser using DataLab