Sample size calculation for linear mixed model in R

Question

I am trying to estimate appropriate sample size for an analysis using a linear mixed model.

In this analysis, I have my outcome, pain intensity ranging from 0-10, regressed against 4 covariates that are sedentary time, standing time, physical activity time, and other behaviors time.

For this calculation, I am using a cohort that has pain measurements from 14 times over one year, and the covariates are measured only at baseline. I would like to investigate how baseline covariates are associated with future development of pain.

The standard deviation of pain intensity at baseline is 2.94. A clinical effect size we expect to see is '1' on a scale of 0-10.

Also, the intra-class correlation between measurements of pain is 0.05.

Can somebody tell me how I can perform the sample size calculation in R software, at a statistical power of 80% and a p value of .05?

Answer 1

Assume the model

$$ Y_{ij} = \beta_0 + \beta X_i + \mu_i + \varepsilon_{ij} $$

where $X$ is an indicator variable for the group (A or B, say), $\mu_i \sim \mathcal{N}(0 , \sigma_\mu^2)$ and $\varepsilon_{ij} \sim \mathcal{N}(0,\sigma_e^2)$. We have $\text{var}(Y_{ij})= \sigma_\mu^2 + \sigma_e^2 = \sigma^2$ and the intra-class correlation is given by \begin{align*} \rho &= \frac{\text{Cov}(Y_{ij}, Y_{ij'})}{\text{Var}(Y_{ij})} \\ &= \frac{\sigma_\mu^2}{\sigma^2} \\ \end{align*}

Let the hypotheses $H_0$ and $H_1$ be:

$$ H_0 : \beta =0 \quad H_1 : \beta \neq 0 $$

Now say that both groups are composed of $n$ subjects each of them having $k$ measures and let

$$ \bar{X}_{*} = \frac{1}{n} \sum_{i=1}^n \frac{1}{k} \sum_{j=1}^k Y_{*ij} $$

for $* \in \{A,B \}$ be the sample mean for a group.

A standard test statistic can be written as

$$ Z = \frac{\bar{X}_A - \bar{X}_B}{\sqrt{\text{var}(\bar{X}_A - \bar{X}_B)}} $$

In the case of repeated measurement the sample size required should be adapted because $\text{var}(\bar{X}_A - \bar{X}_B)$ is not as a variance for iid observations, indeed :

\begin{align*} \text{var}(\bar{X}_A - \bar{X}_B) &= \text{var}(\bar{X}_A ) + \text{var}(\bar{X}_B) \\ &= 2 \text{var}(\bar{X}_A ) \\ &= 2 \text{var} \Big ( \frac{1}{n} \sum_{i=1}^n \frac{1}{k} \sum_{j=1}^k Y_{Aij} \Big ) \\ &= \frac{2}{nk^2} \text{var} \Big ( \sum^k ( \mu_i + \varepsilon{ij} ) \Big ) \\ &= \frac{2}{nk^2} \text{var} \Big ( k \mu_i + \sum^k \varepsilon{ij} \Big ) \\ &= \frac{2}{nk^2} \Big ( k^2 \sigma_\mu^2 + k \sigma_e^2 \Big )^* \\ & = \frac{2}{n} \Big ( \sigma_\mu^2 + \frac{\sigma_e^2}{k} \Big ) \\ \end{align*}

$^*$ This line is because $\mu_i$ and $\varepsilon_{ij}$ are independent.

We can now compute a sample size but with the adjusted variance and with $k$ known (14 in your case).

For iid observations with variance $\sigma^2_*$ the sample size required per group (assuming balanced groups, $\alpha$ and $\beta$ as type I and II errors, $\Delta$ as a mean difference between both groups and a two-sided test) is:

$$ n= \frac{\big( Z_{1-\alpha/2} + Z_{1-\beta} )^2}{\Delta^2} \times 2\sigma^2_* $$

where $Z_q$ is the $q^{\text{th}}$ quantile of a standard normal distribution.

We can use this formula accounting for the corrected variance, hence the number of subjects per arm is given by :

\begin{align*} n &= \frac{\big( Z_{1-\alpha/2} + Z_{1-\beta} )^2}{\Delta^2} \times 2 \Big ( \sigma_\mu^2 + \frac{\sigma_e^2}{k} \Big ) \\ &= \frac{\big( Z_{1-\alpha/2} + Z_{1-\beta} )^2}{\Delta^2} \times \frac{2 (\sigma_\mu^2 + \sigma_e^2)} {k} \Big ( 1 + (k-1)\rho \Big ) \end{align*}