---
title: "pwr4exp: Power Analysis for Experimental Designs"
author: "R Package Version `r packageVersion('pwr4exp')`"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{pwr4exp: Power Analysis for Experimental Designs}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---

```{r, include = FALSE, message=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(pwr4exp)
```

# Introduction

**pwr4exp** is an R package for performing statistical power analysis for experiments analyzed using linear mixed models (LMM). It supports power calculations for omnibus F-tests (ANOVA-like tests) and specific contrasts (t-tests) for Gaussian response variables, employing the Satterthwaite approximation for degrees of freedom.

**Key features:**

-   **Consistency:** Specify the model formula to ensure consistency between power analysis and data analysis.
-   **Flexibility:** Support a variety of experimental designs through either the general-purpose function (`mkdesign`) or specialized functions (e.g., `designCRD`, `designRCBD`, `designLSD`, `designCOD`, `designSPD`) for standard designs.
-   **Correlation structures:** Allow specification of various residual variance-covariance structures, which is essential for analyzing repeated measures, spatial, or other correlated data.

The workflow involves two main steps:

1.  [Creating a design object](#design): Defining the experimental structure and specifying parameters (fixed effects, (co-)variances).

2.  [Calculating power](#power): Computing the power of F-tests or t-tests for the effects of interest.

# Creating a design object {#design}

A design object in **pwr4exp** is a `list` containing design matrices, fixed effects, variance-covariance components, and other higher-level parameters created by the design-generating functions. It is not an experimental design in the traditional sense, where randomization and unit allocation are performed, but rather a container for all the information necessary to conduct power analysis.

Two approaches are provided to define a design:

-   [**General approach with `mkdesign`:**](#mkdesign) A design object can be created by providing a model formula and a data frame that outlines the experimental layout.
-   [**Standard designs with Specific Functions:**](#standard_design) Functions such as `designCRD`, `designRCBD`, `designLSD`, `designCOD`, or `designSPD` generate design objects for common experimental setups. These functions allow specification of design characteristics (e.g., number of treatments, blocks, or replicates) and automatically handle the creation of the data frame.

## The `mkdesign` function {#mkdesign}

The function `mkdesign` is versatile for designs as long as the data frame is provided.

**Usage**

``` r
mkdesign(
  formula, data,
  beta = NULL, means = NULL,
  vcomp = NULL, sigma2 = NULL,
  correlation = NULL,
  template = FALSE,
  REML = TRUE
)
```

**Arguments**

-   **`formula`:** A right-hand-side model formula that defines fixed and random effects. (Use `lme4::lmer` syntax for specifying random effects.)

-   **`data`:** A data frame containing all independent variables required by the model, consistent with the design's structure.

-   **`means`** or **beta:** Expected fixed effects. Either a vector of model coefficients (`beta`) or expected means (`means`) must be provided. [Templates](#fixeff_template) for parameter ordering can be generated by setting `template = TRUE`.

-   **`vcomp`:** Variance components for any random effects present in the model, provided as a numeric vector. A [template](#vcomp_template) for input ordering can be generated by setting `template = TRUE`.

-   **`sigma2`:** The residual variance.

-   **`correlation`:** (Optional) A correlation structure (`nlme::corClasses`) for repeated measures.

-   [**`template`:**](#template) When set to `TRUE` or when only the `formula` and `data` arguments are provided, templates for `beta`, `means`, and `vcomp` are generated to indicate the required input order.

-   **`REML`:** (logical) Whether to use the REML or ML information matrix. The default is `TRUE` (REML).

## Specific design functions {#standard_design}

While `mkdesign` provides full flexibility, **pwr4exp** also offers specific functions for common experimental designs. These functions allow users to define a design based on treatment structure and replications without manually creating a data frame.

> **Completely randomized design (CRD)**

``` r
designCRD(treatments, label, replicates, formula, beta, means, sigma2, template = FALSE)
```

> **Randomized complete block design (RCBD)**

``` r
designRCBD(treatments, label, blocks, formula, beta, means, vcomp, sigma2, template = FALSE)
```

> **Latin square design (LSD)**

``` r
designLSD(treatments, label, squares, reuse, formula, beta, means, vcomp, sigma2, template = FALSE)
```

> **Crossover design**

This is a special case of LSD where time periods and individuals act as blocks. Period blocks are reused when replicating squares.

``` r
designCOD(treatments, label, squares, formula, beta, means, vcomp, sigma2, template = FALSE)
```

> **Split-plot design (SPD)**

``` r
designSPD(trt.main, trt.sub, label, replicates, formula, beta, means, vcomp, sigma2, template = FALSE)
```

The inputs for these functions are similar to `mkdesign`, except that `data` is replaced by predefined design characteristics such as treatment levels and replications.

**Key inputs**

> -   **`treatments`**: A numeric vector specifying the number of levels for treatment factors. Multiple factors are arranged factorially. For example, `treatments = c(3, 2)` specifies two treatment factors one with 3 levels and another with 2 levels, forming a factorial treatment design.
>
> -   **`trt.main`** and **`trt.sub`**: Define main-plot and subplot treatments for SPD, following the same rule for `treatment`.
>
> -   **`label`**: (optional) A list assigning labels to factor levels. If not specified, default names are assigned. For one treatment factor, the default is `list(trt = c("1", "2", ...))`. For two factors, the default is `list(facA = c("1", "2", ...), facB = c("1", "2", ...))`, where "facA" and "facB" represent the two factors, and "1", "2", etc., represent the levels of each factor.
>
> -   **`replicates`**: Number of experimental units per treatment in CRD or number of main plots (i.e., the number of experimental units per main-plot treatment) in SPD.
>
> -   **`blocks`**: Number of blocks in RCBD.
>
> -   **`squares`**: Number of squares in LSD or COD.
>
> -   **`reuse`**: Specifies how squares are replicated in LSD. One of:
>
> > -   `"row"`: reuse row blocks
> > -   `"col"`: reuse column blocks
> > -   `"both"`: reuse both row and column blocks
>
> -   [**`template`**](#template): When set to `TRUE`, templates for `beta`, `means`, and `vcomp` are generated to indicate the required input order.

Each of these design-generating functions has a default `formula` based on the treatment structure (e.g., one factor or factorial factors). If `formula` is not specified, a default formula with main effects and all interactions (if applicable) is used internally. In RCBD, LSD, COD, and SPD designs, all block factors are fitted as random effects. The `formula` component in the output design object (a `list`) can be inspected.

## Parameter templates {#template}

Templates help ensure the correct ordering of fixed effects (`beta` or `means`) and random effects (`vcomp`).

Below, an exemplary data frame is created to illustrate these templates. Note that this example **does NOT** represent a realistic design. The data frame includes:

> -   Four categorical variables (`fA`, `fB`, `fC`, `fD`),
> -   Two numerical variables (`x`, `z`).

```{r}
 df1 <- expand.grid(
   fA = factor(1:2), # factor A with 2 levels
   fB = factor(1:2), # factor B with 2 levels
   fC = factor(1:3), # factor C with 3 levels
   fD = factor(1:3), # factor D with 3 levels
   subject = factor(1:10)  # 10 subjects
 )
 df1$x <- rnorm(nrow(df1))  # Numerical variable x
 df1$z <- rnorm(nrow(df1))  # Numerical variable z
```

### `beta` Template {#fixeff_template}

> The `beta` template displays the order of model coefficients in the fitted model.

This is particularly useful when specifying expected model coefficients directly as effect size measures and when the model includes complex interactions. For example, the expected values of the model coefficients for `~ fA*fB + x` should be provided in the following sequence:

```{r}
mkdesign( ~ fA * fB + x, df1)$fixeff$beta
```

### `means` Template

For treatment factors, it may be more more convenient to provide `means` instead of `beta`. The `means` template represents either marginal or cell means for factors, depending on the presence of interactions. For numerical variables included in the model, regression coefficients remain necessary alongside alongside the means for factors.

**Factors without interactions**

> If no interactions are present, **marginal means** for each factor level are required.

For example, for the model `~ fA + fB`, means should be reported in the order specified by the template:

```{r}
mkdesign(~ fA + fB, df1, template = TRUE)$fixeff$means
```

**Interactions among factors**

> When factors interact, **conditional (cell) means** must be provided for all possible combinations of their levels.

For example, for a model with three interacting factors (`fA`, `fB`, `fC`), cell means are required for all 12 level combinations in the following order:

```{r}
mkdesign(~ fA * fB * fC, df1)$fixeff$means
```

**Numerical variables**

> For numerical variables, regression coefficients are required. If the model only includes numerical variables, the intercept must also be included. In this case, the values in means are identical to beta:

```{r}
mkdesign(~ x + z, df1)$fixeff$means
```

> When numerical variables interact, regression coefficients for both main effects and interaction terms must be included.

```{r}
mkdesign(~ x * z, df1)$fixeff$means
```

**Factor-by-numerical interactions**

> For interactions between factor and numerical variables, both **Marginal means** for the factor, and **Regression coefficients** for the numerical variable at each level of the factor must be provided.

For instance, in the model `~ fA * x`, the means for levels of `fA` and the regression coefficients for `x` at each level of `fA` must be provided:

```{r}
mkdesign(~ fA * x, df1)$fixeff$means
```

**Combining Multiple Situations**

> The rules outlined for factor and numerical variables apply when multiple types of variables and interactions appear in the same model.

For example, consider a model including interactions among three factors (`fA`, `fB`, `fC`), factor-by-numerical interaction (`fD * x`), and main effects for another numerical variable (`z`).

```{r}
mkdesign(~ fA * fB * fC + fD * x + z, df1)$fixeff$means
```

The required elements and their order in means are:

1\. Means for each level of `fD` (positions 1–3)

2\. Regression coefficient for `z` (position 4)

3\. Regression coefficients for `x` at each level of `fD` (positions 5–7)

4\. Cell means for all combinations of levels of `fA`, `fB`, and `fC` (positions 8–19)

### `vcomp` Template {#vcomp_template}

> The `vcomp` template represents variance-covariance matrices, where integer values indicate the order of variance components in the input vector.

For a single random effect, a template is unnecessary since it corresponds to a single variance value. For multiple random effects, the template specifies the sequence of variance and covariance components, ensuring proper alignment when specifying variance components in the model.

For example, consider a model that includes both a random intercept and a random slope for `fA` by `subject`:

```{r}
mkdesign(~ fA * fB * fC * fD + (1 + fA | subject), df1)$varcov
```

The template specifies the following required inputs:

1.  Variance of the random intercept

2.  Covariance between the random intercept and slop

3.  Variance of the random slop (`fA2`)

## Correlation structures

Although specifying complex variance-covariance structures, such as unstructured `nlme::corSymm`, may seem impractical for power analysis, various correlation structures from the **nlme** package can be applies.

The available correlation structures, as outlined in [`nlme::corClasses`](https://www.rdocumentation.org/packages/nlme/versions/3.1-167/topics/corClasses), including

-   [`corAR1`](https://www.rdocumentation.org/packages/nlme/versions/3.1-163/topics/corAR1) – First-order autoregressive

-   [`corARMA`](https://www.rdocumentation.org/packages/nlme/versions/3.1-163/topics/corARMA) – Autoregressive moving average

-   [`corCAR1`](https://www.rdocumentation.org/packages/nlme/versions/3.1-163/topics/corCAR1) – Continuous-time autoregressive

-   [`corCompSymm`](https://www.rdocumentation.org/packages/nlme/versions/3.1-163/topics/corCompSymm) – Compound symmetry

-   [`corExp`](https://www.rdocumentation.org/packages/nlme/versions/3.1-163/topics/corExp) – Exponential spatial correlation

-   [`corGaus`](https://www.rdocumentation.org/packages/nlme/versions/3.1-163/topics/corGaus) – Gaussian spatial correlation

-   [`corLin`](https://www.rdocumentation.org/packages/nlme/versions/3.1-163/topics/corLin) – Linear spatial correlation

-   [`corSymm`](https://www.rdocumentation.org/packages/nlme/versions/3.1-163/topics/corSymm) – Unstructured correlation

-   [`corRatio`](https://www.rdocumentation.org/packages/nlme/versions/3.1-163/topics/corRatio) – Ratio-based spatial correlation

-   [`corSpher`](https://www.rdocumentation.org/packages/nlme/versions/3.1-163/topics/corSpher) – Spherical spatial correlation

Note: In `nlme::corAR1()` and `nlme::corARMA()` when `p=1` and `q=0`, the time variable must be an integer. However, in `pwr4exp`, this restriction has been released, factors are also supported.

# Power calculation {#power}

Once a design object has been created, calculating statistical power is straightforward.

## Power Calculation for Omnibus F-Tests

The statistical power for **F-tests** (ANOVA-like tests) can be computed using the `pwr.anova` function.

**Usage**

``` r
pwr.anova(object, sig.level = 0.05, type = 3)
```

**Arguments**

-   **`object`**: the design object created in the previous step.

-   **`sig.level`**: significance level, default `0.05`.

-   **`type`**: the type of ANOVA table (default: `3`).

## Power Calculation for Specific Contrasts

To evaluate statistical power for **t-tests** on specific contrasts, use the `pwr.contrast` function.

**Usage**

``` r
pwr.contrast(object, which, by = NULL, 
  contrast = c("pairwise", "poly", "trt.vs.ctrl"),
  sig.level = 0.05, p.adj = FALSE, 
  alternative = c("two.sided", "one.sided"), strict = TRUE
)
```

**Arguments**

-   **`object`**: The design object.

-   **`which`**: The factor of interest. Multiple factors can be combined using `:` or `*` (e.g., `"facA*facB"` treats combined factor levels as a single factor).

-   **`by`**: The variable to condition on.

-   **`contrast`**: Contrast method, on of

    -   `"pairwise"` - pairwise comparisons

    -   `"poly"` - polynomial contrasts

    -   `"trt.vs.ctrl"` - treatment vs. control comparison

    -   Alternatively, a **numeric vector** for a single custom contrast, or a **named list** of contrast vectors. If a list is provided, each vector must match the number of factor levels in each `by` group. In multi-factor scenarios, factor levels are combined and treated as a single factor.

-   **`sig.level`**: significance level (default: `0.05`).

-   **`p.adj`**: whether the `sig.level` should be adjusted using the Bonferroni method (default: FALSE).

-   **`alternative`**: `"two.sided"` (default) or `"one.sided"`.

-   **`strict`**: Controls how power is calculated in two-sided tests (default: `TRUE`). If `TRUE`, includes the probability of rejection in the opposite direction of the true effect. If `FALSE`, power equals half the significance level when the true difference is zero.

## Power Calculation for Model Coefficients

The `pwr.summary` function computes the statistical power for testing model coefficients (t-tests).

**Usage**

``` r
pwr.summary(object, sig.level = 0.05)
```

**Arguments**

-   **`object`**: A design object
-   **`sig.level`**: The significance level (default: 0.05)

# Practical Examples

## Example 1. Completely Randomized Design

In this example, we create a CRD with one treatment factor.

**Design parameters**

1.  **Treatments**: 1 treatment factor with 4 levels.
2.  **Replicates**: 8 experimental units per treatment.
3.  **Mean and effect size**: Expected means for four levels: 35, 30, 37, 38
4.  **Error variance**: 15

**Step 1: Creating the CRD**

```{r}
crd <- designCRD(
  treatments = 4,
  replicates = 8,
  means = c(35, 30, 37, 38),
  sigma2 = 15
)
```

**Step 2: Power for omnibus test**

To compute the power for the **overall F-test (ANOVA)**:

```{r}
pwr.anova(crd)
```

**Step 3: Power for specific contrasts**

> **a) Pairwise Comparisons**

To compute power for pairwise differences:

```{r}
pwr.contrast(crd, which =  "trt", contrast = "pairwise")
```

> **b) Polynomial Contrasts**

To test for trends across treatment levels:

```{r}
pwr.contrast(crd, which =  "trt", contrast = "poly")
```

> **c) Treatment vs. Control Comparison**

The power for detecting differences between treatments and the control:

```{r}
pwr.contrast(crd, which =  "trt", contrast = "trt.vs.ctrl")
```

> **d) Custom Contrast: All Treatments vs. Control**

In addition to pre-defined contrast method, custom contrast vectors can be specified. For example, to compare the average of all treatments against the control, use the contrast vector `c(-1, 1/3, 1/3, 1/3)`. The power for this test can be computed as:

```{r}
pwr.contrast(crd, which =  "trt", contrast = list(trts.vs.ctrl = c(-1, 1/3, 1/3, 1/3)))
```

**Adjusting for Multiple Comparisons**

In real-world analysis, **P-values** often need **adjustment for multiple comparisons (MCP)**. While most **post-hoc** methods cannot be applied **directly in power analysis**, we can **lower the significance level** to approximate MCP adjustments.

> **Using a More Conservative Significance Level**

```{r}
pwr.contrast(crd, which =  "trt", contrast = "pairwise", sig.level = 0.01)
```

> **Using the Bonferroni Correction**
>
> The `pwr.contrast` function includes an option for **Bonferroni correction**, though it may be overly conservative:

```{r}
pwr.contrast(crd, which =  "trt", contrast = "pairwise", sig.level = 0.05, p.adj = TRUE)
```

## Example 2. Randomized complete block design

In this example, we create a **Randomized Complete Block Design (RCBD)** with **two treatment factors** in a **2×2 factorial design**.

**Design Parameters**

1.  **Treatments**: A 2x2 factorial design (factors: `A` and `B`).
2.  **Replicates**: 8 blocks.
3.  **Expected Means**: $$
        \begin{array}{c|c|c}
     & B1 & B2 \\
        \hline
        A1 & 35 & 38 \\
        A2 & 40  & 41 \\
        \end{array}
        $$ Corresponding `beta` values are as follows (Optional):

-   *Intercept (`A1B1`): 35*
-   *Effect of `A2` alone: 5 units*
-   *Effect of `B2` alone: 3 units*
-   *Interaction (`A2:B2`): -2 units (i.e., the combined effect of `A2` and `B2` is 2 units below the additive effects).*

4.  **Variance among blocks**: 11.
5.  **Error variance**: 4.
6.  The total variance of the response variable (15) is decomposed into variance between blocks (11) and variance within blocks (4).

**Step 1: Creating the RCBD**

-   **Generate the templates** to determine the correct parameter order. Provide `beta` or `means`, and `vcomp` according to the order below.

```{r}
designRCBD(treatments = c(2, 2), blocks = 8, template = TRUE)
```

-   **Create the design object**

```{r}
rcbd <- designRCBD(
  treatments = c(2, 2),
  blocks = 8,
  # beta = c(35, 5, 3, -2), # identical to means
  means = c(35, 40, 38, 41),
  vcomp = 11,
  sigma2 = 4
)
```

The default factor names and formula can be inspected in the design object:

```{r}
unique(rcbd$deStruct$fxTrms$fixedfr)
rcbd$deStruct$formula
```

Treatment factors and factor levels can be renamed and relabeled, and the model formula can be changed. For example, if interactions are removed, marginal means for both factors must be provided instead of cell means.

```{r}
designRCBD(treatments = c(2, 2), 
           label = list(factorA = c("A1", "A2"), factorB = c("B1", "B2")), 
           blocks = 8, 
           formula = ~ factorA + factorB + (1|block), 
           template = TRUE)
```

**Step 2: Evaluate statistical power**

To test overall effects of `facA` and `facB`

```{r}
pwr.anova(rcbd)
```

To test differences between levels of `facA`, conditioned on `facB`:

```{r}
pwr.contrast(rcbd, which = "facA", by = "facB")
```

## Example 3. Latin square design

In this example, we extend *Example 2* by introducing **another blocking factor**, transforming the design into an LSD.

-   The **treatment structure and effect sizes remain the same** as in *Example 2*.

-   The LSD controls for **two sources of variability** (row and column blocks) while evaluating treatment effects.

-   The **total variance (15)** is further decomposed into:

    -   **Variance between row blocks**: `11`

    -   **Variance between column blocks**: `2`

    -   **Residual error variance**: `2`

**Step 1: Creating the LSD**

Generate the templates to determine the correct order of inputs.

```{r}
designLSD(
  treatments = c(2, 2),
  squares = 4,
  reuse = "both",
  template = TRUE
)
```

Either `beta` or `means` can be provided, as in *Example 2*. Note that although the variance component template sets the order for `row` and `col` variances in `vcomp`, these values don’t impact power as long as `sigma2` is fixed because treatment effects are tested against residual error. However, in designs where the error terms for treatment factors include variance components beyond residual error (like the next example), variance components can affect power for main plot factors. It's also recommended to use variance components as a guide to make an informed estimate of `sigma2`.

**Create the LSD**

```{r}
lsd <- designLSD(
  treatments = c(2, 2),
  label = list(temp = c("T1", "T2"), dosage = c("D1", "D2")),
  squares = 4,
  reuse = "both",
  means = c(35, 40, 38, 41),
  vcomp = c(11, 2),
  sigma2 = 2
)
```

Once the design has been created, `pwr.anova` and `pwr.contrast` can be used to evaluate statistical power as demonstrated above.

## Example 4: Split-plot Design

In this example, we create an SPD with two treatment factors, one at the main plot level and another at the subplot level. The design parameters are as follows:

1.  **Treatments:**

    -   **Main plot factor** with **2 levels**

    -   **Subplot factor** with **3 levels**

2.  **Replicates**:

    -   Each main plot treatment has **5 plots**, for a total of **10 plots**.

    -   This is a standard SPD, where the plots are the blocks at the subplot level and the subplot design follows an RCBD structure.

3.  **Expected cell means** are: $$
        \begin{array}{c|c|c}
          & trt.sub1 & trt.sub2 & trt.sub3 \\
        \hline
        trt.main1 & 20 & 22 & 24\\
        trt.main2 & 22 & 24 & 28 \\
        \end{array}
        $$

4.  **Variance Components**:

    -   **Main-plot error**: `4`

    -   **Subplot (residual) error**: `11`

    -   **Total variance**: `15`

**Generate input templates**

```{r}
designSPD(
  trt.main = 2,
  trt.sub = 3, 
  replicates = 10, 
  template = TRUE
)
```

**Create the SPD**

```{r}
spd <- designSPD(
  trt.main = 2,
  trt.sub = 3, 
  replicates = 10, 
  means = c(20, 22, 22, 24, 24, 28),
  vcomp = 4,
  sigma2 = 11
)
```

One can verify that different values of **main-plot error (`vcomp`)** affect the power of the **main plot factor**, unlike in *Example 3*, where variance components did not influence power.

## Example 5: Repeated measures

This example illustrates a **repeated measures design,** where **three treatments** (`CON`, `TRT1`, `TRT2`) are measured **hourly over 8 hours**. **Within-subject correlations** are modeled using an **AR(1) structure** with $\rho = 0.6$ and $\sigma^2 = 2$.

**Design Details**

1\. **Subjects**: 6 per treatment group (total: 18 subjects).\
2. **Treatments**: `CON`, `TRT1`, and `TRT2`.\
3. **Time Points**: 8 hourly measurements.

4\. **Means**:$$
\begin{array}{c|cccccccc}
\textbf{Treatment} & \textbf{1} & \textbf{2} & \textbf{3} & \textbf{4} & \textbf{5} & \textbf{6} & \textbf{7} & \textbf{8} \\
\hline
\text{CON} & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 \\
\text{TRT1} & 2.50 & 3.50 & 3.98 & 4.03 & 3.68 & 3.35 & 3.02 & 2.94 \\
\text{TRT2} & 3.50 & 4.54 & 5.80 & 5.84 & 5.49 & 4.71 & 4.08 & 3.78 \\
\end{array}
$$

**Step 1: Creating the design**

**Create a data frame for the design**

```{r}
n_subject = 6 # Subjects per treatment
n_trt = 3 # Number of treatments
n_hour = 8 # Number of repeated measures (time points)
trt = c("CON", "TRT1", "TRT2")

df.rep <- data.frame(
  subject = as.factor(rep(seq_len(n_trt*n_subject), each = n_hour)),
  hour = as.factor(rep(seq_len(n_hour), n_subject*n_trt)),
  trt = rep(trt, each = n_subject*n_hour)
)
```

**Input templates**

Either values of `beta` or `means` are required in the following order:

```{r}
mkdesign(formula = ~ trt*hour, data = df.rep)
```

**Create the design**

```{r}
design.rep <- mkdesign(
formula = ~ trt*hour,
data = df.rep,
means =  c(1, 2.50, 3.5,
           1, 3.50, 4.54,
           1, 3.98, 5.80,
           1, 4.03, 5.4,
           1, 3.68, 5.49,
           1, 3.35, 4.71,
           1, 3.02, 4.08,
           1, 2.94, 3.78),
sigma2 = 2,
correlation = corAR1(value = 0.6, form = ~ hour|subject)
)
```

**Step 2: Power calculation**

Statistical power for main effects of treatment and time, and their interaction:

```{r}
pwr.anova(design.rep)
```

Statistical power for treatment differences at each time point:

```{r}
pwr.contrast(design.rep, which = "trt", by = "hour", contrast = "trt.vs.ctrl", p.adj = TRUE)[1:2]
```

As shown in *Example 5*, any design can be created using the `mkdesign` function as long as the appropriate data frame is provided. However, the package currently does not support non-normal response variables.

# Fundamental concepts

`pwr4exp` is developed based on LMM theory. The general form of an LMM can be expressed as:

$$
y = X\beta + Zu + \varepsilon
$$

where: $y$ represents the observations of the response variable, $\beta$ represents the fixed effect coefficients, $u$ denotes the random effects, where $u \sim N_q(0, G)$, $\varepsilon$ represents the random errors, where $\varepsilon \sim N_n(0, R)$, $X_{(n \times p)}$ and $Z_{(n \times q)}$ are the design matrices for the fixed and random effects, respectively.

It is assumed that $u$ and $\varepsilon$ are independent, and the marginal distribution of $y$ follows a normal distribution $y \sim N_n(X\beta, V)$, where:

$$
V = ZGZ^T + R
$$

## Inference on Treatment Effects

Inference on treatment effects often involves testing omnibus hypotheses and contrasts. These can be formulated using the general linear hypothesis:

$$
H_0: K\beta = 0
$$

where $K$ is a contrast matrix. When the variance-covariance parameters in $G$ and $R$ are known, the estimate of $\beta$ is:

$$
\hat{\beta} = (X^TV^{-1}X)^{-1}X^TV^{-1}y
$$

And its variance is:

$$
C = (X^TV^{-1}X)^{-1}
$$

The sampling distribution of $K'\hat{\beta}$ is:

$$
K'\hat{\beta} \sim N(0, K'CK)
$$

However, in practical situations, the matrices $G$ and $R$ are unknown and must be estimated using methods like Maximum Likelihood (ML) or Restricted ML (REML). The estimate of $\beta$ is obtained by plugging in the estimated covariance matrices $\hat{V}$, where:

$$
\hat{V} = Z\hat{G}Z^T + \hat{R}
$$

The resulting estimate of $\beta$ is:

$$
\hat{\beta} = (X^T\hat{V}^{-1}X)^{-1}X^T\hat{V}^{-1}y
$$

And its estimated variance is:

$$
\hat{C} = (X^T\hat{V}^{-1}X)^{-1}
$$

When testing the null hypothesis $H_0: K\beta = 0$, an approximate F-statistic is used. The F-statistic is given by:

$$
F = \frac{(K\hat{\beta})^T [K\hat{C}K^T]^{-1} (K\hat{\beta})}{v_1}
$$

$F$ follows an approximate F-distribution $F(v_1, v_2)$ under $H_0$, where $v_1 = \text{rank}(K) \geq 1$ represents the numerator degrees of freedom (df), $v_2$ is the denominator df.

When $\text{rank}(K) = 1$, the F-statistic simplifies to the square of the t-statistic:

$$
F = t^2
$$where $t = \frac{k'\hat{\beta}}{\sqrt{k'\hat{C}K}}$ with $v_2$ df.

In balanced designs, where data is analyzed using a variance components model, commonly applied in experimental animal research, $v_2$​ can be precisely determined through degrees of freedom decomposition, as applied in analysis of variance (ANOVA).

However, for more general cases, $v_2$ must be approximated using methods.

The Satterthwaite approximation (Satterthwaite, 1946) for DF in t-tests can be calculated as outlined by Giesbrecht and Burns (1985):

$$
v_2 = \frac{2(k^T \hat{C} k)^2}{g^T A g}
$$

where: $g$ is the gradient of $k^T C(\hat{\theta}) k$ with respect to $\theta$, the variance-covariance parameters in $V$, evaluated at $\hat{\theta}$. Matrix $A$ is the asymptotic variance-covariance matrix of $\hat{\theta}$, obtained from the information matrix of ML or REML estimation of $\hat{\theta}$ (Stroup, 2012).

In F-tests, $v_2$ can be calculated by following the procedures described by Fai and Cornelius (1996). First, $KC\hat{K}^T$ is decomposed to yield $KC\hat{K}^T = P^T D P$, where $P$ is an orthogonal matrix of eigen vectors, and $D$ is a diagonal matrix of eigenvalues. Define $k_m Ck_m^T$, where $k_m$ is the $m$-th row of $P K$, and let:

$$
v_m = \frac{2(D_m)^2}{g_m^T A g_m}
$$

where: $D_m$is the $m$th diagonal element of $D$, $g_m$ is the gradient of $k_m C k_m^T$ with respect to $\theta$, evaluated at $\hat{\theta}$.

Then let:

$$
E = \sum_{m=1}^{v_1} \frac{v_m}{v_m - 2} I(v_m > 2)
$$

where $I(v_m > 2)$ denotes the indicator function. The denominator DF $v_2$ is calculated as:

$$
v_2 = \frac{2E}{E - v_1}
$$ The Satterthwaite approximation can be applied in power analysis by plugging in the assumed values of variance parameters (Stroup, 2002).

## Power Calculation Under the Alternative Hypothesis

Under the alternative hypothesis $H_A: K'\beta \neq 0$, the F-statistic follows a non-central distribution $F(v_1, v_2, \phi)$, where $\phi$ is the non-centrality parameter that measures the departure from the null hypothesis $H_0$. The non-centrality parameter $\phi$ is given by:

$$
\phi = (K\hat{\beta})^T [K\hat{C}K^T]^{-1} (K\hat{\beta})
$$

Once the distribution of the F-statistic under $H_A$ is known, the power of the test can be calculated as the conditional probability of rejecting $H_0$ when $H_A$ is true:

$$
\text{Power} = P(\text{reject } H_0: F > F_{\text{crit}} \mid H_A)
$$

Where: $F_{\text{crit}}$ is the critical value of the F-statistic used to reject $H_0$, determined such that $P(F > F_{\text{crit}} \mid H_0) = \alpha$, where $\alpha$ is the type I error rate.

The determination of the degrees of freedom $v_1$ and $v_2$, as well as the non-centrality parameter $\phi$, are critical steps for power calculation.

Generally, power analysis requires specifying the following components:

-   The **Design** to be evaluated, which determines the matrices $X$ (fixed effects) and $Z$ (random effects).
-   The **Treatment Effects**, which determine $\beta$ (the fixed effect coefficients).
-   The **Variance Components**, which determine $G$ (the covariance matrix of the random effects) and $R$ (the covariance matrix of the residuals).

A key aspect of conducting a valid power analysis is obtaining **reasonable estimates for the magnitude of the parameters** that will be used in the model. This includes:

-   **Treatment Effects** ($\beta$): The size of the treatment effect(s) you expect to detect. This can be obtained from previous studies, pilot experiments, or subject-matter expertise.
-   **Variance Components** ($G$ and $R$): These represent the variability in the random effects and residual errors. Variance estimates are often obtained from:
    -   **Pilot studies** or preliminary data, which can provide initial estimates of variability in random effects (e.g., subject-to-subject variability or group-level variability).
    -   **Literature** on similar experiments, where variance components have been reported.
    -   **Subject-matter expertise**, where researchers provide estimates based on their knowledge of the system being studied.

Performing a power analysis with unrealistic parameter magnitudes can lead to incorrect conclusions, either overestimating the likelihood of detecting a treatment effect or requiring an unnecessarily large sample size.

# References

Fai, A. H., & Cornelius, P. L. (1996). Approximate F-tests of multiple degree of freedom hypotheses in generalized least squares analyses of unbalanced split-plot experiments. *Journal of Statistical Computation and Simulation, 54*(4), 363–378.

Giesbrecht, F. G., & Burns, J. C. (1985). Two-stage analysis based on a mixed model: Large-sample asymptotic theory and small-sample simulation results. *Biometrics, 41*(2), 477–486.

Satterthwaite, F. E. (1946). An approximate distribution of estimates of variance components. *Biometrics Bulletin, 2*(6), 110–114.

Stroup, W. W. (2002). Power analysis based on spatial effects mixed models: A tool for comparing design and analysis strategies in the presence of spatial variability. *Journal of Agricultural, Biological, and Environmental Statistics, 7*(4), 491–511.

Stroup, W. W. (2012). *Generalized linear mixed models: Modern concepts, methods, and applications* (1st ed.). CRC Press.