AndriSignorell/DescTools

Incorrect results derived with MultinomCI with Goodman method

yanka3 opened this issue · 6 comments

MultinomCI(x = c(91,49,37,43),conf.level = 0.95,method="goodman")
est lwr.ci upr.ci
[1,] 0.4136364 0.3253373 0.5078604
[2,] 0.2227273 0.1545920 0.3098852
[3,] 0.1681818 0.1093614 0.2497670
[4,] 0.1954545 0.1317168 0.2800860

However, the original paper indicates the outputs should be
image

@yanka3, please,

  1. format R code and its output correctly, e.g., use reprex. Otherwise, it is sufficient to write the code) between backtick as in this example:
```r
# Your code
```
  1. Provide SAS output as formatted code, not as a print screen.

These things would make debugging easier.

  1. Please, elaborate more on TYPE=2 in SAS code. Does it define that Goodman's CI is used here?

I checked Python's statsmodels implementation:

import statsmodels.stats.api as sms

counts = [91, 49, 37, 43]
sms.multinomial_proportions_confint(counts, method="goodman")
#> array([[0.33420248, 0.49783321],
#>        [0.16085875, 0.2998874 ],
#>        [0.11455134, 0.24011208],
#>        [0.13746901, 0.27023577]])

Here are only CIs without the estimate of proportions.
And the results are the same as the SAS output reported by @yanka3.

@yanka3, thank you for reporting this.

From R package scimple:

scimple::scimple_ci(c(91,49,37,43), alpha=0.05, methods="goodman")
#>    method lower_limit upper_limit    adj_ll    adj_ul     volume inpmat alpha
#> 1 goodman   0.3342025   0.4978332 0.3342025 0.4978332 0.00037924     91  0.05
#> 2 goodman   0.1608587   0.2998874 0.1608587 0.2998874 0.00037924     49  0.05
#> 3 goodman   0.1145513   0.2401121 0.1145513 0.2401121 0.00037924     37  0.05
#> 4 goodman   0.1374690   0.2702358 0.1374690 0.2702358 0.00037924     43  0.05

Created on 2023-04-16 with reprex v2.0.2

Here is the source:
https://github.com/hrbrmstr/scimple/blob/948d6d7fe4a9de67c3009d68e5546645ff0e127c/R/goodman.r

Thank you all for your analysis. This was an error in how I had defined the chisquare quantile. I fixed that (v. 0.99.48.3) and the function now yields the same results as SAS mentioned by yanka3.
Please close the issue if you can confirm the solution.

It seems OK now.

DescTools::MultinomCI(x = c(91, 49, 37, 43), conf.level = 0.95, method = "goodman")
#>            est    lwr.ci    upr.ci
#> [1,] 0.4136364 0.3342025 0.4978332
#> [2,] 0.2227273 0.1608587 0.2998874
#> [3,] 0.1681818 0.1145513 0.2401121
#> [4,] 0.1954545 0.1374690 0.2702358

Created on 2023-04-20 with reprex v2.0.2

Session info
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.2.3 (2023-03-15 ucrt)
#>  os       Windows 10 x64 (build 19045)
#>  system   x86_64, mingw32
#>  ui       RTerm
#>  language (EN)
#>  collate  English_United States.utf8
#>  ctype    English_United States.utf8
#>  tz       Europe/Helsinki
#>  date     2023-04-20
#>  pandoc   2.19.2 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version   date (UTC) lib source
#>  boot          1.3-28.1  2022-11-22 [1] CRAN (R 4.2.2)
#>  cellranger    1.1.0     2016-07-27 [1] CRAN (R 4.2.0)
#>  class         7.3-21    2023-01-23 [2] CRAN (R 4.2.3)
#>  cli           3.6.1     2023-03-23 [1] CRAN (R 4.2.3)
#>  data.table    1.14.8    2023-02-17 [1] CRAN (R 4.2.2)
#>  DescTools     0.99.48.3 2023-04-20 [1] Github (AndriSignorell/DescTools@46a6e43)
#>  digest        0.6.31    2022-12-11 [1] CRAN (R 4.2.2)
#>  e1071         1.7-13    2023-02-01 [1] CRAN (R 4.2.2)
#>  evaluate      0.20      2023-01-17 [1] CRAN (R 4.2.2)
#>  Exact         3.2       2022-09-25 [1] CRAN (R 4.2.1)
#>  expm          0.999-7   2023-01-09 [1] CRAN (R 4.2.2)
#>  fastmap       1.1.1     2023-02-24 [1] CRAN (R 4.2.2)
#>  fs            1.6.1     2023-02-06 [1] CRAN (R 4.2.2)
#>  gld           2.6.6     2022-10-23 [1] CRAN (R 4.2.1)
#>  glue          1.6.2     2022-02-24 [1] CRAN (R 4.2.0)
#>  htmltools     0.5.5     2023-03-23 [1] CRAN (R 4.2.3)
#>  httr          1.4.5     2023-02-24 [1] CRAN (R 4.2.2)
#>  knitr         1.42      2023-01-25 [1] CRAN (R 4.2.2)
#>  lattice       0.20-45   2021-09-22 [2] CRAN (R 4.2.3)
#>  lifecycle     1.0.3     2022-10-07 [1] CRAN (R 4.2.1)
#>  lmom          2.9       2022-05-29 [1] CRAN (R 4.2.0)
#>  magrittr      2.0.3     2022-03-30 [1] CRAN (R 4.2.0)
#>  MASS          7.3-58.3  2023-03-07 [1] CRAN (R 4.2.2)
#>  Matrix        1.5-4     2023-04-04 [1] CRAN (R 4.2.3)
#>  mvtnorm       1.1-3     2021-10-08 [1] CRAN (R 4.2.0)
#>  proxy         0.4-27    2022-06-09 [1] CRAN (R 4.2.0)
#>  purrr         1.0.1     2023-01-10 [1] CRAN (R 4.2.2)
#>  R.cache       0.16.0    2022-07-21 [1] CRAN (R 4.2.1)
#>  R.methodsS3   1.8.2     2022-06-13 [1] CRAN (R 4.2.0)
#>  R.oo          1.25.0    2022-06-12 [1] CRAN (R 4.2.0)
#>  R.utils       2.12.2    2022-11-11 [1] CRAN (R 4.2.2)
#>  R6            2.5.1     2021-08-19 [1] CRAN (R 4.2.0)
#>  Rcpp          1.0.10    2023-01-22 [1] CRAN (R 4.2.2)
#>  readxl        1.4.2     2023-02-09 [1] CRAN (R 4.2.2)
#>  reprex        2.0.2     2022-08-17 [1] CRAN (R 4.2.1)
#>  rlang         1.1.0     2023-03-14 [1] CRAN (R 4.2.2)
#>  rmarkdown     2.21      2023-03-26 [1] CRAN (R 4.2.3)
#>  rootSolve     1.8.2.3   2021-09-29 [1] CRAN (R 4.2.0)
#>  rstudioapi    0.14      2022-08-22 [1] CRAN (R 4.2.1)
#>  sessioninfo   1.2.2     2021-12-06 [1] CRAN (R 4.2.0)
#>  styler        1.9.1     2023-03-04 [1] CRAN (R 4.2.2)
#>  vctrs         0.6.2     2023-04-19 [1] CRAN (R 4.2.3)
#>  withr         2.5.0     2022-03-03 [1] CRAN (R 4.2.0)
#>  xfun          0.38      2023-03-24 [1] CRAN (R 4.2.3)
#>  yaml          2.3.7     2023-01-23 [1] CRAN (R 4.2.2)
#> 
#>  [1] C:/Users/user/AppData/Local/R/win-library/4.2
#>  [2] C:/Program Files/R/R-4.2.3/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────

Thanks, I'll close then...