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
@yanka3, please,
- 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
```
- Provide SAS output as formatted code, not as a print screen.
These things would make debugging easier.
- Please, elaborate more on
TYPE=2in 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.05Created 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.2702358Created 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...