Hello. I’m new to the “Bayesian” world and wanted to ask how hypothesis testing is done properly since p-values are frowned upon amongst Bayesians. I have heard about the Bayes Factor but didn’t see a practical example yet.

Let’s say, I have the model:

y = alpha + beta * x

And I want to know if predictor x has a significant outcome on y. How do I go on from here?

Also should I do this with an R function or directly code it into my Stan file?

Here’s a minimum viable code:

R part:

```
## Load Data
data("radon", package = "rstanarm")
data <- radon
# Seperate variables
logradon <- data$log_radon
floorind <- data$floor
N <- length(logradon)
# Lists for STAN
cplist <- list(Y = logradon,
X = floorind,
N = N)
```

Stan part

```
data {
int<lower=1> N;
vector[N] X;
vector[N] Y;
}
parameters {
real alpha;
real beta;
real<lower=0> sigma_y;
}
model {
vector[N] mu_y;
for (i in 1:N)
mu_y[i] = alpha + beta * X[i];
Y ~ normal(mu_y, sigma_y);
sigma_y ~ cauchy(0, 2.5);
alpha ~ normal(0, 5);
beta ~ normal(0, 5);
}
//vector[N] log_lik;
//for (i in 1:N) log_lik[i] = normal_lpdf(Y[i]| alpha + beta * X[i], sigma_y);
generated quantities {
vector[N] y_pred;
vector[N] indiv_squared_errors;
real <lower = 0> sum_of_squares;
real <lower = 0> root_mean_squared_error;
for (i in 1:N) {
y_pred[i] = alpha + beta * X[i];
indiv_squared_errors[i] = (Y[i] - y_pred[i])^2;
}
sum_of_squares = sum(indiv_squared_errors);
root_mean_squared_error = sqrt(sum_of_squares/N);
}
```