variancePartition

Estimating contributions to expression variation

In 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
# using r2glmm
r2beta(fit)
  Effect   Rsq upper.CL lower.CL
1  Model 0.896    0.904    0.886
2      X 0.896    0.904    0.886

So 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.

dream

Hypothesis testing

While 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.

Session Info

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.1 LTS

Matrix products: default
BLAS:   /media/volume/teran2_disk/biocbuild/bbs-3.20-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.68.0         
 [7] graph_1.84.0             annotate_1.84.0          XML_3.99-0.17           
[10] AnnotationDbi_1.68.0     IRanges_2.40.0           S4Vectors_0.44.0        
[13] Biobase_2.66.0           BiocGenerics_0.52.0      zenith_1.8.0            
[16] tximportData_1.34.0      tximport_1.34.0          readr_2.1.5             
[19] edgeR_4.4.0              pander_0.6.5             variancePartition_1.36.2
[22] BiocParallel_1.40.0      limma_3.62.1             ggplot2_3.5.1           
[25] 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.52.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.6.0             
 [13] htmltools_0.5.8.1           curl_6.0.0                 
 [15] broom_1.0.7                 SparseArray_1.6.0          
 [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.18.0      
 [31] digest_0.6.37               numDeriv_2016.8-1.1        
 [33] colorspace_2.1-1            GenomicRanges_1.58.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.4.1             
 [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.32.0        
 [51] corpcor_1.6.10              gtools_3.9.5               
 [53] caTools_1.18.3              tools_4.4.1                
 [55] remaCor_0.0.18              glue_1.8.0                 
 [57] nlme_3.1-166                grid_4.4.1                 
 [59] reshape2_1.4.4              generics_0.1.3             
 [61] snow_0.4-4                  gtable_0.3.6               
 [63] tzdb_0.4.0                  tidyr_1.3.1                
 [65] hms_1.1.3                   utf8_1.2.4                 
 [67] XVector_0.46.0              pillar_1.9.0               
 [69] stringr_1.5.1               babelgene_22.9             
 [71] vroom_1.6.5                 splines_4.4.1              
 [73] dplyr_1.1.4                 BiocFileCache_2.14.0       
 [75] lattice_0.22-6              bit_4.5.0                  
 [77] tidyselect_1.2.1            locfit_1.5-9.10            
 [79] Biostrings_2.74.0           SummarizedExperiment_1.36.0
 [81] RhpcBLASctl_0.23-42         xfun_0.49                  
 [83] statmod_1.5.0               matrixStats_1.4.1          
 [85] KEGGgraph_1.66.0            stringi_1.8.4              
 [87] UCSC.utils_1.2.0            yaml_2.3.10                
 [89] boot_1.3-31                 evaluate_1.0.1             
 [91] codetools_0.2-20            archive_1.1.10             
 [93] tibble_3.2.1                Rgraphviz_2.50.0           
 [95] cli_3.6.3                   RcppParallel_5.1.9         
 [97] xtable_1.8-4                Rdpack_2.6.1               
 [99] munsell_0.5.1               jquerylib_0.1.4            
[101] Rcpp_1.0.13-1               GenomeInfoDb_1.42.0        
[103] EnvStats_3.0.0              dbplyr_2.5.0               
[105] png_0.1-8                   Rfast_2.1.0                
[107] parallel_4.4.1              blob_1.2.4                 
[109] prettyunits_1.2.0           bitops_1.0-9               
[111] mvtnorm_1.3-2               lmerTest_3.1-3             
[113] scales_1.3.0                purrr_1.0.2                
[115] crayon_1.5.3                fANCOVA_0.6-1              
[117] rlang_1.1.4                 EnrichmentBrowser_2.36.0   
[119] KEGGREST_1.46.0