BESTpower {BEST}R Documentation

Estimating statistical power

Description

Estimation of the probability of meeting the goals of a study given initial information or assumptions about the population parameters. For prospective power estimation, the sequence
makeData -> BESTmcmc -> BESTpower
is recommended: see makeData.

Usage

BESTpower(BESTobj, N1, N2, credMass=0.95,
 ROPEm, ROPEsd, ROPEeff, 
 maxHDIWm, maxHDIWsd, maxHDIWeff,
 compValm = 0, nRep = 200, mcmcLength = 10000,
 saveName = NULL, showFirstNrep = 0, verbose = 2, rnd.seed=NULL, parallel=NULL)

Arguments

BESTobj

an object of class BEST produced by BESTmcmc.

N1

planned sample size for the first (or only) group of observations; may be a scalar if sample size is fixed, or a vector if sample size varies; values will be recycled if necessary.

N2

planned sample size for the second group of observations; ignored if BESTobj concerns only one group.

credMass

the probability mass to include in HDIs when checking criteria.

ROPEm

a two element vector, such as c(-1, 1), specifying the limit of the ROPE on the difference of means (for 2 groups) or the mean (for 1 group).

ROPEsd

a two element vector, such as c(-1, 1), specifying the limit of the ROPE on the (difference of) standard deviations.

ROPEeff

a two element vector, such as c(-1, 1), specifying the limit of the ROPE on the effect size.

maxHDIWm

the maximum acceptable width for the HDI for the difference in means (for 2 groups) or for the mean (for a single group).

maxHDIWsd

the maximum acceptable width for the HDI for the (difference of) standard deviation.

maxHDIWeff

the maximum acceptable width for the HDI for the effect size.

compValm

for a single group, the value of the mean which represents no effect; used to calculate the effect size. Ignored for 2 groups.

nRep

number of simulations to carry out.

mcmcLength

length of the MCMC chains to use for each simulation.

saveName

if required, the results may saved to a file after each iteration and saveName specifies the file name (or path relative to the current working directory) to use. The power object can be loaded with load. Set to NULL (the default) to disable saving.

showFirstNrep

the number of results to display as plots at the beginning of the simulation run. (This uses dev.new(), which does not work in Rstudio. The plots will appear sequentially in the plot window and you will have to use the back arrow to review them.)

verbose

controls output to the R Console: 0 suppresses all output; 1 gives just a progress bar; 2 gives maximum detail.

rnd.seed

a positive integer (or NULL): the seed for the random number generator, used to obtain reproducible samples if required.

parallel

if NULL or TRUE and > 3 cores are available, the MCMC chains are run in parallel. (If TRUE and < 4 cores are available, a warning is given.)

Details

For each of the parameters of interest - (difference in) mean, (difference in) standard deviation and effect size - we consider 4 criteria and the probability that each will be met:

1. The HDI of the posterior density of the parameter lies entirely outside the ROPE and is greater than the ROPE.

2. The HDI of the posterior density of the parameter lies entirely outside the ROPE and is less than the ROPE.

3. The HDI of the posterior density of the parameter lies entirely inside the ROPE.

4. The width of the HDI is less than the specified maxHDIWx.

The mass inside the above HDIs depends on the credMass argument.

A uniform beta prior is used for each of these probabilities and combined with the results of the simulations to give a conjugate beta posterior distribution. The means and 95% HDI credible intervals are returned.

Value

A matrix with a row for each criterion and columns for the mean and lower and upper limits of a 95% credible interval for the posterior probability of meeting the criterion.

Note that this matrix always has 12 rows. Rows corresponding to criteria which are not specified will have NAs.

Note

At least 1000 simulations are needed to get good estimates of power and these can take a long time. If the run is interrupted, the results so far can be recovered from the file specified in saveName.

The chains in BESTobj must have at least nRep values. To allow for some degree of autocorrelation among values, it would be prudent to make these chains at least 10 * nRep in length.

Author(s)

Original code by John Kruschke, modified by Mike Meredith.

References

Kruschke, J. K. 2013. Bayesian estimation supersedes the t test. Journal of Experimental Psychology: General 142(2):573-603. doi: 10.1037/a0029146

Kruschke, J. K. 2011. Doing Bayesian data analysis: a tutorial with R and BUGS. Elsevier, Amsterdam, Chapter 13.

See Also

makeData for details of preparing a BESTobj for a prospective power analysis.

Examples


## For retrospective power analysis, see the example in BEST-package.

# 1. Generate idealized data set:
proData <- makeData(mu1=108, sd1=17, mu2=100, sd2=15, nPerGrp=20, 
                         pcntOut=10, sdOutMult=2.0, rnd.seed=NULL)

# 2. Generate credible parameter values from the idealized data:
proMCMC <- BESTmcmc(proData$y1, proData$y2, numSavedSteps=2000, parallel=FALSE)  

# 3. Compute the prospective power for planned sample sizes:
# We'll  do just 5 simulations to show it works; should be several hundred.
N1plan <- N2plan <- 50
powerPro <- BESTpower(proMCMC, N1=N1plan, N2=N2plan,
               ROPEm=c(-1.5,1.5), ROPEsd=c(-2,2), ROPEeff=c(-0.5,0.5), 
               maxHDIWm=15.0, maxHDIWsd=10.0, maxHDIWeff=1.0, nRep=5, parallel=FALSE)
powerPro


[Package BEST version 0.5.4 Index]