What type of negative binomial does pscl fit for hurdle models?

46 Views Asked by At

I am working with the pscl package to fit hurdle models, and the negative binomial looks adequate to model the count part. As far as I know, the negative binomial can be implemented in two ways: linear relationship between variance and mean (type 1) or quadratic relationship between variance and mean (type 2).

The pscl::hurdle() function does not allow the user to select the type of negative binomial. It appears to implement type 2 by default, based on comparisons with glmmTMB::glmmTMB(), which allows the user to specify the type of negative binomial (see reproducible example below).

I am not particularly concerned about this lack of user control, but I do feel rather curious about whether there is a reason for this. This question has been asked before, but I haven't found an answer.

Cheers!

library(glmmTMB)
library(MuMIn)
#> Warning: package 'MuMIn' was built under R version 4.2.3
library(pscl)
#> Warning: package 'pscl' was built under R version 4.2.3
#> Classes and Methods for R developed in the
#> Political Science Computational Laboratory
#> Department of Political Science
#> Stanford University
#> Simon Jackman
#> hurdle and zeroinfl functions by Achim Zeileis

data("bioChemists", package = "pscl")

list.mods <-
  list(pscl::hurdle(art ~ fem + mar, data = bioChemists, dist = "negbin", zero.dist = "binomial", link = "logit"),
       glmmTMB::glmmTMB(art ~ fem + mar, data = bioChemists, family = truncated_nbinom2(), ziformula = ~ .),
       glmmTMB::glmmTMB(art ~ fem + mar, data = bioChemists, family = truncated_nbinom1(), ziformula = ~ .)
  )


lapply(list.mods, MuMIn::AICc)
#> [[1]]
#> [1] 3220.249
#> 
#> [[2]]
#> [1] 3220.249
#> 
#> [[3]]
#> [1] 3223.829
lapply(list.mods, function(x) coef(summary(x)))
#> [[1]]
#> [[1]]$count
#>                Estimate Std. Error    z value     Pr(>|z|)
#> (Intercept)  0.55487952  0.1192149  4.6544466 3.248520e-06
#> femWomen    -0.27871829  0.1001454 -2.7831354 5.383633e-03
#> marMarried   0.01663515  0.1053245  0.1579419 8.745026e-01
#> Log(theta)   0.32251568  0.2279771  1.4146847 1.571609e-01
#> 
#> [[1]]$zero
#>                Estimate Std. Error   z value     Pr(>|z|)
#> (Intercept)  0.90794254  0.1563982  5.805326 6.424077e-09
#> femWomen    -0.24195114  0.1492268 -1.621365 1.049394e-01
#> marMarried   0.07802262  0.1563561  0.499006 6.177752e-01
#> 
#> 
#> [[2]]
#> [[2]]$cond
#>                Estimate Std. Error    z value     Pr(>|z|)
#> (Intercept)  0.55488673  0.1192145  4.6545246 3.247291e-06
#> femWomen    -0.27872608  0.1001454 -2.7832139 5.382331e-03
#> marMarried   0.01662937  0.1053244  0.1578872 8.745457e-01
#> 
#> [[2]]$zi
#>                Estimate Std. Error    z value     Pr(>|z|)
#> (Intercept) -0.90793049  0.1563979 -5.8052611 6.426577e-09
#> femWomen     0.24195695  0.1492267  1.6214055 1.049307e-01
#> marMarried  -0.07803608  0.1563558 -0.4990931 6.177138e-01
#> 
#> [[2]]$disp
#> NULL
#> 
#> 
#> [[3]]
#> [[3]]$cond
#>                Estimate Std. Error     z value     Pr(>|z|)
#> (Intercept)  0.56636984  0.1443617  3.92326931 8.735541e-05
#> femWomen    -0.28952867  0.1420362 -2.03841500 4.150845e-02
#> marMarried  -0.00213153  0.1446531 -0.01473546 9.882432e-01
#> 
#> [[3]]$zi
#>                Estimate Std. Error    z value     Pr(>|z|)
#> (Intercept) -0.90793602  0.1563981 -5.8052868 6.425591e-09
#> femWomen     0.24194597  0.1492270  1.6213289 1.049471e-01
#> marMarried  -0.07803492  0.1563561 -0.4990846 6.177198e-01
#> 
#> [[3]]$disp
#> NULL

sessionInfo()
#> R version 4.2.0 (2022-04-22 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 22000)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8   
#> [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
#> [5] LC_TIME=Spanish_Spain.utf8    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] pscl_1.5.5.1  MuMIn_1.47.5  glmmTMB_1.1.3
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.10           compiler_4.2.0        nloptr_2.0.1         
#>  [4] highr_0.9             TMB_1.8.1             tools_4.2.0          
#>  [7] boot_1.3-28           digest_0.6.29         lme4_1.1-29          
#> [10] evaluate_0.15         lifecycle_1.0.3       nlme_3.1-157         
#> [13] lattice_0.20-45       rlang_1.1.1           reprex_2.0.2         
#> [16] Matrix_1.4-1          cli_3.6.1             rstudioapi_0.13      
#> [19] yaml_2.3.5            mvtnorm_1.1-3         xfun_0.40            
#> [22] fastmap_1.1.0         coda_0.19-4           withr_2.5.0          
#> [25] stringr_1.5.0         knitr_1.39            fs_1.5.2             
#> [28] vctrs_0.6.3           stats4_4.2.0          grid_4.2.0           
#> [31] glue_1.6.2            survival_3.3-1        rmarkdown_2.14       
#> [34] multcomp_1.4-19       TH.data_1.1-1         minqa_1.2.4          
#> [37] magrittr_2.0.3        codetools_0.2-18      htmltools_0.5.6      
#> [40] emmeans_1.8.5-9000004 splines_4.2.0         MASS_7.3-57          
#> [43] xtable_1.8-4          numDeriv_2016.8-1.1   sandwich_3.0-1       
#> [46] stringi_1.7.6         estimability_1.4.1    zoo_1.8-10

Created on 2023-09-07 with reprex v2.0.2

0

There are 0 best solutions below