R语言 用估计的边际平均值解释栅栏模型

iyzzxitl  于 2023-09-27  发布在  其他
关注(0)|答案(1)|浏览(117)

我的目标是通过估计的边际均值来解释障碍模型的系数。我更喜欢解释概率（从logit量表反向转换），而不是对数优势（模型系数）或比值比（exp（对数优势））。我想使用emmeans()来实现这个目标，因为它与许多模型兼容，并且我有在线性模型和二项式GLM中使用它的经验。

栅栏模型的栅栏部分产生与二项式GLM相同的系数，这在其他地方已经评论过（here或here）。
然而，我不完全理解为什么，根据设置，从跨栏估计的平均值不匹配的二项式GLM。特别是

  • lin.pred = FALSE。这从二项式GLM产生不同的概率。在emmeansdocumentation之后，我认为它们是在概率尺度上平均的。
  • lin.pred = TRUE, type = "response"。此设置产生与二项式GLM相同的概率。根据emmeans文档，我认为它们是在logit标度上平均的。

我有两个问题：

  1. 1.如果我们的目标是解释概率，而不是logits，那么哪种设置更可取？
  2. 1.Inf自由度的含义是什么？
    下面是一个可重复的例子，使用生物化学学生在至少3年的博士学位期间发表的论文数量，重点是两个因素水平（单身与已婚）之间的成对比较。
    我的首选方法是与二项式GLM的估计概率相匹配的方法。我的理解是，已婚博士生比单身博士生更有可能发表至少一篇论文。已婚学生有74%的机会至少完成一篇论文，而单身学生有67%的机会：这7%的变化是一个很小的差异，并不显著。
    提前感谢！
  1. library(emmeans)
  2. library(pscl)
  3. #> Warning: package 'pscl' was built under R version 4.2.3
  4. #> Classes and Methods for R developed in the
  5. #> Political Science Computational Laboratory
  6. #> Department of Political Science
  7. #> Stanford University
  8. #> Simon Jackman
  9. #> hurdle and zeroinfl functions by Achim Zeileis
  10. data("bioChemists", package = "pscl")
  11. # str(bioChemists)
  12. hurdle <- hurdle(art ~ ., data = bioChemists, dist = "poisson", zero.dist = "binomial", link = "logit")
  13. glm <- glm(art > 0 ~ ., data = bioChemists, family = binomial(link = "logit"))
  15. # Same coefficients
  16. coef(summary(hurdle))$zero
  17. #>                Estimate Std. Error    z value     Pr(>|z|)
  18. #> (Intercept)  0.23679601 0.29551883  0.8012891 4.229643e-01
  19. #> femWomen    -0.25115113 0.15910522 -1.5785222 1.144457e-01
  20. #> marMarried   0.32623358 0.18081823  1.8042074 7.119880e-02
  21. #> kid5        -0.28524872 0.11113043 -2.5667921 1.026441e-02
  22. #> phd          0.02221940 0.07955710  0.2792887 7.800233e-01
  23. #> ment         0.08012135 0.01301763  6.1548321 7.515710e-10
  24. coef(glm)
  25. #> (Intercept)    femWomen  marMarried        kid5         phd        ment 
  26. #>  0.23679601 -0.25115113  0.32623358 -0.28524872  0.02221940  0.08012135
  28. ####### 
  29. ## lin.pred = FALSE
  30. emmeans::emmeans(hurdle, specs = pairwise ~ mar, mode = "zero", lin.pred = FALSE)
  31. #> $emmeans
  32. #>  mar     emmean    SE  df lower.CL upper.CL
  33. #>  Single   0.806 0.167 903    0.478     1.13
  34. #>  Married  0.858 0.122 903    0.619     1.10
  35. #> 
  36. #> Results are averaged over the levels of: fem 
  37. #> Confidence level used: 0.95 
  38. #> 
  39. #> $contrasts
  40. #>  contrast         estimate    SE  df t.ratio p.value
  41. #>  Single - Married  -0.0522 0.213 903  -0.245  0.8066
  42. #> 
  43. #> Results are averaged over the levels of: fem
  45. ####### 
  46. ## lin.pred = TRUE, type = "response"
  47. emmeans::emmeans(hurdle, specs = pairwise ~ mar, mode = "zero", lin.pred = TRUE, type = "response")
  48. #> $emmeans
  49. #>  mar      prob     SE  df lower.CL upper.CL
  50. #>  Single  0.677 0.0306 903    0.615    0.734
  51. #>  Married 0.744 0.0200 903    0.703    0.781
  52. #> 
  53. #> Results are averaged over the levels of: fem 
  54. #> Confidence level used: 0.95 
  55. #> Intervals are back-transformed from the logit scale 
  56. #> 
  57. #> $contrasts
  58. #>  contrast         odds.ratio   SE  df null t.ratio p.value
  59. #>  Single / Married      0.722 0.13 903    1  -1.804  0.0715
  60. #> 
  61. #> Results are averaged over the levels of: fem 
  62. #> Tests are performed on the log odds ratio scale
  64. emmeans::emmeans(glm, specs = pairwise ~ mar, regrid = "response")
  65. #> $emmeans
  66. #>  mar     response     SE  df asymp.LCL asymp.UCL
  67. #>  Single     0.677 0.0305 Inf     0.617     0.736
  68. #>  Married    0.743 0.0201 Inf     0.704     0.783
  69. #> 
  70. #> Results are averaged over the levels of: fem 
  71. #> Confidence level used: 0.95 
  72. #> 
  73. #> $contrasts
  74. #>  contrast         estimate     SE  df z.ratio p.value
  75. #>  Single - Married  -0.0667 0.0377 Inf  -1.772  0.0764
  76. #> 
  77. #> Results are averaged over the levels of: fem
  79. #######
  80. # Estimates and contrasts at the logit scale appear to match
  81. emmeans::emmeans(hurdle, specs = pairwise ~ mar, mode = "zero", lin.pred = TRUE)
  82. #> $emmeans
  83. #>  mar     emmean    SE  df lower.CL upper.CL
  84. #>  Single   0.741 0.140 903    0.466     1.02
  85. #>  Married  1.068 0.105 903    0.862     1.27
  86. #> 
  87. #> Results are averaged over the levels of: fem 
  88. #> Results are given on the logit (not the response) scale. 
  89. #> Confidence level used: 0.95 
  90. #> 
  91. #> $contrasts
  92. #>  contrast         estimate    SE  df t.ratio p.value
  93. #>  Single - Married   -0.326 0.181 903  -1.804  0.0715
  94. #> 
  95. #> Results are averaged over the levels of: fem 
  96. #> Results are given on the log odds ratio (not the response) scale.
  98. emmeans::emmeans(glm, specs = pairwise ~ mar)
  99. #> $emmeans
  100. #>  mar     emmean    SE  df asymp.LCL asymp.UCL
  101. #>  Single   0.741 0.140 Inf     0.467      1.02
  102. #>  Married  1.068 0.105 Inf     0.862      1.27
  103. #> 
  104. #> Results are averaged over the levels of: fem 
  105. #> Results are given on the logit (not the response) scale. 
  106. #> Confidence level used: 0.95 
  107. #> 
  108. #> $contrasts
  109. #>  contrast         estimate    SE  df z.ratio p.value
  110. #>  Single - Married   -0.326 0.181 Inf  -1.804  0.0712
  111. #> 
  112. #> Results are averaged over the levels of: fem 
  113. #> Results are given on the log odds ratio (not the response) scale.
  115. sessionInfo()
  116. #> R version 4.2.0 (2022-04-22 ucrt)
  117. #> Platform: x86_64-w64-mingw32/x64 (64-bit)
  118. #> Running under: Windows 10 x64 (build 22000)
  119. #> 
  120. #> Matrix products: default
  121. #> 
  122. #> locale:
  123. #> [1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8   
  124. #> [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
  125. #> [5] LC_TIME=Spanish_Spain.utf8    
  126. #> 
  127. #> attached base packages:
  128. #> [1] stats     graphics  grDevices utils     datasets  methods   base     
  129. #> 
  130. #> other attached packages:
  131. #> [1] pscl_1.5.5.1          emmeans_1.8.5-9000004
  132. #> 
  133. #> loaded via a namespace (and not attached):
  134. #>  [1] compiler_4.2.0      highr_0.9           tools_4.2.0        
  135. #>  [4] digest_0.6.29       evaluate_0.15       lifecycle_1.0.3    
  136. #>  [7] lattice_0.20-45     rlang_1.1.1         reprex_2.0.2       
  137. #> [10] Matrix_1.4-1        cli_3.6.1           rstudioapi_0.13    
  138. #> [13] yaml_2.3.5          mvtnorm_1.1-3       xfun_0.40          
  139. #> [16] fastmap_1.1.0       coda_0.19-4         withr_2.5.0        
  140. #> [19] stringr_1.5.0       knitr_1.39          fs_1.5.2           
  141. #> [22] vctrs_0.6.3         grid_4.2.0          glue_1.6.2         
  142. #> [25] survival_3.3-1      rmarkdown_2.14      multcomp_1.4-19    
  143. #> [28] TH.data_1.1-1       magrittr_2.0.3      codetools_0.2-18   
  144. #> [31] htmltools_0.5.6     splines_4.2.0       MASS_7.3-57        
  145. #> [34] xtable_1.8-4        numDeriv_2016.8-1.1 sandwich_3.0-1     
  146. #> [37] estimability_1.4.1  stringi_1.7.6       zoo_1.8-10

创建于2023-09-06使用reprex v2.0.2

r

1条答案

按热度按时间
wn9m85ua

wn9m85ua1#

Russell Lenth（emmeans包的开发者）在GitHub上提供了一个答案。我把它贴在这里，用emmeans和glmmTMB拟合的障碍模型之间的比较，显示了一致的结果。

  1. library(emmeans)
  2. #> Warning: package 'emmeans' was built under R version 4.2.3
  3. library(glmmTMB)
  4. #> Warning: package 'glmmTMB' was built under R version 4.2.3
  5. #> Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
  6. #> TMB was built with Matrix version 1.6.1
  7. #> Current Matrix version is 1.4.1
  8. #> Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
  9. library(pscl)
  10. #> Warning: package 'pscl' was built under R version 4.2.3
  11. #> Classes and Methods for R developed in the
  12. #> Political Science Computational Laboratory
  13. #> Department of Political Science
  14. #> Stanford University
  15. #> Simon Jackman
  16. #> hurdle and zeroinfl functions by Achim Zeileis
  18. # Create datasets: number of papers during a Phd
  19. data("bioChemists", package = "pscl")
  21. # Declare models
  22. pscl.hurdle <- pscl::hurdle(art ~ fem + mar, data = bioChemists, dist = "poisson", zero.dist = "binomial", link = "logit")
  23. glmmtmb.hurdle <- glmmTMB::glmmTMB(art ~ fem + mar, data = bioChemists, family = truncated_poisson(link = "log"), zi = ~ fem + mar)
  25. # Truncated count
  26. emmeans(pscl.hurdle, ~ mar, mode = "count", lin.pred = TRUE, type = "response")
  27. #>  mar     count     SE  df lower.CL upper.CL
  28. #>  Single   2.07 0.1122 909     1.86     2.31
  29. #>  Married  2.09 0.0807 909     1.94     2.26
  30. #> 
  31. #> Results are averaged over the levels of: fem 
  32. #> Confidence level used: 0.95 
  33. #> Intervals are back-transformed from the log scale
  34. emmeans(glmmtmb.hurdle, ~ mar, comp = "cond", type = "response")
  35. #>  mar     rate     SE  df asymp.LCL asymp.UCL
  36. #>  Single  2.07 0.1122 Inf      1.86      2.30
  37. #>  Married 2.09 0.0807 Inf      1.94      2.26
  38. #> 
  39. #> Results are averaged over the levels of: fem 
  40. #> Confidence level used: 0.95 
  41. #> Intervals are back-transformed from the log scale
  43. # Binomial
  44. emmeans(pscl.hurdle, ~ mar, mode = "prob0")
  45. #>  mar     emmean     SE  df lower.CL upper.CL
  46. #>  Single   0.313 0.0265 909    0.261    0.365
  47. #>  Married  0.297 0.0191 909    0.259    0.334
  48. #> 
  49. #> Results are averaged over the levels of: fem 
  50. #> Confidence level used: 0.95
  51. emmeans(glmmtmb.hurdle, ~ mar, comp = "zi", type = "response")
  52. #>  mar     response     SE  df asymp.LCL asymp.UCL
  53. #>  Single     0.313 0.0267 Inf     0.263     0.367
  54. #>  Married    0.296 0.0190 Inf     0.261     0.335
  55. #> 
  56. #> Results are averaged over the levels of: fem 
  57. #> Confidence level used: 0.95 
  58. #> Intervals are back-transformed from the logit scale
  60. # Response (binomial * truncated count)
  61. emmeans(pscl.hurdle, ~ mar, mode = "response")
  62. #>  mar     emmean     SE  df lower.CL upper.CL
  63. #>  Single    1.64 0.0900 909     1.47     1.82
  64. #>  Married   1.69 0.0633 909     1.57     1.82
  65. #> 
  66. #> Results are averaged over the levels of: fem 
  67. #> Confidence level used: 0.95
  68. emmeans(glmmtmb.hurdle, ~ mar, comp = "response")
  69. #>  mar     emmean     SE  df asymp.LCL asymp.UCL
  70. #>  Single    1.64 0.0900 Inf      1.47      1.82
  71. #>  Married   1.69 0.0633 Inf      1.57      1.82
  72. #> 
  73. #> Results are averaged over the levels of: fem 
  74. #> Confidence level used: 0.95
  76. sessionInfo()
  77. #> R version 4.2.0 (2022-04-22 ucrt)
  78. #> Platform: x86_64-w64-mingw32/x64 (64-bit)
  79. #> Running under: Windows 10 x64 (build 22000)
  80. #> 
  81. #> Matrix products: default
  82. #> 
  83. #> locale:
  84. #> [1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8   
  85. #> [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
  86. #> [5] LC_TIME=Spanish_Spain.utf8    
  87. #> 
  88. #> attached base packages:
  89. #> [1] stats     graphics  grDevices utils     datasets  methods   base     
  90. #> 
  91. #> other attached packages:
  92. #> [1] pscl_1.5.5.1  glmmTMB_1.1.7 emmeans_1.8.8
  93. #> 
  94. #> loaded via a namespace (and not attached):
  95. #>  [1] Rcpp_1.0.10         nloptr_2.0.1        compiler_4.2.0     
  96. #>  [4] highr_0.9           TMB_1.9.6           tools_4.2.0        
  97. #>  [7] boot_1.3-28         lme4_1.1-29         digest_0.6.29      
  98. #> [10] nlme_3.1-157        evaluate_0.15       lifecycle_1.0.3    
  99. #> [13] lattice_0.20-45     rlang_1.1.1         reprex_2.0.2       
  100. #> [16] Matrix_1.4-1        cli_3.6.1           rstudioapi_0.13    
  101. #> [19] yaml_2.3.5          mvtnorm_1.1-3       xfun_0.40          
  102. #> [22] fastmap_1.1.0       coda_0.19-4         withr_2.5.0        
  103. #> [25] stringr_1.5.0       knitr_1.39          fs_1.5.2           
  104. #> [28] vctrs_0.6.3         grid_4.2.0          glue_1.6.2         
  105. #> [31] survival_3.3-1      rmarkdown_2.14      multcomp_1.4-19    
  106. #> [34] minqa_1.2.4         TH.data_1.1-1       magrittr_2.0.3     
  107. #> [37] codetools_0.2-18    htmltools_0.5.6     splines_4.2.0      
  108. #> [40] MASS_7.3-57         xtable_1.8-4        numDeriv_2016.8-1.1
  109. #> [43] sandwich_3.0-1      estimability_1.4.1  stringi_1.7.6      
  110. #> [46] zoo_1.8-10

创建于2023-09-14带有reprex v2.0.2

展开查看全部
赞(0）分享  回复(0）举报  2023-09-27
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备13028337号-1 大数据知识库 https://www.saoniuhuo.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com