implement SPI

[strengejacke]

#521

This does not yet change all the docs, I think this would result in too many conflicts, probably. So we might mention the "SPI" option in the docs later.

[bwiernik]

Do we want to copy the function from SPIn rather than soft-importing, considering that the package hasn't been updated since 2013?

[strengejacke]

I wasn't sure about that, but this can be done easily. We just need the code in .spi(). It's more a moral question to me about copying the code 1:1 ;-)

[bwiernik]

SPIn is licensed GPL>=2, so we are contractually in the clear so long as we attribute the source.

I think we should do this considering that the CRAN version of SPIn generates several NOTEs based on changes to CRAN policies that have happened since it was last updated.

https://www.r-project.org/nosvn/R.check/r-devel-linux-x86_64-debian-clang/SPIn-00check.html

[strengejacke]

That method or the code doesn't seem to be very robust:

library(bayestestR)
x <- rnorm(1000)
spi(x) |> as.data.frame()
#>     CI    CI_low  CI_high
#> 1 0.95 -1.851012 2.079607
SPIn::SPIn(x)$spin
#> [1] -1.851012  2.079607

m <- insight::download_model("brms_1")
d <- as.data.frame(m)
sapply(d, spi)
#> Error in quadprog::solve.QP(D.u, d.u, A.u, c(1, rep(0, u.u - u.l + 2)), : constraints are inconsistent, no solution!

sapply(d, SPIn::SPIn)
#> Error in solve.QP(D.u, d.u, A.u, c(1, rep(0, u.u - u.l + 2)), u.u - u.l): constraints are inconsistent, no solution!

^{Created on 2022-04-04 by the reprex package (v2.0.1)}

[bwiernik]

The constraints are the A.l matrix. I would need to look more into the code for how this is generated and the paper this is based on to understand what's going on

[DominiqueMakowski]

should I submit without this PR or wait for a bit?

[DominiqueMakowski]

should I submit without this PR or wait for a bit?

bump 😬

[strengejacke]

w/o this PR, I'd say we must first take a closer look and see how this function can be made more stable. The original SPIn() in the SPIn package has some more arguments, which might help adjusting the algorithm to not fail in the above example - but it doesn't help user for whom it is probably still too opaque at the moment what they could do.

[strengejacke]

But you could elaborate a bit on #519 :-)

[bwiernik]

@strengejacke The issue that's causing it to fail is the bandwidth argument bw is somewhat too large for the intercept and lp_ parameters.

Let's try this:

1. Use tryCatch() with the function using the default bw = round((sqrt(n.sims) - 1)/2)
2. If that errors, try again multiplying that bw by .9, .8, then .7

That should make it safe for most distributions.

[bwiernik]

(I'll take care of that)

[bwiernik]

library(bayestestR)
m <- insight::download_model("brms_1")
d <- as.data.frame(m)
sapply(d, spi)
# > sapply(d, spi)
#         b_Intercept b_wt      b_cyl      sigma    lp__     
# CI      0.95        0.95      0.95       0.95     0.95     
# CI_low  36.25617    -4.703458 -2.386442  1.985231 -83.11161
# CI_high 43.31241    -1.567123 -0.6765242 3.385691 -77.98328

[bwiernik]

Assuming the checks pass, we can include this

[strengejacke]

@DominiqueMakowski have you already submitted? If yes, it's no big deal, but if not, we may include this PR as well - we still need to add this to the docs as well.

[strengejacke]

I don't know, this test seems to me like the results are a bit random...

https://github.com/easystats/bayestestR/runs/5856422169?check_suite_focus=true#step:10:336

[bwiernik]

That's seems very random. We might want to try to stabilize the seed on that test? But seems unrelated to the current changes

[strengejacke]

@bwiernik There are some edge cases, where this if-statement is TRUE, and then x.l was not initialized. I set it to NA, is that ok?

[strengejacke]

Same here. Try spi(distribution_normal(1000), ci = 0.0000001), in which case this if-statement is TRUE, and x.u wasn't initialized later when we return the data frame.

[bwiernik]

What if statement do you mean?

[strengejacke]

We have

  # upper bound
  if (!is.na(ui) && all(ui == n.sims)) {
    x2 <- u
    w.u <- 1
    x.u <- NA # <- this here is new
  } else {

and at the bottom of the function:

  # output
  data.frame(
    "CI" = ci,
    "CI_low" = x.l,
    "CI_high" = x.u
  )

when if (!is.na(ui) && all(ui == n.sims)) { was TRUE, x.u wasn't initialized, so I added x.u <- NA # <- this here is new to the if-statement. I also did some minor fixes based on the tests I added, and fixed all errors/issues I found. Mostly rare edge cases, as you can see from the tests. But in general, it now works well with models.

[strengejacke]

See my last changes here: a22371d
Based on the tests I copied from the test-hdi.R file. This showed some issues I fixed in that commit.

[strengejacke]

@IndrajeetPatil why does test-coverage always fail? Is there something wrong with the script/action?

[bwiernik]

Oh, no. The x1 and x2 variables assigned above should just be changed to x.l and x.u. Missed updating those

[strengejacke]

Cool, I think we're set.

[bwiernik]

Not sure what's up with those tests