Robustness analysis for informed Bayesian t-tests

Adventures in WASM and Quarto

Author

Affiliation

Lincoln Colling

School of Psychology, University of Sussex

tldr

Currently, robustness analysis for informed Bayesian t-tests is not implemented in JASP. So I thought I’d just build it! This is also a great excuse to play around with Rust WASM and Observable blocks in Quarto.

import {wasm} from "@ljcolling/robust"

Motivation

Feel free to skip this section if you’re not interested in the technical details. You can just jump straight to the calculator. I recently released the bayesplay Web-app and R Package for computing Bayes factors. At the prompting of a colleague, I decided to include Robustness analysis as a feature of this package. This was fairly straightforward to do, and I was able to quickly include it in the R package. However, when it came to including it in the web-app I hit a snag. The snag was that it was very slow. The common way to perform robustness analysis is to generate thousands of Bayes factor value for slightly different values of the prior parameters. The web-app currently performed all the computations in Go compiled to WASM. I was running into memory issues and fighting against the GC, and this was slowing it down. So I thought to myself, “this is a great excuse to write some RUST!” (Because that’s what all the cool kids think!). I didn’t want to build out a complete app just for benchmarking, so I had a second thought. I love ObservableHQ for building quick mock-ups of things. And RMarkdown (or rather Quarto) now has native support for Observable, so I thought I’d just make a blog post instead. And to keep things even simpler, I thought that instead of implementing robustness analysis for all the models that bayesplay supports, I’d just do it for two-sample informed Bayesian t-tests (e.g., Gronau et al., 2020). I picked this, because you can run these models in JASP. Well, almost, when you click the Bayes factor robustness check plot tick-box you get a message telling you: “Bayes factor robustness check plot currently not supported for informed prior”. So that’s the aim, implement this missing feature of JASP in WASM, and Observable.

Some design choices

For the calculator, I’ve aimed to replicate the JASP interface as closely as possible. JASP asks for a t value, and two sample sizes. This is slightly different to how you’d do it in bayesplay where you’d enter a d value and two sample sizes.

For alternative the prior, JASP gives you the choice between a normal distribution, a cauchy distribution, and a student t distribution. Since this is just a quick proof of concept, I have decided against implementing all three. And because a cauchy distribution is just a student t distribution where the degrees of freedom is equal to 1 I’ve decided to just implement the student t prior. For this prior you need a location, a scale and a df value.

For the robustness analysis itself, I had to decide which values you can vary. Varying all three at once (location, scale, and df) get’s tricky because then you have to visualise that in three dimensions. So instead, I’ve made it possible to vary the location and the scale together, the location by itself, and the scale by itself. Hopefully everything down below is self explanatory. So on to the calculator and plots themselves! I’ve already pre-filled it with the data from Gronau et al. (2020).

The calculator

The data

viewof tvalue = Inputs.range([-10, 10], {label: md`*t*`, value: -.90, step:0.01})

viewof n1 = Inputs.range([1, 100], {label: md`*n1*`, step: 1, value: 53})

viewof n2 = Inputs.range([1, 100], {label: md`*n2*`, step: 1, value: 57})

Alternative hypothesis

viewof tail = Inputs.radio(new Map([ ["Group 1 > Group 2", 1], ["Group 1 < Group 2", -1], ["Group 1 ≠ Group 2", 0]]), 
{value: 1, label: html`Alt. Hypothesis`, format: ([x, _]) => html`<span class="test" ${x}; margin-bottom: -2px;">${x}`})

Bayes factor

viewof bf_type = Inputs.radio(["BF10", "BF01", "Log(BF10)"], 
{value: "BF01", label: html`Bayes factor`, format: x => html`<span ${x}; margin-bottom: -2px;">${x}`})

The prior

viewof location = Inputs.range([-10, 10], {label: md`*location*`, step: 0.001, value: 0.350})
viewof scale = Inputs.range([0.001, 100], {label: md`*scale*`, step: 0.001, value: 0.103})
viewof df = Inputs.range([1, 30], {label: md`*df*`, step: 1, value: 3})

bf = {
  return wasm.bf_twosample_t(tvalue, n1, n2, location, scale, df, tail)
}

{
  if (bf != "") {
  let bf_fmt = format_bf(bf, bf_type)
  let p_value = Math.round(wasm.ttest(tvalue, n1+n2, tail) * 1000) / 1000;
  let bf_res = bf_fmt > 1 ? Math.round(bf_fmt * 10)/ 10 : Math.round(bf_fmt * 100) / 100
  return html`
  <style type="text/css">
.tg .tg-0lax{text-align:center;vertical-align:top}
.tg .tg-0pky{border-color:inherit;text-align:center;vertical-align:top}
</style>
<p style="margin-top:2em;">Results</p>
<table class="tg" width=500px>
<thead>
  <tr style="border-top-style:solid;border-top-width:medium;border-bottom-width:thin;border-color:black;border-right-style:none;border-left-style:none;">
    <th class="tg-0lax" style="text-transform:none;"><i>t</i></th>
    <th class="tg-0lax" style="text-transform:none;">n1</th>
    <th class="tg-0lax" style="text-transform:none;">n2</th>
    <th class="tg-0lax" style="text-transform:none;">${bf_type}</th>
    <th class="tg-0lax" style="text-transform:none;"><i>p</i></th>
  </tr>
</thead>
<tbody>
  <tr style="border-bottom-style:solid;border-bottom-width:medium;border-bottom-color:black">
    <td class="tg-0pky">${tvalue}</td>
    <td class="tg-0pky">${n1}</td>
    <td class="tg-0pky">${n2}</td>
    <td class="tg-0pky">${bf_res}</td>
    <td class="tg-0pky">${p_value}</td>
  </tr>
</tbody>
</table>
  `
}
return md`Click compute to generate results`
}

Robustness Check

Location

viewof location_min = Inputs.range([-10, 10], {label: md`*location minimum*`, step: 0.001, value: 0.350 - 1})
viewof location_max = Inputs.range([-10, 10], {label: md`*location maximum*`, step: 0.001, value: 0.350 + 2})
viewof location_step = Inputs.range([0, 10], {label: md`*location precision*`, step: 0.001, value: 0.05})

viewof locations = md`Location steps: ${Math.round((location_max - location_min) / location_step)}`

Scale

viewof scale_min = Inputs.range([0, 10], {label: md`*scale minimum*`, step: 0.001, value: 0.103 - 1})
viewof scale_max = Inputs.range([0, 10], {label: md`*scale maximum*`, step: 0.001, value: 0.103 + 2})
viewof scale_step = Inputs.range([0, 10], {label: md`*scale precision*`, step: 0.001, value: 0.05})

viewof scales = md`Scale steps: ${Math.round((scale_max - scale_min) / scale_step)}`

viewof total_steps = md`Total Bfs to compute: ${Math.round((scale_max - scale_min) / scale_step) *  Math.round((location_max - location_min) / location_step)  }`

Robustness plot

mutable robust = {}
viewof computerobut = Inputs.button("Compute", { reduce: () => {
  var start_time = performance.now()
  let rr = wasm.bf_twosample_t_robust(tvalue, n1, n2, location, scale, df, tail, location_min, location_max, location_step, scale_min, scale_max, scale_step)
  let rr_location = wasm.bf_twosample_t_robust_location(tvalue, n1, n2, location, scale, df, tail, location_min, location_max, location_step, scale_min, scale_max, scale_step)
  let rr_scale = wasm.bf_twosample_t_robust_scale(tvalue, n1, n2, location, scale, df, tail, location_min, location_max, location_step, scale_min, scale_max, scale_step)
  let end_time = performance.now()
  let compute_time = (end_time - start_time)
  let data = {rr: rr, rr_scale: rr_scale, rr_location: rr_location, compute: compute_time}
  window.data = data
  mutable robust = data
  } 
  })

md`*compute time: ${Math.round(robust.compute)}ms*`

RectPlot = (data, xlab, ylab, xstep, ystep, obx, oby, thres) => {
  const xmin = Math.min(...data.map((x) => x[xlab]));
  const xmax = Math.max(...data.map((x) => x[xlab]));
  const ymin = Math.min(...data.map((x) => x[ylab]));
  const ymax = Math.max(...data.map((x) => x[ylab]));

  const thresher = (bf, thres) => {
    if (bf >= thres) {
      return "H1";
    }
    if (bf <= (1 / thres)) {
      return "H0";
    }
    return "Inc";
  };

 data.map((x) => (x.hyp = thresher(x.bf, thres)));
  const grid = vl
    .markRect({ clip: false })
    .data(data)
    .encode(
      vl
        .x({ bin: { maxbins: (xmax - xmin) / xstep } })
        .fieldQ(xlab)
        .axis({
          values: generate_range(xmin, xmax, 5),
          grid: false
        })
        .title(xlab),
      vl
        .y({ bin: { maxbins: (ymax - ymin) / ystep } })
        .fieldQ(ylab)
        .axis({
          values: generate_range(ymin, ymax, 5),
          grid: false
        })
        .title(ylab),
      vl
        .color()
        .fieldO("hyp")
        .title("Hypothesis")
        .scale({
          domain: ["H1", "H1", "Inc"],
          range: ["#bd6855", "green", "orange"],
          ticks: ["a", "b", "c"]
        }),
      vl
        .fill()
        .fieldO("hyp")
        .title("")
        .scale({
          domain: ["H0", "H1", "Inc"],
          range: ["#bd6855", "green", "orange"]
        })
        .legend({
          labelExpr:
            "{'H1': 'Evidence for H1', 'Inc': 'No evidence','H0': 'Evidence for H0'}[datum.label]"
        })
    );


 const point = vl
    .markPoint({ color: "black", fill: "black" })
    .data([{ x: obx, y: oby }])
    .encode(vl.x().fieldQ("x"), vl.y().fieldQ("y"));

 return vl.layer(grid, point).title(`Regions marked according to BF > ${thres}`).width(400).height(400);

}

viewof thres = Inputs.range([1, 20], {value: 10, step: 1, label: "Thershold"})
RectPlot(robust.rr, "scale", "location", scale_step, location_step, scale, location, thres).render()

ScalePlot = (data) => {
  
  data.map((x) => x.bf = format_bf(x.bf, bf_type))
  const line = vl
    .markLine({color: "black", clip: true})
    .data(data)
    .encode(
      vl
        .x()
        .fieldQ("scale")
        .title("scale"),
      vl.y().fieldQ("bf").title(`Bayes factor [${bf_type}]`)
        
    ).title(`Scale robustness plot for location = ${location}`).width(510).height(300).render()

    return line;

  }
ScalePlot(robust.rr_location)

LocationPlot = (data) => {
  data.map((x) => x.bf = format_bf(x.bf, bf_type))
  const line = vl
    .markLine({color: "black", clip: true})
    .data(data)
    .encode(
      vl
        .x()
        .fieldQ("location")
        .title("location"),
      vl.y().fieldQ("bf").title(`Bayes factor [${bf_type}]`)
        
    ).title(`Location robustness plot for scale = ${scale}`).width(510).height(300).render()

    return line;

  }
LocationPlot(robust.rr_scale)

format_bf = (bf, format) => {
    if (format == "BF10") {
        return 1 / bf
      }

    if (format == "BF01") {
        return bf
      }

    return  Math.log10(1/bf)
  }

generate_range = (min, max, steps) => {
  var step_size = (max - min) / steps
  if (step_size !== 0) {
    const steps = Math.abs(min - max) / step_size;
    const range = Array(steps + 1)
      .fill(min)
      .map((e, i) => e + i * step_size);
    return range;
  }
  return [min];
}

Final thoughts

First and foremost, I just can’t believe how fast Rust is. I knew Rust was fast, but I didn’t expect there to be such a huge difference between my Rust implementation and the equivalent Go implementation. Since the code is very simple I didn’t run into any of the usual headaches with the borrow checker that can sometimes make writing Rust difficult. The only real Rust gotcha was that two closures in Rust are always different types. So you can’t have a function that returns one of two possible closures as you might be able to do in many other languages. This not really a problem. Just something I found surprising.

The other thought is that I really wish there were some great numerics libraries in Rust. Go has Gonum, which is really good. In Rust, there’s the option of peroxide. However, interestingly, both these libraries really lacked something that was crucial for my task here. That is, they lack an implementation of the non-central t distribution. With Rust, there are some libraries that wrap C/C++ libraries, but getting these to compile to a WASM target didn’t work out the box, so I just went and implemented it myself. Thankfully I already have a Go implementation based on SciPy’s implementation (with the addition of a normal approximation when it fails). So I could just translate this. To throw a slight spanner in the works I also had to implement the type-(p,q) generalised hypergeometric function, because the non-central t implementation relies on it. All this just made me wish that there were great libraries like the GNU Scientific Library, the Math.NET Numerics library and the Boost numeric libraries available in Rust (or Go).

References

Gronau, Q. F., Ly, A., & Wagenmakers, E.-J. (2020). Informed bayesian t-tests. The American Statistician, 74(2), 137–143. https://doi.org/10.1080/00031305.2018.1562983

Reuse

https://creativecommons.org/licenses/by/4.0/

Citation

BibTeX citation:

@online{jcolling2022,
  author = {Lincoln J Colling},
  editor = {},
  title = {Robustness Analysis for Informed {Bayesian} t-Tests},
  date = {2022-04-14},
  langid = {en},
  abstract = {A Rust library for informated Bayesian t-tests.}
}

For attribution, please cite this work as:

Lincoln J Colling. (2022, April 14). Robustness analysis for informed Bayesian t-tests.