Incorrect results derived with MultinomCI with Goodman method

[yanka3]

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

[GegznaV]

@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

2) Provide SAS output as formatted code, not as a print screen.

These things would make debugging easier.

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

[GegznaV]

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.

[GegznaV]

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

[AndriSignorell]

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.

[GegznaV]

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

``` r 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 #> #> ────────────────────────────────────────────────────────────────────────────── ```

[AndriSignorell]

Thanks, I'll close then...