# Bayes Testing

**URL:** <https://discourse.julialang.org/t/bayes-testing/60373>\
**Category:** Statistics\
**Tags:** question, package\
**Created:** [May 1, 2021, 10:16am UTC](https://discourse.julialang.org/t/bayes-testing/60373 "2021-05-01T10:16:34Z")\
**Posts on this page:** 8\
**Page:** 2

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [May 3, 2021, 7:36am UTC](https://discourse.julialang.org/t/bayes-testing/60373/21 "2021-05-03T07:36:01Z")

</div>

Here is a version of the JZS t-test, I need to wirte out the fromula of Student’s t-test, or find a pacakge, that only expresses the t-value, and not a full description.

[http://pcl.missouri.edu/bf-one-sample](http://pcl.missouri.edu/bf-one-sample)  
[http://pcl.missouri.edu/sites/default/files/Rouder.bf\_.pdf](http://pcl.missouri.edu/sites/default/files/Rouder.bf_.pdf)

```julia
mu=mean(group1) # the group that is not manipulated
t=(mean(x)*null_hypothesis)/(std(x)/sqrt(N))
N=length(x) 
B=(1+(t^2/v)^inv(N/2))/
quadgk(sqrt(1+N*g)*((1+(t^2/((1+N*g)*v)^(v*1/2))*
inv(sqrt(2*pi*g))*g^(-3/2)*exp(-1/(2*g)))),0,inf)	

```

but I don’t know what g means

---

<div class="post-metadata">

**Author:** ![Storopoli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/storopoli/32/209278_2.png) [@Storopoli](https://discourse.julialang.org/u/Storopoli)\
**Post date:** [May 3, 2021, 9:14am UTC](https://discourse.julialang.org/t/bayes-testing/60373/22 "2021-05-03T09:14:56Z")

</div>

Bayesian t-test is more easy, you can try to convert a Stan implementation of John Kruschke’s Bayesian Estimation Supersedes the t-test (BEST), in John K. Kruschke, “Bayesian Estimation Supersedes the t test,” _Journal of Experimental Psychology_ 142, no. 2 (May 2013): 573-603, doi:10.1037/a0029146.

This should be easy to do in Turing

```stan
// Adapted from code by Michael Clark
// https://github.com/m-clark/Miscellaneous-R-Code/blob/master/ModelFitting/Bayesian/rstant_testBEST.R

// Stuff coming in from R
data {
  int<lower=1> N; // Sample size
  int<lower=2> n_groups; // Number of groups
  vector[N] y; // Outcome variable
  int<lower=1, upper=n_groups> group_id[N]; // Group variable
}

// Stuff to transform in Stan
transformed data {
  real mean_y;
  mean_y = mean(y); 
}

// Stuff to estimate
parameters {
  vector[2] mu; // Estimated group means 
  vector<lower=0>[2] sigma; // Estimated group sd
  real<lower=0, upper=100> nu; // df for t distribution
}

// Models and distributions
model {
  // Priors
  // curve(expr = dnorm(mean_y, 2), from = -5, to = 5)
  mu ~ normal(mean_y, 2);
  
  // curve(expr = dcauchy(x, location = 0, scale = 1), from = 0, to = 40)
  sigma ~ cauchy(0, 1);
  
  // Kruschke uses a nu of exponential(1/29)
  // curve(expr = dexp(x, 1/29), from = 0, to = 200)
  nu ~ exponential(1.0/29);
  
  
  // Likelihood
  for (n in 1:N){
    y[n] ~ student_t(nu, mu[group_id[n]], sigma[group_id[n]]);
  }
}

// Stuff to calculate with Stan
generated quantities {
  // Mean difference
  real mu_diff;
  
  // Effect size; see footnote 1 in Kruschke:2013
  // Standardized difference between two means
  // See https://en.wikipedia.org/wiki/Effect_size#Cohen's_d
  real cohen_d;
  
  // Common language effect size
  // The probability that a score sampled at random from one distribution will 
  // be greater than a score sampled from some other distribution
  // See https://janhove.github.io/reporting/2016/11/16/common-language-effect-sizes
  real cles;

  mu_diff = mu[1] - mu[2];
  cohen_d = mu_diff / sqrt(sum(sigma)/2);
  cles = normal_cdf(mu_diff / sqrt(sum(sigma)), 0, 1);
}

```

Stan Code from [https://www.andrewheiss.com/blog/2019/01/29/diff-means-half-dozen-ways/imdb\_best.stan](https://www.andrewheiss.com/blog/2019/01/29/diff-means-half-dozen-ways/imdb_best.stan).

One more thing? Why not fit two models that represents the competing hypotheses and do a LOO-CV using Pareto Smoothing Importance Sampling (PSIS-LOO)? The LOO calculates the Leave-one-out Cross-Validation of the elpd (_expected log pointwise predictive density_):

![Screen Shot 2021-05-03 at 06.10.35](https://global.discourse-cdn.com/julialang/original/3X/2/0/20cc5c8ad95ae636fc01eece8c1127bda385df55.png)

and see which model performs better? You wouldn’t even be restricted to only comparing two competing hypotheses at once. You could compare 10, 20, … all together.

It’s hard to convert frequentist approximations shortcuts (as E.T. Jaynes says “adhockeries”) to Bayesian Equivalent. As Allen Downey says ([Bayesian and frequentist results are not the same, ever – Probably Overthinking It](https://www.allendowney.com/blog/2021/04/25/bayesian-and-frequentist-results-are-not-the-same-ever/)):

> The posterior distribution represents everything you know about the parameters; if you reduce it to a single number, an interval, or a probability, you lose useful information. In fact, you lose exactly the information that makes the posterior distribution useful in the first place.

> It’s like comparing a car and an airplane by driving the airplane on the road. You would conclude that the airplane is complicated, expensive, and not particularly good as a car. But that would be a silly conclusion because it’s a silly comparison. The whole point of an airplane is that it can fly.

---

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [May 3, 2021, 4:03pm UTC](https://discourse.julialang.org/t/bayes-testing/60373/23 "2021-05-03T16:03:50Z")

</div>

Is this in r ? I did see a BEST algorithm written in Python

---

<div class="post-metadata">

**Author:** ![Storopoli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/storopoli/32/209278_2.png) [@Storopoli](https://discourse.julialang.org/u/Storopoli)\
**Post date:** [May 3, 2021, 5:21pm UTC](https://discourse.julialang.org/t/bayes-testing/60373/24 "2021-05-03T17:21:29Z")

</div>

Stan: [https://mc-stan.org/](https://mc-stan.org/)

---

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [May 3, 2021, 5:27pm UTC](https://discourse.julialang.org/t/bayes-testing/60373/25 "2021-05-03T17:27:27Z")

</div>

Here is BEST in Python:  
[https://github.com/strawlab/best/blob/master/best/\_\_init\_\_.py](https://github.com/strawlab/best/blob/master/best/ __init__.py)

I think I can translate most of this pretty easily, but I need to go over Gaussian KDE, I did find KernelDensity.jl

Also, is NonCentralT, some kind of t-test?

---

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [May 3, 2021, 5:28pm UTC](https://discourse.julialang.org/t/bayes-testing/60373/26 "2021-05-03T17:28:40Z")

</div>

OK, there’s a stan package for Julia, but it still involves downloading STAN. I think the python one is easier to change. I’ll get as far as I can, and message it to you.

---

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [May 5, 2021, 5:00am UTC](https://discourse.julialang.org/t/bayes-testing/60373/27 "2021-05-05T05:00:15Z")

</div>

Is “cles” the final Bayesian Power (term used in Krushke, 2017)?

---

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [May 12, 2021, 3:05am UTC](https://discourse.julialang.org/t/bayes-testing/60373/28 "2021-05-12T03:05:17Z")

</div>

Is this a type of Gaussian Mixed Model around a Non-central T distribution?

[Previous page](https://discourse.julialang.org/t/bayes-testing/60373.md?page=1)