Required number of levels for a random effects

statistics

mixed-models

Author

Thierry Onkelinx

Published

September 3, 2018

When we analyse a mixed models, the question often arises whether a covariate should be used as a random effect or as a fixed effects. Let’s assume a simple design. Two types of fertilizer are tested on a number of fields (\(n_s\)). Each field is split in two and the fertilizers are assigned at random to these halves while making sure that each field has both treatments. From a conceptual point of view, we are only interested in the effect of the fertilizers. However, we need to take the potential effect of the field into account. Hence we use fertiliser as a fixed effect and field as a random effect.

Example design
field	Part 1	Part2
1	A	B
2	A	B
3	A	B
4	B	A
5	B	A
6	B	A
7	B	A
8	A	B
9	A	B
10	B	A

While this makes conceptually sense, we might run into computational problems. Instead of estimating each individual field effect, we want to estimate the variability due to field effect. This is the variance of the random effect, rather than the estimates for the individual random effect levels. But how precise is the estimate of this random effect variance? The sample variance \(s^2\) follows a scaled \(\chi^2\) distribution.

\[(n_s - 1)\frac{s^2}{\sigma^2} \sim \chi^2_{n_s-1}\]

This equation makes it straightforward to calculate the distribution of the \(r = \sigma^2/s^2\) ratio for a given \(n_s\). We can rewrite this as \(\sigma^2=rs^2\) or the true variance is a multiple \(r\) of the estimated variance \(s^2\). The distribution of \(r\) is derived below.

\[(n_s - 1)\frac{s^2}{\sigma^2} \sim \chi^2_{n_s-1}\]

\[\frac{s^2}{\sigma^2} \sim \frac{\chi^2_{n_s-1}}{(n_s - 1)}\]

\[r = \frac{\sigma^2}{s^2} \sim \frac{(n_s - 1)}{\chi^2_{n_s-1}}\]

Even when with a large number of levels (\(n_s = 1000\)), there still is a reasonable amount of uncertainty around \(\sigma^2\) since the 95% interval of \(r\) ranges from 0.918 to 1.094. The full distribution for \(r\) when \(n_s = 1000\) is shown in the figure below.

Figure 1: Distribution of the relative true variance (\(\frac{\sigma^2}{s^2}\)) when \(n_s = 1000\). Dashed line: \(\frac{\sigma^2}{s^2} = 1\) Dotted lines: 2.5% and 97.5% quantiles of \(\frac{\sigma^2}{s^2}\). Dash-dotted line: median of \(\frac{\sigma^2}{s^2}\).

A thousand random effect levels is not always feasible. So what will happen if we use a more realistic number of random effect levels, e.g. \(n_s = 20\). The distribution of \(r\) will be wider as \(n_s\) decreases. At \(n_s = 20\), the 95% interval of \(r\) will range from 0.578 to 2.133.

Figure 2: Distribution of the relative true variance (\(\frac{\sigma^2}{s^2}\)) when \(n_s = 20\). Dashed line: \(\frac{\sigma^2}{s^2} = 1\) Dotted lines: 2.5% and 97.5% quantiles of \(\frac{\sigma^2}{s^2}\). Dash-dotted line: median of \(\frac{\sigma^2}{s^2}\).

How low can we go? The figure below depicts the density of \(r\) when \(n_s = 4\). The 95% interval of \(r\) ranges from 0.321 to 13.902, so we can hardly do any inference on the estimated variance.

Figure 3: Distribution of the relative true variance (\(\frac{\sigma^2}{s^2}\)) when \(n_s = 4\). Dashed line: \(\frac{\sigma^2}{s^2} = 1\) Dotted lines: 2.5% and 97.5% quantiles of \(\frac{\sigma^2}{s^2}\). Dash-dotted line: median of \(\frac{\sigma^2}{s^2}\).

Recommendations

We see that the number of random effect levels has a strong impact on the uncertainty of the estimate variance. The figure below displays the 97.5% quantile of \(r\) for a large range of \(n_s\) values. Since the right tail of \(r\) is heavier than the left tail, it is a good idea to focus on the 97.5% quantile. Based on this figure we can give the following recommendations:

get \(n_s > 1000\) levels when an accurate estimate of the random effect variance is crucial. E.g. when a single number will be use for power calculations.
get \(n_s > 100\) levels when a reasonable estimate of the random effect variance is sufficient. E.g. power calculations with sensitivity analysis of the random effect variance.
get \(n_s > 20\) levels for an experimental study
in case \(10 < n_s < 20\) you should validate the model very cautious before using the output
in case \(n_s < 10\) it is safer to use the variable as a fixed effect.

Figure 4: Uncertainty on the random effect variance in function of the number of random effect levels.

Session info

These R packages were used to create this post.

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.1 (2023-06-16)
 os       Ubuntu 22.04.3 LTS
 system   x86_64, linux-gnu
 ui       X11
 language nl_BE:nl
 collate  nl_BE.UTF-8
 ctype    nl_BE.UTF-8
 tz       Europe/Brussels
 date     2023-09-02
 pandoc   3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package     * version    date (UTC) lib source
 cli           3.6.1      2023-03-23 [1] CRAN (R 4.3.0)
 colorspace    2.1-0      2023-01-23 [1] CRAN (R 4.3.0)
 digest        0.6.32     2023-06-26 [1] CRAN (R 4.3.1)
 dplyr       * 1.1.2      2023-04-20 [1] CRAN (R 4.3.0)
 evaluate      0.21       2023-05-05 [1] CRAN (R 4.3.0)
 fansi         1.0.4      2023-01-22 [1] CRAN (R 4.3.0)
 farver        2.1.1      2022-07-06 [1] CRAN (R 4.3.0)
 fastmap       1.1.1      2023-02-24 [1] CRAN (R 4.3.0)
 forcats     * 1.0.0      2023-01-29 [1] CRAN (R 4.3.0)
 generics      0.1.3      2022-07-05 [1] CRAN (R 4.3.0)
 ggplot2     * 3.4.2      2023-04-03 [1] CRAN (R 4.3.0)
 glue          1.6.2      2022-02-24 [1] CRAN (R 4.3.0)
 gtable        0.3.3      2023-03-21 [1] CRAN (R 4.3.0)
 hms           1.1.3      2023-03-21 [1] CRAN (R 4.3.0)
 htmltools     0.5.5      2023-03-23 [1] CRAN (R 4.3.0)
 htmlwidgets   1.6.2      2023-03-17 [1] CRAN (R 4.3.0)
 jsonlite      1.8.7      2023-06-29 [1] CRAN (R 4.3.1)
 knitr       * 1.43       2023-05-25 [1] CRAN (R 4.3.0)
 labeling      0.4.2      2020-10-20 [1] CRAN (R 4.3.0)
 lifecycle     1.0.3      2022-10-07 [1] CRAN (R 4.3.0)
 lubridate   * 1.9.2.9000 2023-05-15 [1] https://inbo.r-universe.dev (R 4.3.0)
 magrittr      2.0.3      2022-03-30 [1] CRAN (R 4.3.0)
 munsell       0.5.0      2018-06-12 [1] CRAN (R 4.3.0)
 pillar        1.9.0      2023-03-22 [1] CRAN (R 4.3.0)
 pkgconfig     2.0.3      2019-09-22 [1] CRAN (R 4.3.0)
 purrr       * 1.0.1      2023-01-10 [1] CRAN (R 4.3.0)
 R6            2.5.1      2021-08-19 [1] CRAN (R 4.3.0)
 readr       * 2.1.4      2023-02-10 [1] CRAN (R 4.3.0)
 rlang         1.1.1      2023-04-28 [1] CRAN (R 4.3.0)
 rmarkdown     2.23       2023-07-01 [1] CRAN (R 4.3.1)
 rstudioapi    0.14       2022-08-22 [1] CRAN (R 4.3.0)
 scales        1.2.1      2022-08-20 [1] CRAN (R 4.3.0)
 sessioninfo   1.2.2      2021-12-06 [1] CRAN (R 4.3.0)
 stringi       1.7.12     2023-01-11 [1] CRAN (R 4.3.0)
 stringr     * 1.5.0      2022-12-02 [1] CRAN (R 4.3.0)
 tibble      * 3.2.1      2023-03-20 [1] CRAN (R 4.3.0)
 tidyr       * 1.3.0      2023-01-24 [1] CRAN (R 4.3.0)
 tidyselect    1.2.0      2022-10-10 [1] CRAN (R 4.3.0)
 tidyverse   * 2.0.0      2023-02-22 [1] CRAN (R 4.3.0)
 timechange    0.2.0      2023-01-11 [1] CRAN (R 4.3.0)
 tzdb          0.4.0      2023-05-12 [1] CRAN (R 4.3.0)
 utf8          1.2.3      2023-01-31 [1] CRAN (R 4.3.0)
 vctrs         0.6.3      2023-06-14 [1] CRAN (R 4.3.0)
 withr         2.5.0      2022-03-03 [1] CRAN (R 4.3.0)
 xfun          0.39       2023-04-20 [1] CRAN (R 4.3.0)
 yaml          2.3.7      2023-01-23 [1] CRAN (R 4.3.0)

 [1] /home/thierry/R/x86_64-pc-linux-gnu-library/4.3
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────