Sample size calculation for linear mixed model in R

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?

• statistical power of 80% and a p value of .05? please explain briefly – Subhash C. Davar Jul 6 at 11:09

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*}

• @usεr11852 $\Delta$ is a mean difference between both groups under an alternative hypothesis. I edited my answer. – winperikle Jul 5 at 16:00
• Cool. Nice answer. +1 – usεr11852 says Reinstate Monic Jul 5 at 18:56
• thanks a lot for your answer. although i understand roughly your answer, it is kind of difficult for me to undertstand all these annotations. is there a way by which i can translate this in r? also, i see my exposure variables are continous in nature, so i dont have two groups A and B. i this case, how do i translate this in my case? – nidhi gupta Jul 8 at 13:14
• @nidhigupta R is not a requirement, a basic calculator can carry this out, but you can use R. In order to do so, you need to replace the different variables in the formula by the values corresponding to what you want (like $Z_{1-\alpha/2}$ by $1.96$). If you only have one group and want to compare its mean with a theoretical one, just replace $\bar{X}_B$ by this constant and see that the variance of the sample mean minus this constant is half the variance of the difference between two groups and then adapt the sample size. – winperikle Jul 9 at 9:48
• perfect!! i understand now. thanks once again! – nidhi gupta Jul 9 at 13:03