---
title: "Sample Size Planning"
output:
  html_document:
    df_print: paged
    toc: true
---

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

## Load Packages

```{r load-pkg, message=FALSE}
# To install a package, run the following ONCE (and only once on your computer)
# install.packages("simr")  
library(tidyverse)
library(lme4)
library(simr)
library(MuMIn)
theme_set(theme_bw())  # Theme; just my personal preference
```

See also an example from the vignette of the `simr` package: https://cran.r-project.org/web/packages/simr/vignettes/fromscratch.html

Note: In the examples of this note I used 100 simulation samples only to save
time. In practice, you should use at least 1,000 to obtain a reasonable 
precision on the power estimation. 

## Binary Predictor at Level 2

$$
  \begin{aligned}
    Y_{ij} & = \beta_{0j} + e_{ij} \\
    \beta_{0j} & = \gamma_{00} + {\color{red}{\gamma_{01}}} T_j + u_{0j}
  \end{aligned}
$$

### Preparing simulated data

First, generate the predictor matrix and the cluster ID

```{r sim_dat}
# Cluster ID (start with 10 clusters)
num_clus <- 10
id <- 1:num_clus
# Binary lv-2 predictor
treat <- factor(c("c", "t"))
# Level-2 data
lv2_dat <- tibble(id, treat = rep_len(treat, length.out = num_clus))
# Cluster size
clus_size <- 5
# Expand each cluster to include more rows
(sim_dat <- 
    lv2_dat %>% 
    slice(rep(1:n(), each = clus_size)))
```

Then specify the fixed effect coefficients and the random effect variance:

```{r m1-params}
gams <- c(0, 0.625)  # gamma_00, gamma_01
taus <- matrix(0.25)  # conditional ICC = 0.25 / (1 + 0.25) = 0.2
sigma <- 1  # sigma, usually fixed to be 1.0 for easier scaling
```

Use `simr::makeLmer()` to build an artificial `lmer` object:

```{r sim_m1}
sim_m1 <- makeLmer(y ~ treat + (1 | id), 
                   fixef = gams, VarCorr = taus, sigma = sigma, 
                   data = sim_dat)
# Check R^2
r.squaredGLMM(sim_m1)
# Cohen's d
gams[2] / sqrt(taus + sigma^2)
```

### Obtain simulation-based power

Now, check the power:

```{r power_m1, message=FALSE}
power_m1 <- powerSim(sim_m1, 
                     # number of simulated samples; use a small number for 
                     # trial run
                     nsim = 20,
                     progress = FALSE)  # suppress the progress bar
power_m1
```

Power was obviously not sufficient as it is far from 80%, so we'll increase it. 

#### Increasing number of clusters

Increase to 50 clusters

```{r sim_m1_50J, message=FALSE, cache=TRUE}
sim_m1_50J <- extend(sim_m1, 
                     along = "id",  # add more levels of `id` (i.e., clusters)
                     n = 50)  # increase to 50 clusters
power_m1_50J <- powerSim(sim_m1_50J, nsim = 20, progress = FALSE)
power_m1_50J
```

This looks close to 80%. Let's try several numbers of clusters (30, 35, 40, 45, 50) with a larger `nsim`. To my knowledge `simr` does not directly handle multiple values for number of clusters, so I provided you a handy `power_extend()` function, and an associated function for plotting the results.  

```{r power_extend}
# Wrap to a function
power_extend <- function(fit, along, within, new_size, time_points, ...) {
  if (!missing(time_points)) {
    power_lst <- vector("list", length(time_points))
    for (j in seq_along(time_points)) {
      sim_m <- extend(fit, along = along, within = within, 
                      values = (1:time_points[j]) - 1)
      power_lst[[j]] <- powerSim(sim_m, ...)
    }
    setNames(power_lst, time_points)
  } else {
    power_lst <- vector("list", length(new_size))
    for (j in seq_along(new_size)) {
      sim_m <- extend(fit, along = along, within = within, n = new_size[j])
      power_lst[[j]] <- powerSim(sim_m, ...)
    }
    setNames(power_lst, new_size)
  }
}
# Compute proportions and CI
powerinterval <- function(object, alpha = object$alpha, level = 0.95, 
                          method = getSimrOption("binom"), ...) {
    x <- sum(object$pval < alpha, na.rm = TRUE)
    n <- object$n
    if (is.na(x) || is.na(n) || (n == 0)) {
        return(rep(NA, 3))
    }
    binom::binom.confint(x, n, level, method, ...)[c("mean", "lower", "upper")]
}
# Add a plotting function to get power curve
plot_power_extend <- function(object, target = .80, ...) {
  nsim <- object[[1]]$n
  df <- purrr::map_dfr(object, powerinterval, .id = "n", ...)
  df$n <- as.integer(df$n)
  ggplot(df, aes(x = n, y = mean, group = 1)) + 
    geom_hline(yintercept = target) + 
    geom_smooth(method = "glm", 
                formula = cbind(y * nsim, (1 - y) * nsim) ~ x, 
                method.args = list(family = binomial("probit"))) + 
    geom_pointrange(aes(ymin = lower, ymax = upper)) + 
    labs(y = "power") + 
    ylim(0, 1)
}
```

```{r powers_m1, cache=TRUE, message=FALSE}
powers_m1 <- power_extend(sim_m1, along = "id", 
                          new_size = seq.int(30, to = 50, by = 5), 
                          nsim = 100)
# Plot
plot_power_extend(powers_m1)
```

This suggests that about 40 clusters are needed, if each cluster contains 5 observations. 

#### Increasing cluster size

Instead of increasing the number of cluster, one can also increase the cluster
size:

```{r power_m1_25n, message=FALSE, cache=TRUE}
sim_m1_25n <- extend(sim_m1, 
                     # need to specify the combination of id and lv-2 predictors
                     within = "id + treat",
                     n = 25)  # cluster size
power_m1_25n <- powerSim(sim_m1_25n, nsim = 100, progress = FALSE)
power_m1_25n
```

Note that the power of increasing cluster size by 5 times is much smaller than
that of increasing number of clusters by 5 times. Here is the power curve:

```{r powers_m1-2, cache=TRUE, message=FALSE}
# Use the previously defined function
powers_m1 <- power_extend(sim_m1, within = "id + treat", 
                          new_size = c(5, 10, 25, 50), 
                          nsim = 100)
plot_power_extend(powers_m1)
```

In this case, increasing the cluster size only increases the power to close to 50%. There is an upper limit in power improvement with increasing the cluster size, when the number of cluster is kept constant. 

## Continuous Predictor at Level 2

$$
  \begin{aligned}
    Y_{ij} & = \beta_{0j} + e_{ij} \\
    \beta_{0j} & = \gamma_{00} + {\color{red}{\gamma_{01}}} W_j + u_{0j}
  \end{aligned}
$$

### Preparing Simulated Data

First, generate the predictor matrix and the cluster ID

```{r sim_dat-m2}
# Cluster ID (start with 10 clusters)
num_clus <- 10
id <- 1:num_clus
# Binary lv-2 predictor
w <- rnorm(num_clus)
# Force SD = 1
w <- (w - mean(w)) / sd(w) * (num_clus - 1) / num_clus
# Level-2 data
lv2_dat <- tibble(id, w = w)
# Cluster size
clus_size <- 5
# Expand each cluster to include more rows
(sim_dat <- 
    lv2_dat %>% 
    slice(rep(1:n(), each = clus_size)))
```

Then specify the fixed effect coefficients and the random effect variance:

```{r m2-params}
gams <- c(1, 0.3)  # gamma_00, gamma_01
taus <- matrix(0.5)  # conditional ICC = 0.5 / (1 + 0.5) = 0.33
sigma <- 1  # sigma, usually fixed to be 1.0 for easier scaling
```

Use `simr::makeLmer()` to build an artificial `lmer` object:

```{r sim_m2}
sim_m2 <- makeLmer(y ~ w + (1 | id), 
                   fixef = gams, VarCorr = taus, sigma = sigma, 
                   data = sim_dat)
# Check R^2
r.squaredGLMM(sim_m2)
```

### Obtain simulation-based power

Use various numbers of clusters

```{r powers_m2, cache=TRUE, message=FALSE}
# Use it on our example:
powers_m2 <- power_extend(sim_m2, along = "id", 
                          new_size = c(10, 20, 50, 100), 
                          nsim = 100)
plot_power_extend(powers_m2)
```

Use various cluster sizes

```{r powers_m2-2, cache=TRUE, message=FALSE}
# Use it on our example:
powers_m2 <- power_extend(sim_m2, 
                          within = "id + w", 
                          new_size = c(5, 10, 25, 50), 
                          nsim = 100)
plot_power_extend(powers_m2)
```

## Binary Predictor at Level 2 With a Continuous Level-1 Covariate (and Random Slope)

Cluster-mean centering needs to be considerd

$$
  \begin{aligned}
    Y_{ij} & = \beta_{0j} + \beta_{1j} X\text{_cmc}_{ij} + e_{ij} \\
    \beta_{0j} & = \gamma_{00} + \gamma_{01} W_j + \gamma_{02} X\text{_cm}_{j} + u_{0j} \\
    \beta_{1j} & = \gamma_{10} + {\color{red}{\gamma_{11}}} W_j + u_{1j}
  \end{aligned}
$$  

### Preparing Simulated Data

First, generate the predictor matrix and the cluster ID

```{r sim_dat-m4}
# Cluster ID (start with 10 clusters)
num_clus <- 10
id <- 1:num_clus
# Binary lv-2 predictor
treat <- factor(c("c", "t"))
# Cluster means of X, with ICC(X) = 0.2
x_cm <- rnorm(num_clus, mean = 0, sd = sqrt(0.2))
# Force SD = sqrt(0.2)
x_cm <- (x_cm - mean(x_cm)) / sd(x_cm) * sqrt(0.2) * (num_clus - 1) / num_clus
# Level-2 data
lv2_dat <- tibble(id, treat = rep_len(treat, length.out = num_clus), 
                  x_cm)
# Expand each cluster to include more rows
clus_size <- 20  # Cluster size
lv2_dat <- lv2_dat %>% 
    slice(rep(1:n(), each = clus_size))
# Within-cluster component of X, with sigma^2(X) = 0.8
num_obs <- num_clus * clus_size
x_cmc <- rnorm(num_obs, mean = 0, sd = 1)
x_cmc <- x_cmc - ave(x_cmc, lv2_dat$id, FUN = mean)
# Force SD = sqrt(0.8)
x_cmc <- x_cmc / sd(x_cmc) * sqrt(0.8) * (num_obs - 1) / num_obs
# Expand each cluster to include more rows
(sim_dat <- mutate(lv2_dat, x_cmc = x_cmc))
```

Then specify the fixed effect coefficients and the random effect variance:

```{r m4-params}
# gamma_00, gamma_01, gamma_02, gamma_10, gamma_11
gams <- c(0, 0.3, 0.2, 0.1, 0.2)  
taus <- matrix(c(0.25, 0, 
                 0, 0.10), nrow = 2)  # tau^2_0 = 0.25, tau^2_1 = 0.1
sigma <- sqrt(.96)  # sigma
```

Use `simr::makeLmer()` to build an artificial `lmer` object:

```{r sim_m4}
sim_m4 <- makeLmer(y ~ treat + x_cm + x_cmc + treat:x_cmc + (x_cmc | id), 
                   fixef = gams, VarCorr = taus, sigma = sigma, 
                   data = sim_dat)
# Check R^2
r.squaredGLMM(sim_m4)
# Compare to a model without the cross-level interaction
sim_m4n <- makeLmer(y ~ treat + x_cm + x_cmc + (x_cmc | id), 
                    fixef = gams[1:4], VarCorr = taus, sigma = sigma, 
                    data = sim_dat)
# Check R^2
r.squaredGLMM(sim_m4n)
# So the cross-level interaction corresponds to an R^2 change of about 2.5%
```

### Obtain simulation-based power

Use various numbers of clusters

Note: by default the first term in the fixed effect (other than the intercept)
will be tested. We'll specify the interaction instead.

```{r powers_m4, cache=TRUE, message=FALSE}
# Use it on our example:
powers_m4 <- power_extend(sim_m4, 
                          along = "id", 
                          new_size = c(10, 20, 50, 100), 
                          nsim = 100, 
                          test = fixed("treat:x_cmc", 
                                       method = "anova"))
plot_power_extend(powers_m4)
```

Use various cluster sizes (results not shown)

```{r powers_m4-2, eval=FALSE}
# Use it on our example:
powers_m4 <- power_extend(sim_m4, 
                          within = "id + treat + x_cm",
                          new_size = c(40, 80, 200), 
                          nsim = 100, 
                          test = fixed("treat:x_cmc", 
                                       method = "anova"))
plot_power_extend(powers_m4)
```

## Growth Modeling

With a binary person-level predictor

$$
  \begin{aligned}
    Y_{ti} & = \beta_{0j} + \beta_{1i} \text{time}_{ti} + e_{ti} \\
    \beta_{0i} & = \gamma_{00} + {\color{red}{\gamma_{01}}} W_i + u_{0i} \\
    \beta_{1i} & = \gamma_{10} + \gamma_{11} W_i + u_{1i}
  \end{aligned}
$$

### Preparing Simulated Data

First, generate the predictor matrix and the cluster ID

```{r sim_dat-m5}
# Person ID (start with 10 clusters)
num_clus <- 10
id <- 1:num_clus
# Binary lv-2 predictor
treat <- factor(c("c", "t"))
# Level-2 data
lv2_dat <- tibble(id, treat = rep_len(treat, length.out = num_clus))
# Expand each cluster to include more rows
clus_size <- 4  # Number of time points
# Within-cluster component of X, with sigma^2(X) = 0.8
time <- 0:(clus_size - 1)
# Expand each cluster to include more rows
(sim_dat <- expand_grid(lv2_dat, time))
```

Then specify the fixed effect coefficients and the random effect variance:

```{r m5-params}
# gamma_00, gamma_01, gamma_10, gamma_11
gams <- c(1, 0.5, -0.2, 0.3)  
taus <- matrix(c(0.5, 0, 
                 0, 0.10), nrow = 2)  # tau^2_0 = 0.5, tau^2_1 = 0.1
sigma <- sqrt(1)  # sigma
```

Use `simr::makeLmer()` to build an artificial `lmer` object:

```{r sim_m5}
sim_m5 <- makeLmer(y ~ treat + time + treat:time + (time | id), 
                   fixef = gams, VarCorr = taus, sigma = sigma, 
                   data = sim_dat)
# Check R^2
r.squaredGLMM(sim_m5)
# Compare to a model without the binary treatment predictor
sim_m5n <- makeLmer(y ~ time + (time | id), 
                    fixef = gams[c(1, 3)], VarCorr = taus, sigma = sigma, 
                    data = sim_dat)
# Check R^2
r.squaredGLMM(sim_m5n)
# So treat corresponds to an R^2 change of about 10%
```

### Obtain simulation-based power

Use various numbers of clusters

Note: by default the first term in the fixed effect (other than the intercept)
will be tested. We'll first specify the effect of treat.

```{r powers_m5, cache=TRUE, message=FALSE}
# Use it on our example:
powers_m5 <- power_extend(sim_m5, 
                          along = "id", 
                          new_size = c(10, 20, 50, 100), 
                          nsim = 100)
plot_power_extend(powers_m5)
```

And then the slope difference

```{r powers_m5-2, cache=TRUE, message=FALSE}
# Use it on our example:
powers_m5 <- power_extend(sim_m5, 
                          along = "id", 
                          new_size = c(10, 20, 50, 100), 
                          nsim = 100, 
                          test = fixed("treat:time", 
                                       method = "anova"))
plot_power_extend(powers_m5)
```

Use more time points

```{r powers_m5-3, cache=TRUE, message=FALSE}
# Use the `time_points` argument
powers_m5 <- power_extend(sim_m5, 
                          along = "time", 
                          time_points = c(4, 8), 
                          nsim = 100)
plot_power_extend(powers_m5)
```

Increasing number of time points does not help the power with the intercept difference. Now let's look at the slope difference. 

```{r powers_m5-4, cache=TRUE, message=FALSE}
# Use the `time_points` argument
powers_m5 <- power_extend(sim_m5, 
                          along = "time", 
                          time_points = c(4, 8), 
                          nsim = 100, 
                          test = fixed("treat:time", 
                                       method = "anova"))
plot_power_extend(powers_m5)
```
