Factors with missing levels produce criptic error message

[william-hutchison]

Hello, I have found sccomp produces an error if any factors in the sccomp_estimate() formula have levels which are not represented in the data.

This error can be reproduced with the following code:

library(dplyr)
library(tidyr)
library(sccomp)
library(forcats)
data("seurat_obj")
data("sce_obj")
data("counts_obj")
 
counts_obj = 
  counts_obj |>
  mutate(count = count+1) |> 
  with_groups(sample, ~ .x |> mutate(proportion = count/sum(count))) |>
  mutate(type = factor(type, levels = c("benign", "cancer", "hi")))

counts_obj_results <- 
    counts_obj |>
        sccomp_estimate(
        formula_composition = ~ 0 + type , 
        .sample = sample,  
        .cell_group = cell_group, 
        .count = proportion
    )

and produces the error:

Chain 6: Rejecting initial value:
Chain 6:   Error evaluating the log probability at the initial value.
Chain 6: Exception: Exception: normal_lpdf: Location parameter[1] is -nan, but must be finite! (in 'glm_multi_beta_binomial', line 66, column 6 to column 91) (in 'glm_multi_beta_binomial', line 336, column 6 to line 343, column 10)
Chain 6: 
Chain 6: Initialization between (-2, 2) failed after 100 attempts. 
Chain 6:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error : Initialization failed."
error occurred during calling the sampler; sampling not done
here are whatever error messages were returned
[[1]]
Stan model 'glm_multi_beta_binomial' does not contain samples.

[[2]]
Stan model 'glm_multi_beta_binomial' does not contain samples.

[[3]]
Stan model 'glm_multi_beta_binomial' does not contain samples.

[[4]]
Stan model 'glm_multi_beta_binomial' does not contain samples.

[[5]]
Stan model 'glm_multi_beta_binomial' does not contain samples.

[[6]]
Stan model 'glm_multi_beta_binomial' does not contain samples.

Error in dimnames(x)$parameters : 
  $ operator is invalid for atomic vectors

sessionInfo()
R version 4.3.3 (2024-02-29)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /stornext/Home/data/allstaff/h/hutchison.w/Environments/micromamba/lib/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [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: Australia/Melbourne
tzcode source: system (glibc)

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

other attached packages:
[1] forcats_1.0.0 sccomp_1.7.15 tidyr_1.3.1   dplyr_1.1.4  

loaded via a namespace (and not attached):
 [1] dotCall64_1.1-1             SummarizedExperiment_1.32.0
 [3] gtable_0.3.5                spam_2.10-0                
 [5] QuickJSR_1.1.3              ggplot2_3.5.1              
 [7] ggrepel_0.9.5               inline_0.3.19              
 [9] Biobase_2.62.0              lattice_0.22-6             
[11] tzdb_0.4.0                  vctrs_0.6.5                
[13] tools_4.3.3                 bitops_1.0-7               
[15] generics_0.1.3              stats4_4.3.3               
[17] parallel_4.3.3              tibble_3.2.1               
[19] fansi_1.0.6                 pkgconfig_2.0.3            
[21] Matrix_1.6-5                S4Vectors_0.40.2           
[23] RcppParallel_5.1.7          lifecycle_1.0.4            
[25] GenomeInfoDbData_1.2.11     stringr_1.5.1              
[27] compiler_4.3.3              munsell_0.5.1              
[29] codetools_0.2-20            SeuratObject_5.0.1         
[31] GenomeInfoDb_1.38.8         RCurl_1.98-1.14            
[33] pillar_1.9.0                crayon_1.5.2               
[35] SingleCellExperiment_1.24.0 DelayedArray_0.28.0        
[37] StanHeaders_2.32.6          boot_1.3-30                
[39] abind_1.4-5                 parallelly_1.37.1          
[41] rstan_2.32.6                tidyselect_1.2.1           
[43] digest_0.6.35               stringi_1.8.4              
[45] future_1.33.2               purrr_1.0.2                
[47] listenv_0.9.1               grid_4.3.3                 
[49] colorspace_2.1-0            cli_3.6.2                  
[51] SparseArray_1.2.4           magrittr_2.0.3             
[53] loo_2.7.0                   patchwork_1.2.0            
[55] S4Arrays_1.2.1              pkgbuild_1.4.4             
[57] utf8_1.2.4                  future.apply_1.11.2        
[59] withr_3.0.0                 readr_2.1.5                
[61] scales_1.3.0                sp_2.1-3                   
[63] XVector_0.42.0              matrixStats_1.3.0          
[65] globals_0.16.3              gridExtra_2.3              
[67] progressr_0.14.0            hms_1.1.3                  
[69] GenomicRanges_1.54.1        IRanges_2.36.0             
[71] rstantools_2.4.0            rlang_1.1.3                
[73] Rcpp_1.0.12                 glue_1.7.0                 
[75] BiocGenerics_0.48.1         jsonlite_1.8.8             
[77] R6_2.5.1                    MatrixGenerics_1.14.0      
[79] zlibbioc_1.48.2   

Thanks!