Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Can CIs be "non-significant" when the p-value is "significant"? #584

Closed
shirdekel opened this issue Apr 8, 2023 · 2 comments
Closed

Can CIs be "non-significant" when the p-value is "significant"? #584

shirdekel opened this issue Apr 8, 2023 · 2 comments
Labels
bug 🐜 Something isn't working

Comments

@shirdekel
Copy link

I'm working with 2x2 data whose table looks like this:

Treatment Control
A 13 16
B 6 25

I want to get a risk ratio with CIs but when using effectsize::riskratio() this overlaps with 1, whereas the chi-square test p value was < 0.05. Shouldn't these correspond? PropCIs::riskscoreci(), which is showcased in Agresti's An Introduction to Categorical Data Analysis, 3rd Edition, gives a CI that doesn't overlap with 1 (as expected). Similarly, phi CIs overlap with 0.

Doesn't the calculation in PropCIs::riskscoreci() seem more consistent with the chi-square results?

library(effectsize)
library(PropCIs)

## chi-square test shows that p < 0.05

c(13, 6, 16, 25) |>
  matrix(nrow = 2) |>
  chisq.test(correct = FALSE)
#> 
#>  Pearson's Chi-squared test
#> 
#> data:  matrix(c(13, 6, 16, 25), nrow = 2)
#> X-squared = 4.4929, df = 1, p-value = 0.03404

## effectsize::riskratio() shows that the CI overlaps with 1

c(13, 6, 16, 25) |>
  matrix(nrow = 2) |>
  riskratio()
#> Risk ratio |       95% CI
#> -------------------------
#> 1.75       | [0.71, 4.34]

## PropCIs::riskscoreci() shows that the CI doesn't overlap with 1

riskscoreci(13, 19, 16, 41, conf.level = .95)
#> 
#> 
#> 
#> data:  
#> 
#> 95 percent confidence interval:
#>  1.046215 2.860102

## effectsize::phi() shows that the CI overlaps with 0

c(13, 6, 16, 25) |>
  matrix(nrow = 2) |>
  phi(alternative = "two.sided")
#> Phi (adj.) |       95% CI
#> -------------------------
#> 0.24       | [0.00, 0.51]

Created on 2023-04-08 by the reprex package (v2.0.1)

Session info
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.2.1 (2022-06-23)
#>  os       macOS Monterey 12.6
#>  system   aarch64, darwin20
#>  ui       X11
#>  language (EN)
#>  collate  en_AU.UTF-8
#>  ctype    en_AU.UTF-8
#>  tz       Australia/Sydney
#>  date     2023-04-08
#>  pandoc   2.19.2 @ /opt/homebrew/bin/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version date (UTC) lib source
#>  bayestestR    0.13.0  2022-09-18 [1] CRAN (R 4.2.0)
#>  cli           3.4.1   2022-09-23 [1] CRAN (R 4.2.0)
#>  datawizard    0.6.5   2022-12-14 [1] CRAN (R 4.2.0)
#>  digest        0.6.29  2021-12-01 [1] CRAN (R 4.2.0)
#>  effectsize  * 0.8.3.6 2023-04-08 [1] Github (easystats/effectsize@4c60e28)
#>  evaluate      0.16    2022-08-09 [1] CRAN (R 4.2.0)
#>  fansi         1.0.3   2022-03-24 [1] CRAN (R 4.2.0)
#>  fastmap       1.1.0   2021-01-25 [1] CRAN (R 4.2.0)
#>  fs            1.5.2   2021-12-08 [1] CRAN (R 4.2.0)
#>  glue          1.6.2   2022-02-24 [1] CRAN (R 4.2.0)
#>  highr         0.9     2021-04-16 [1] CRAN (R 4.2.0)
#>  htmltools     0.5.3   2022-07-18 [1] CRAN (R 4.2.0)
#>  insight       0.19.1  2023-03-18 [1] CRAN (R 4.2.0)
#>  knitr         1.40    2022-08-24 [1] CRAN (R 4.2.0)
#>  lifecycle     1.0.3   2022-10-07 [1] CRAN (R 4.2.0)
#>  magrittr      2.0.3   2022-03-30 [1] CRAN (R 4.2.0)
#>  parameters    0.20.2  2023-01-27 [1] CRAN (R 4.2.0)
#>  pillar        1.8.1   2022-08-19 [1] CRAN (R 4.2.0)
#>  pkgconfig     2.0.3   2019-09-22 [1] CRAN (R 4.2.0)
#>  PropCIs     * 0.3-0   2018-02-23 [1] CRAN (R 4.2.0)
#>  purrr         0.3.4   2020-04-17 [1] CRAN (R 4.2.0)
#>  R.cache       0.16.0  2022-07-21 [1] CRAN (R 4.2.0)
#>  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.0  2022-06-28 [1] CRAN (R 4.2.0)
#>  reprex        2.0.1   2021-08-05 [1] CRAN (R 4.2.0)
#>  rlang         1.1.0   2023-03-14 [1] CRAN (R 4.2.0)
#>  rmarkdown     2.14    2022-04-25 [1] CRAN (R 4.2.0)
#>  sessioninfo   1.2.2   2021-12-06 [1] CRAN (R 4.2.0)
#>  stringi       1.7.8   2022-07-11 [1] CRAN (R 4.2.0)
#>  stringr       1.4.1   2022-08-20 [1] CRAN (R 4.2.0)
#>  styler        1.7.0   2022-03-13 [1] CRAN (R 4.2.0)
#>  tibble        3.1.8   2022-07-22 [1] CRAN (R 4.2.0)
#>  utf8          1.2.2   2021-07-24 [1] CRAN (R 4.2.0)
#>  vctrs         0.6.0   2023-03-16 [1] CRAN (R 4.2.1)
#>  withr         2.5.0   2022-03-03 [1] CRAN (R 4.2.0)
#>  xfun          0.33    2022-09-12 [1] CRAN (R 4.2.0)
#>  yaml          2.3.5   2022-02-21 [1] CRAN (R 4.2.0)
#> 
#>  [1] /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────
@mattansb
Copy link
Member

mattansb commented Apr 8, 2023

I think we might have a mistake in the CI method for RR.

@bwiernik we currently have:

SE_logRR <- sqrt(p1 / ((1 - p1) * n1)) + sqrt(p2 / ((1 - p2) * n2))
Z_logRR <- stats::qnorm(alpha / 2, lower.tail = FALSE)
confs <- exp(log(RR) + c(-1, 1) * SE_logRR * Z_logRR)

But I think that first row should be:

SE_logRR <- sqrt((1 - p1) / (n1 * p1) + (1 - p2) / (n2 * p2))

This gives identical results to {PropCIs} as well as {emmeans}:

tab <- c(13, 6, 16, 25) |>
  matrix(nrow = 2)

effectsize::riskratio(tab)
#> Risk ratio |       95% CI
#> -------------------------
#> 1.75       | [1.07, 2.86]

PropCIs::riskscoreci(13, 19, 16, 41, conf.level = .95)
#> 
#> 
#> 
#> data:  
#> 
#> 95 percent confidence interval:
#>  1.046215 2.860102
#>  

d <- tab |> t() |> 
  as.data.frame() |> 
  transform(g = 0:1)

m <- glm(cbind(V1, V2) ~ g, data = d, family = binomial())

emmeans::emmeans(m, ~ g) |> 
  emmeans::regrid(transform = "log") |> 
  emmeans::contrast(method = "pairwise", type = "resp", infer = TRUE)
#>  contrast ratio    SE  df asymp.LCL asymp.UCL null z.ratio p.value
#>  g0 / g1   1.75 0.438 Inf      1.07      2.86    1   2.248  0.0246
#> 
#> Confidence level used: 0.95 
#> Intervals are back-transformed from the log scale 
#> Tests are performed on the log scale

As a side note, since $\chi^2$ test is one tailed, a significance test should match the one-sided CI:

chisq.test(tab, correct = FALSE)
#> 
#> 	Pearson's Chi-squared test
#> 
#> data:  tab
#> X-squared = 4.4929, df = 1, p-value = 0.03404
#> 

effectsize::phi(tab, adjust = FALSE)
#> Phi  |       95% CI
#> -------------------
#> 0.27 | [0.05, 1.00]
#> 
#> - One-sided CIs: upper bound fixed at [1.00].

Also note that Fisher's exact test is not significant here:

fisher.test(tab)
#> 
#> 	Fisher's Exact Test for Count Data
#> 
#> data:  tab
#> p-value = 0.05172
#> alternative hypothesis: true odds ratio is not equal to 1
#> 95 percent confidence interval:
#>   0.9412296 12.9916305
#> sample estimates:
#> odds ratio 
#>   3.314406 
#>   

mattansb added a commit that referenced this issue Apr 8, 2023
@mattansb mattansb added the bug 🐜 Something isn't working label Apr 8, 2023
@mattansb
Copy link
Member

Fixed in #585

mattansb added a commit that referenced this issue Apr 11, 2023
* fix to #584

* news + tests
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug 🐜 Something isn't working
Projects
None yet
Development

No branches or pull requests

2 participants