---
title: "BayesianPower"
author: "Fayette Klaassen"
# date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{BayesianPower}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
## Introduction
*BayesianPower* can be used for sample size determination (using `bayes_sampsize`) and power calculation (using `bayes_power`) when Bayes factors are used to compare an inequality constrained hypothesis $H_i$ to its complement $H_c$, another inequality constrained hypothesis $H_j$ or the unconstrained hypothesis $H_u$. Power is defined as a combination of controlled error probabilities. The unconditional or conditional error probabilities can be controlled. Four approaches to control these probabilities are available in the methods of this package.
**Users are advised to read this vignette and the paper available at <a href = "https://doi.org/10.17605/OSF.IO/D9EAJ">10.17605/OSF.IO/D9EAJ</a> where the four available approaches are presented in detail (Klaassen, Hoijtink \& Gu, unpublished)).**
## Power calculation with `bayes_power()`
`bayes_power(n, h1, h2, m1, m2, sd1=1, sd2=1,
scale = 1000,
bound1 = 1, bound2 = 1/bound1,
datasets = 1000, nsamp = 1000,
seed = 31)`
### Arguments
` n` A number. The sample size for which the error probabilities must be computed.
` h1` A constraint matrix defining H1, see below for more details.
` h2` A constraint matrix defining H2, or a character `'u'` or `'c'` for the unconstrained or
complement hypothesis.
`m1` A vector of expected population means under H1 (standardized), see below for more details.
` m2` A vector of expected populations means under H2 (standardized).
` m2` must be of same length as ` m1`.
` sd1` A vector of standard deviations under H1. Must be a single number (equal
standard deviation under all populations), or a vector of the same length as `m1`
` sd2` A vector of standard deviations under H2. Must be a single number (equal
standard deviation under all populations), or a vector of the same length as `m2`
` scale` A number or use the default `1000` to set the prior scale.
` bound1` A number. The boundary above which BF12 favors H1, see below for more details.
` bound2` A number. The boundary below which BF12 favors H2.
` datasets` A number. The number of datasets to simulate to compute the error probabilities
` nsamp` A number. The number of prior or posterior samples to determine the
complexity or fit.
` seed` A number. The random seed to be set.
### Details
#### Specifying hypotheses
Hypotheses are defined by means of a constraint matrix, that specifies the ordered constraints between the means $\boldsymbol\mu$ using a constraint matrix $R$, such that $R \boldsymbol{\mu} > \bf{0}$, where $R$ is a matrix with $J$ columns and $K$ rows, where $J$ is the number of group means and $K$ is the number of constraints between the means, $\boldsymbol\mu$ is a vector of $J$ means and $\bf{0}$ is a vector of $K$ zeros.
The constraint matrix $R$ contains a set of linear inequality constraints.
Consider
```{r}
R <- matrix(c(1,-1,0,0,1,-1), nrow = 2, byrow = TRUE)
mu <- c(.4, .2, 0)
R
mu
R %*% mu
(R %*% mu) > 0
```
The matrix $R$ specifies that the sum of $1 \times \mu_1$ and $-1 \times \mu_2$ and $0 \times \mu_3$ is larger than $0$, **and** the sum of $0 \times \mu_1$ and $1 \times \mu_2$ and $-1 \times \mu_3$ is larger than $0$. This can also be written as: $\mu_1 > \mu_2 > \mu_3$. For more information about the specification of constraint matrices, see for example [@hoijtink12book].
The argument `h1` has to be a constraint matrix as specified above.
The argument `h2` can be either a constraint matrix, or the character `'u'` or `'c'` if the goal is to compare $H_1$ with $H_u$, the unconstrained hypothesis, or $H_c$ the complement hypothesis.
#### Specifying population means
Hypothesized population means have to be defined under $H_1$ and $H_2$, also if $H_u$ or $H_c$ are considered as $H_2$. The group specific standard deviations can be set under `sd1` and `sd2`, by default, all group standard deviations are $1$.
#### Prior scale
The prior scale can be set using `scale`. By default, a scale of `1000` is used. This implies that the prior covariance matrix is proportional to the standard errors of the sampled data, by a factor of `1000`.
#### Setting bounds
`bound1` and `bound2` describe the boundary used for interpreting a Bayes factor. If `bound1 = 1`, all $BF_{12} > 1$ are considered to express evidence in favor of $H_1$, if `bound1 = 3`, all $BF_{12} > 3$ are considered to express evidence in favor of $H_1$.
Similarly, `bound2` is the boundary *below* which $BF_{12}$ is considered to express evidence in favor of $H_2$.
### Examples
#### Example 1. $H_1$ vs $H_c$
An example where three group means are ordered in $H_1: \mu_1 > \mu_2 > \mu_3$ which is compared to its complement. The power is determined for $n = 40$
```{r, eval = FALSE}
h1 <- matrix(c(1,-1,0,0,1,-1), nrow= 2, byrow= TRUE)
h2 <- 'c'
m1 <- c(.4,.2,0)
m2 <- c(.2,0,.1)
bayes_power(40, m1, m2, h1, h2)
```
#### Example 2. H1 vs H2
An example where four group means are ordered in $H_1: \mu_1 > \mu_2 > \mu_3 > \mu_4$ and in $H_2: \mu_3 > \mu_2 > \ mu_4 > \mu_1$.
Only Bayes factors larger than $3$ are considered evidence in favor of $H_1$ and only Bayes factors smaller than $1/3$ are considered evidence in favor of $H_2$.
```{r, eval = FALSE}
h1 <- matrix(c(1,-1,0,0,0,1,-1,0,0,0,1,-1), nrow= 3, byrow= TRUE)
h2 <- matrix(c(0,-1,1,0,0,1,0,-1,-1,0,0,1), nrow = 3, byrow= TRUE)
m1 <- c(.7,.3,.1,0)
m2 <- c(0,.4,.5,.1)
bayes_power(34, h1, h2, m1, m2, bound1 = 3, bound2 = 1/3)
```
## Sample size determination with `bayes_sampsize()`
`
bayes_sampsize(h1, h2, m1, m2, sd1 = 1, sd2 = 1, scale = 1000,
type = 1, cutoff, bound1 = 1, bound2 = 1 / bound1,
datasets = 1000, nsamp = 1000,
minss = 2, maxss = 1000, seed = 31)
`
### Arguments
The arguments are the same as for `bayes_power()` with the addition of:
`type`A character. The type of error to be controlled. The options are: `"1", "2", "de", "aoi", "med.1", "med.2"`. See below for more details.
`cutoff` A number. The cutoff criterion for type.
If `type` is `"1", "2", "de", "aoi"`, `cutoff` must be between $0$ and $1$.
If `type` is `"med.1"` or `"med.2"`, `cutoff` must be larger than $1$. See below for more details.
`minss` A number. The minimum sample size.
`maxss` A number. The maximum sample size.
### Details
`bayes_sampsize()` iteratively uses `bayes_power()` to determine the error probabilities for a sample size, evaluates whether the chosen error is below the cutoff, and adjusts the sample size.
#### `type`
[@klaassenPIH] describes in detail the different types of controlling error probabilities that can be considered. Specifying `"1"` or `"2"` indicates that the Type 1 or Type 2 error probability has to be controlled, respectively the probability of concluding $H_2$ is the best hypothesis when $H_1$ is true or concluding that $H_1$ is the best hypothesis when $H_2$ is true. Note that when $H_1$ or $H_2$ is considered the best hypothesis depends on the values chosen for `bound1` and `bound2`. Specifying `"de"` or `"aoi" ` indicates that the Decision error probability (average of Type 1 and Type 2) or the probability of Indecision has to be controlled. Finally, specifying `" med.1"` or `"med.2"` indicates the minimum desired median $BF_{12}$ when $H_1$ is true, or the minimum desired median $BF_{21}$ when $H_2$ is true.
### Examples
#### Example 1. $H_1$ versus $H_c$, controlling decision error
```{r, eval = FALSE}
h1 <- matrix(c(1, -1, 0,
0, 1, -1),
nrow= 2, byrow= TRUE)
h2 <- 'c'
m1 <- c(.4, .2, 0)
m2 <- c(.2, 0, .1)
bayes_sampsize(h1, h2, m1, m2, type = "de", cutoff = .125)
```
#### Example 2. H1 versus H2, controlling indecision error
```{r, eval = FALSE}
h1 <- matrix(c(1, -1, 0, 0,
0, 1, -1, 0,
0, 0, 1, -1),
nrow= 3, byrow= TRUE)
h2 <- matrix(c(0, -1, 1, 0,
0, 1, 0, -1,
-1, 0, 0, 1),
nrow = 3, byrow= TRUE)
m1 <- c(.7, .3, .1, 0)
m2 <- c(0, .4, .5, .1)
bayes_sampsize(h1, h2, m1, m2, type = "aoi", cutoff = .2, minss = 2, maxss = 500)
```
#### Example 3. $H_1$ versus $H_u$, controlling median Bayes factor
```{r, eval = FALSE}
h1 <- matrix(c(1, -1, 0, 0,
0, 1, -1, 0,
0, 0, 1, -1),
nrow= 3, byrow= TRUE)
h2 <- 'u'
m1 <- c(.3, .2, 0)
m2 <- c(0, 0, 0)
bayes_sampsize(h1, h2, m1, m2, type = "med.1", cutoff = 3, minss = 2, maxss = 500)
```
## References
Hoijtink, H. (2012). *Informative hypotheses. Theory and practice for behavioral and social scientists.* Boca Raton: Chapman Hall/CRC.
Klaassen, F., Hoijtink, H., Gu, X. (unpublished). *The power of informative hypotheses.* Pre-print available at <a href = "https://doi.org/10.17605/OSF.IO/D9EAJ">https://doi.org/10.17605/OSF.IO/D9EAJ</a>