Issue with weights in voomLmFit()
1
0
Entering edit mode
@gabrielhoffman-8391
Last seen 28 days ago
United States

Hi Gordon Smyth, I found a subtle error in voomLmFit(). Here is a reproducible example and a fix.

Best, Gabriel

library(edgeR)
library(variancePartition)
library(limma)
packageVersion("limma")
# 3.56.2

# example data is from variancePartition
# but I could have used another example
data(varPartDEdata)

dge = DGEList(counts = countMatrix)
dge = calcNormFactors(dge)

form <- ~ Disease 
dsgn = model.matrix(form, metadata)

# set prior.weights to be 1
w = rep(1, ncol(dge))

# Since all prior weights are 1, these should give 
# identical coefficients
fit1 = voomLmFit( dge, dsgn)
fit2 = voomLmFit( dge, dsgn, prior.weights=w)

# But the coefs are different
coef(fit1)[6,]
#> (Intercept)    Disease1 
#>    4.030300    1.006082
coef(fit2)[6,]
#> (Intercept)    Disease1 
#>   4.0520653   0.9855356

# Since the weights from the first gene are used, the coefs for the first gene are identical, but all other are different

It looks like the weights from the first gene are being used here

In edgeR::voomLmFit(), the function asMatrixWeights() is used which sets attr(weights, "arrayweights") <- TRUE code here

When this is TRUE, this condition is also TRUE, and only the weights from the first gene are used code here

If the code is modified to set attr(weights, "arrayweights") to FALSE, then this loop is used to evaluate one gene at a time using the proper row of weights and the code gives the expected results code here

> sessionInfo( )
R version 4.3.0 (2023-04-21)
Platform: x86_64-apple-darwin22.4.0 (64-bit)
Running under: macOS Ventura 13.5

Matrix products: default
BLAS:   /Users/gabrielhoffman/prog/R-4.3.0/lib/libRblas.dylib 
LAPACK: /usr/local/Cellar/r/4.3.0_1/lib/R/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
[1] variancePartition_1.31.21 BiocParallel_1.34.2      
[3] ggplot2_3.4.2             RUnit_0.4.32             
[5] edgeR_3.42.4              limma_3.56.2             

loaded via a namespace (and not attached):
 [1] gtable_0.3.3        caTools_1.18.2      Biobase_2.60.0     
 [4] lattice_0.21-8      numDeriv_2016.8-1.1 vctrs_0.6.3        
 [7] tools_4.3.0         Rdpack_2.4          bitops_1.0-7       
[10] generics_0.1.3      parallel_4.3.0      pbkrtest_0.5.2     
[13] tibble_3.2.1        fansi_1.0.4         pkgconfig_2.0.3    
[16] Matrix_1.5-4.1      KernSmooth_2.23-21  lifecycle_1.0.3    
[19] stringr_1.5.0       compiler_4.3.0      gplots_3.1.3       
[22] munsell_0.5.0       RhpcBLASctl_0.23-42 codetools_0.2-19   
[25] lmerTest_3.1-3      pillar_1.9.0        nloptr_2.0.3       
[28] tidyr_1.3.0         MASS_7.3-60         aod_1.3.2          
[31] iterators_1.0.14    boot_1.3-28.1       nlme_3.1-162       
[34] gtools_3.9.4        tidyselect_1.2.0    locfit_1.5-9.7     
[37] stringi_1.7.12      fANCOVA_0.6-1       mvtnorm_1.2-2      
[40] reshape2_1.4.4      dplyr_1.1.2         purrr_1.0.1        
[43] splines_4.3.0       grid_4.3.0          colorspace_2.1-0   
[46] cli_3.6.1           magrittr_2.0.3      utf8_1.2.3         
[49] broom_1.0.5         corpcor_1.6.10      withr_2.5.0        
[52] scales_1.2.1        backports_1.4.1     remaCor_0.0.17     
[55] lme4_1.1-33         rbibutils_2.2.13    EnvStats_2.7.0     
[58] rlang_1.1.1         Rcpp_1.0.11         glue_1.6.2         
[61] BiocGenerics_0.46.0 minqa_1.2.5         plyr_1.8.8         
[64] R6_2.5.1
limma • 453 views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 45 minutes ago
WEHI, Melbourne, Australia

Thanks for alerting me to the bug. I've fixed it in edgeR 3.99.2.

It is an usual use-case of course. I've never come across a circumstance where one would give prior sample weights to voom() or voomLmFit(). Nevertheless, the function does allow for that and the code should work correctly.

ADD COMMENT

Login before adding your answer.

Traffic: 745 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6