[Bug]: incidence rate when zero events in one group

[lchensigma]

What happened?

A bug happened!
when zero events in one group, CI of the incidence, incidence ratio, CI and P value are not correct with estimate_incidence_rate and summarize_glm_count.

sessionInfo()

R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 26100)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/Los_Angeles
tzcode source: internal

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

other attached packages:
 [1] tern.mmrm_0.3.2     scda_0.1.6          tern_0.9.6          haven_2.5.4         rtables_0.6.10     
 [6] magrittr_2.0.3      formatters_0.5.9    rstan_2.32.6        StanHeaders_2.32.10 tidybayes_3.0.7    
[11] bayesplot_1.11.1    brms_2.22.0         Rcpp_1.0.13-1       dreamer_3.1.0       DoseFinding_1.2-1  
[16] labelled_2.13.0     gt_0.11.1           flextable_0.9.7     broom_1.0.7         readxl_1.4.3       
[21] boxr_0.3.6          pander_0.6.5        patchwork_1.3.0     ggprism_1.0.5       gtsummary_2.0.3    
[26] ggthemes_5.1.0      knitr_1.49          lubridate_1.9.3     forcats_1.0.0       stringr_1.5.1      
[31] dplyr_1.1.4         purrr_1.0.2         readr_2.1.5         tidyr_1.3.1         tibble_3.2.1       
[36] ggplot2_3.5.1       tidyverse_2.0.0    

loaded via a namespace (and not attached):
  [1] svUnit_1.0.6            shinythemes_1.2.0       splines_4.3.1           later_1.3.2            
  [5] cellranger_1.1.0        xts_0.14.1              forestplot_3.1.5        lifecycle_1.0.4        
  [9] Rdpack_2.6.2            rprojroot_2.0.4         lattice_0.22-6          MASS_7.3-60.0.1        
 [13] crosstalk_1.2.1         ggdist_3.3.2            backports_1.5.0         metafor_4.6-0          
 [17] rmarkdown_2.29          yaml_2.3.10             httpuv_1.6.15           zip_2.3.1              
 [21] askpass_1.2.1           pkgbuild_1.4.5          cowplot_1.1.3           minqa_1.2.8            
 [25] pkgload_1.4.0           multcomp_1.4-26         abind_1.4-8             TH.data_1.1-2          
 [29] tensorA_0.36.2.1        sandwich_3.1-1          gdtools_0.4.1           inline_0.3.20          
 [33] xslt_1.4.6              RBesT_1.7-3             parallelly_1.39.0       bridgesampling_1.1-2   
 [37] codetools_0.2-19        DT_0.33                 xml2_1.3.6              tidyselect_1.2.1       
 [41] rstanarm_2.32.1         farver_2.1.2            lme4_1.1-35.5           matrixStats_1.4.1      
 [45] stats4_4.3.1            base64enc_0.1-3         mathjaxr_1.6-0          jsonlite_1.8.9         
 [49] Formula_1.2-6           emmeans_1.10.5          survival_3.7-0          systemfonts_1.1.0      
 [53] tools_4.3.1             ragg_1.3.3              mmrm_0.3.14             glue_1.8.0             
 [57] gridExtra_2.3           xfun_0.48               here_1.0.1              distributional_0.5.0   
 [61] loo_2.8.0               withr_3.0.2             numDeriv_2016.8-1.1     fastmap_1.2.0          
 [65] boot_1.3-28.1           fansi_1.0.6             shinyjs_2.1.0           openssl_2.2.2          
 [69] digest_0.6.37           estimability_1.5.1      timechange_0.3.0        R6_2.5.1               
 [73] mime_0.12               textshaping_0.4.0       colorspace_2.1-1        gtools_3.9.5           
 [77] markdown_1.13           threejs_0.3.3           utf8_1.2.4              generics_0.1.3         
 [81] fontLiberation_0.1.0    data.table_1.16.2       htmlwidgets_1.6.4       pkgconfig_2.0.3        
 [85] dygraphs_1.1.1.6        gtable_0.3.6            htmltools_0.5.8.1       fontBitstreamVera_0.1.1
 [89] scales_1.3.0            posterior_1.6.0         rstudioapi_0.17.1       tzdb_0.4.0             
 [93] reshape2_1.4.4          uuid_1.2-1              coda_0.19-4.1           checkmate_2.3.2        
 [97] nlme_3.1-166            curl_6.0.1              nloptr_2.1.1            zoo_1.8-12             
[101] parallel_4.3.1          miniUI_0.1.1.1          metadat_1.2-0           pillar_1.9.0           
[105] grid_4.3.1              vctrs_0.6.5             MetaStan_1.0.0          shinystan_2.6.0        
[109] promises_1.3.0          arrayhelpers_1.1-0      xtable_1.8-4            evaluate_1.0.1         
[113] mvtnorm_1.3-2           cli_3.6.3               compiler_4.3.1          rlang_1.1.4            
[117] rstantools_2.4.0        plyr_1.8.9              stringi_1.8.4           QuickJSR_1.4.0         
[121] assertthat_0.2.1        HDInterval_0.2.4        munsell_0.5.1           colourpicker_1.3.0     
[125] Brobdingnag_1.2-9       V8_6.0.0                fontquiver_0.2.1        Matrix_1.6-5           
[129] hms_1.1.3               shiny_1.9.1             rbibutils_2.3           igraph_2.1.1           
[133] RcppParallel_5.1.9      katex_1.5.0             equatags_0.2.1          officer_0.6.7

Relevant log output

Code of Conduct

• I agree to follow this project's Code of Conduct.

Contribution Guidelines

• I agree to follow this project's Contribution Guidelines.

Security Policy

• I agree to follow this project's Security Policy.

[Melkiades]

Hi @lchensigma could you create a minimal example for us to fix your issue? Thanks!!

[lchensigma]

Here is the sample code. No events in group==0. You can test with ref_group ="1" and ref_group ="2"

library(tidyverse)
library(rtables)
library(tern)

dat <- data.frame(TRT01P = rep(1:2, each=50), TRTDURD = rexp(n= 100, rate=1)) |>
  mutate(n_events = ifelse(TRT01P==1,rpois(n=11, lambda = TRTDURD), 0 ),
         lgTRTDURD =log(TRTDURD))
dat



lyt_ir <- basic_table(
  show_colcounts = TRUE,
) %>% 
  split_cols_by(var = "TRT01P", ref_group = "1") %>%
  estimate_incidence_rate(vars = "TRTDURD",
                          n_events = "n_events",
                          table_names = "ir", 
                          .labels = c(rate = "Incidence Rate (per person-year) [2]", n_events = "Total number of events"),
                          .formats = c( person_years= "xx.xxx"),
                          control = control_incidence_rate(
                            input_time_unit = "day",
                            num_pt_year = 1,
                            conf_type = "exact"
                          )
  )  %>%
  summarize_glm_count(
    vars = "n_events",
    variables = list(arm = "TRT01P", offset = "lgTRTDURD", covariates = NULL),
    conf_level = 0.95,
    distribution = "quasipoisson",
    rate_mean_method = "emmeans",
    var_labels = "Unadjusted rate (per person-year) [3]",
    table_names = "unadj",
    .stats = c("rate", "rate_ci", "rate_ratio", "rate_ratio_ci", "pval"),
    .labels = c(
      rate = "Rate", rate_ci = "Rate 95% CI", rate_ratio = "Rate Ratio",
      rate_ratio_ci = "Rate Ratio 95% CI", pval = "p value"
    )
  ) 


result_ir <- build_table(
  lyt = lyt_ir,
  df = dat 
)
result_ir

[shajoezhu]

Hi @lchensigma , your concern is around the p-value right? one would expect a significant difference?

> result_ir
                                               1                 2      
                                             (N=50)           (N=50)    
————————————————————————————————————————————————————————————————————————
Total patient-years at risk                  0.132             0.125    
Total number of events                         18                0      
Incidence Rate (per person-year) [2]         136.06            0.00     
95% CI                                  (80.64, 215.04)    (0.00, 29.43)
Unadjusted rate (per person-year) [3]                                   
  Rate                                       0.3725           0.0000    
    Rate 95% CI                         (0.2455, 0.5654)   (0.0000, Inf)
  Rate Ratio                                                  0.0000    
    Rate Ratio 95% CI                                      (0.0000, Inf)
    p value                                                   0.9931    

[lchensigma]

Also the UL of 95% CI for rate and rate ratio. They can not be infinite.