Abstract
The distinction between modeling a variable as a fixed versus a
random effect depends on the goal of the statistical analysis. While
some theory and software make a strong distinction,
variancePartition and dream take different
approaches based on the goal of each type of analysis. Here we consider
the distinction between fixed and random effects, and the usage of REML
in variancePartition and dream.
variancePartitionIn traditional statistics and biostatistics, there is a strong distinction between modeling categorical variants as fixed and random effects. Random effects correspond to a sample of units from a larger population, while fixed effects correspond to properties of specific individuals. Random effects are typically treated as nuisance variables and integrated out, and hypothesis testing is performed on the fixed effect.
The r2glmm package fits into this traditional framework,
by computing the variance fractions for a given fixed effect as: \[\begin{eqnarray}
\sigma^2_{fixed}/ \left(\sigma^2_{fixed} + \sigma^2_{error}\right)
\end{eqnarray}\]
Importantly, the random effects are not in the denominator. The fraction is only determined by fixed effects and residuals.
In my experience in bioinformatics, this was a problem. Making such distinctions between fixed and random effects seemed arbitrary. Variance in a phenotype could be due to age (fixed) or to variation across subject (random). Including all of the variables in the denominator produced more intuitive results so that 1) the variance fractions sum to one across all components and 2) fixed and random effects could be interpreted on the same scale 3) fractions could be compared across studies with different designs, 4) estimates of variance fractions were most accurate. So in variancePartition the fractions are defined as: \[\begin{eqnarray} \sigma^2_{X}/ \left(\sigma^2_{fixed} + \sigma^2_{random} + \sigma^2_{error}\right) \end{eqnarray}\]
just plugging the each variable in the numerator.
Thus the faction evaluated by variancePartition is different than
r2glmm by definition.
Here is some code explicitly demonstrating this difference:
library("variancePartition")
library("lme4")
library("r2glmm")
set.seed(1)
N <- 1000
beta <- 3
alpha <- c(1, 5, 7)
# generate 1 fixed variable and 1 random variable with 3 levels
data <- data.frame(X = rnorm(N), Subject = sample(c("A", "B", "C"), 100, replace = TRUE))
# simulate variable
# y = X\beta + Subject\alpha + \sigma^2
data$y <- data$X * beta + model.matrix(~ data$Subject) %*% alpha + rnorm(N, 0, 1)
# fit model
fit <- lmer(y ~ X + (1 | Subject), data, REML = FALSE)
# calculate variance fraction using variancePartition
# include the total sum in the denominator
frac <- calcVarPart(fit)
frac Subject X Residuals
0.4505 0.4952 0.0543 # the variance fraction excluding the random effect from the denominator
# is the same as from r2glmm
frac[["X"]] / (frac[["X"]] + frac[["Residuals"]])[1] 0.901 Effect Rsq upper.CL lower.CL
1 Model 0.896 0.904 0.886
2 X 0.896 0.904 0.886So the formulas are different. But why require categorical variables as random effects?
At practical level, categorical variables with too many levels are
problematic. Using a categorical variable with 200 categories as a fixed
effect is statistically unstable. There are so many degrees of freedom
that that variable will absorb a lot of variance even under the null.
Statistically, estimating the variance fraction for a variable with many
categories can be biased if that variable is a fixed effect. Therefore,
variancePartition requires all categorical variables to be
random effects. Modeling this variable as a random effect produces
unbiased estimates of variance fractions in practice. See simulations in
the Supplement (section 1.5) of Hoffman and Schadt
(2016).
The distinction between fixed and random effects is important in the formulation because it affects which variables are put in the denominator. So choosing to model a variable as a fixed versus random effect will definitely change the estimated fraction.
Yet for the variancePartition formulation, all variables
are in the denominator and it isn`t affected by the fixed/random
decision. Moreover, using a random effect empirically reduces the bias
of the estimated fraction.
Finally, why use maximum likelihood to estimate the paramters instead
of the default REML ()? Maximum likelihood fits all parameters jointly
so that it estimates the fixed and random effects together. This is
essential if we want to compare fixed and random effects later.
Conversely, REML estimates the random effects by removing the fixed
effects from the response before estimation. This implicitly removes the
fixed effects from the denominator when evaluating the variance
fraction. REML treats fixed effects as nuisance variables, while
variancePartition considers fixed effects to be a core part
of the analysis.
While REML produced unbiased estimates of the variance components,
the goal of variancePartition is to estimate the variance
fractions for fixed and random effects jointly. In simulations from the
Supplement (section 1.5) of Hoffman and Schadt
(2016), REML produced biased estimates of the variance fractions
while maximum likelihood estimates are unbiased.
dreamWhile dream is also based on a linear mixed model, the
goal of this analysis is to perform hypothesis testing on fixed effects.
Random effects are treated as nuisance variables to be integrated out,
and the approximate null distribution of a t- or F-statistic is
constructed from the model fit.
Since the goal of the analysis is different, the consideration of using REML versus ML is different than above. While is required by called by when , can be used with as either or . Since the Kenward-Roger method gave the best power with an accurate control of false positive rate in our simulations, and since the Satterthwaite method with gives p-values that are slightly closer to the Kenward-Roger p-values, is set as the default.
R Under development (unstable) (2024-10-21 r87258)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.1 LTS
Matrix products: default
BLAS: /home/biocbuild/bbs-3.21-bioc/R/lib/libRblas.so
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_GB LC_COLLATE=C
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: America/New_York
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] r2glmm_0.1.2 lme4_1.1-35.5 Matrix_1.7-1
[4] org.Hs.eg.db_3.20.0 msigdbr_7.5.1 GSEABase_1.69.0
[7] graph_1.85.0 annotate_1.85.0 XML_3.99-0.17
[10] AnnotationDbi_1.69.0 IRanges_2.41.0 S4Vectors_0.45.1
[13] Biobase_2.67.0 BiocGenerics_0.53.2 generics_0.1.3
[16] zenith_1.9.0 tximportData_1.35.0 tximport_1.35.0
[19] readr_2.1.5 edgeR_4.5.0 pander_0.6.5
[22] variancePartition_1.37.1 BiocParallel_1.41.0 limma_3.63.2
[25] ggplot2_3.5.1 knitr_1.49
loaded via a namespace (and not attached):
[1] jsonlite_1.8.9 magrittr_2.0.3
[3] farver_2.1.2 nloptr_2.1.1
[5] rmarkdown_2.29 zlibbioc_1.53.0
[7] vctrs_0.6.5 memoise_2.0.1
[9] minqa_1.2.8 RCurl_1.98-1.16
[11] progress_1.2.3 S4Arrays_1.7.1
[13] htmltools_0.5.8.1 curl_6.0.0
[15] broom_1.0.7 SparseArray_1.7.1
[17] sass_0.4.9 KernSmooth_2.23-24
[19] bslib_0.8.0 pbkrtest_0.5.3
[21] plyr_1.8.9 cachem_1.1.0
[23] lifecycle_1.0.4 iterators_1.0.14
[25] pkgconfig_2.0.3 R6_2.5.1
[27] fastmap_1.2.0 GenomeInfoDbData_1.2.13
[29] rbibutils_2.3 MatrixGenerics_1.19.0
[31] digest_0.6.37 numDeriv_2016.8-1.1
[33] colorspace_2.1-1 GenomicRanges_1.59.0
[35] RSQLite_2.3.7 filelock_1.0.3
[37] labeling_0.4.3 RcppZiggurat_0.1.6
[39] fansi_1.0.6 abind_1.4-8
[41] httr_1.4.7 compiler_4.5.0
[43] bit64_4.5.2 aod_1.3.3
[45] withr_3.0.2 backports_1.5.0
[47] DBI_1.2.3 gplots_3.2.0
[49] MASS_7.3-61 DelayedArray_0.33.1
[51] corpcor_1.6.10 gtools_3.9.5
[53] caTools_1.18.3 tools_4.5.0
[55] remaCor_0.0.18 glue_1.8.0
[57] nlme_3.1-166 grid_4.5.0
[59] reshape2_1.4.4 snow_0.4-4
[61] gtable_0.3.6 tzdb_0.4.0
[63] tidyr_1.3.1 hms_1.1.3
[65] utf8_1.2.4 XVector_0.47.0
[67] pillar_1.9.0 stringr_1.5.1
[69] babelgene_22.9 vroom_1.6.5
[71] splines_4.5.0 dplyr_1.1.4
[73] BiocFileCache_2.15.0 lattice_0.22-6
[75] bit_4.5.0 tidyselect_1.2.1
[77] locfit_1.5-9.10 Biostrings_2.75.1
[79] SummarizedExperiment_1.37.0 RhpcBLASctl_0.23-42
[81] xfun_0.49 statmod_1.5.0
[83] matrixStats_1.4.1 KEGGgraph_1.67.0
[85] stringi_1.8.4 UCSC.utils_1.3.0
[87] yaml_2.3.10 boot_1.3-31
[89] evaluate_1.0.1 codetools_0.2-20
[91] archive_1.1.10 tibble_3.2.1
[93] Rgraphviz_2.51.0 cli_3.6.3
[95] RcppParallel_5.1.9 xtable_1.8-4
[97] Rdpack_2.6.1 munsell_0.5.1
[99] jquerylib_0.1.4 Rcpp_1.0.13-1
[101] GenomeInfoDb_1.43.0 EnvStats_3.0.0
[103] dbplyr_2.5.0 png_0.1-8
[105] Rfast_2.1.0 parallel_4.5.0
[107] blob_1.2.4 prettyunits_1.2.0
[109] bitops_1.0-9 mvtnorm_1.3-2
[111] lmerTest_3.1-3 scales_1.3.0
[113] purrr_1.0.2 crayon_1.5.3
[115] fANCOVA_0.6-1 rlang_1.1.4
[117] EnrichmentBrowser_2.37.0 KEGGREST_1.47.0