---
title: "Stat. 481: Final - BayesTestR"
output: html_notebook
---

This is some of the code from the [bayestestR](https://easystats.github.io/bayestestR/index.html) website.  This is an excellent package for testing hypothesis when applying Bayesian methods.  And for trying out some [Stan](https://mc-stan.org/) code, a more modern software for Bayesian modeling.

Read the [Example1](https://easystats.github.io/bayestestR/articles/example1.html) and [Example2](https://easystats.github.io/bayestestR/articles/example2.html) from the [bayestestR](https://easystats.github.io/bayestestR/index.html) website.

# Example 1

Get the following code to run.

```{r}
library(rstanarm)
library(bayestestR)
library(insight)
```

This code fits a linear regression using Classical least squares for the two variables from the iris dataset.

**Question 1:** Replace the variables in the code below for fitting the regression model to Sepal.Width and Sepal.Length.  That is replace y = Sepal.Length with Sepal.Width and x = Petal.Length with Sepal.Length.

Answer the question below.

```{r}
iris
```


```{r}
model <- lm(Sepal.Length ~ Petal.Length, data=iris)
summary(model)
```

```{r}
library(ggplot2)  # Load the package

# The ggplot function takes the data as argument, and then the variables 
# related to aesthetic features such as the x and y axes.
ggplot(iris, aes(x=Petal.Length, y=Sepal.Length)) +
  geom_point() +  # This adds the points
  geom_smooth(method="lm") # This adds a regression line
```

Now lets try the same model using a Bayesian model.  The function *stan_glm* is a function that uses the Stan language to estimate the Bayesian Regression model using flat nomal priors on the slope and intercept. 

```{r}
model <- stan_glm(Sepal.Length ~ Petal.Length, data=iris)
describe_posterior(model)

prior_summary(model)
plot(model, "hist")
```

```{r}
posteriors <- insight::get_parameters(model)

head(posteriors)  # Show the first 6 rows
```

```{r}
ggplot(posteriors, aes(x = Petal.Length)) +
  geom_density(fill = "orange")

mean(posteriors$Petal.Length)
median(posteriors$Petal.Length)
map_estimate(posteriors$Petal.Length)
```

```{r}
ggplot(posteriors, aes(x = Petal.Length)) +
  geom_density(fill = "orange") +
  # The mean in blue
  geom_vline(xintercept=mean(posteriors$Petal.Length), color="blue", size=1) +
  # The median in red
  geom_vline(xintercept=median(posteriors$Petal.Length), color="red", size=1) +
  # The MAP in purple
  geom_vline(xintercept=map_estimate(posteriors$Petal.Length), color="purple", size=1)
```

```{r}
hdi(posteriors$Petal.Length, ci=0.89)
```


```{r}
describe_posterior(model, test = c("p_direction", "bayesfactor"))
```


**Question:** Is there evidence that the slope and intercept are different from 0?  Recall that if the Bayes Factor is bigger than 3 or less than 1/3 then there is evidence that the slope and intercept are different from zero.

**Answer:** <<< Type your answer here. >>>

# Example 2

Correlations

**Question 2:** Run the code below to test if there is a correlation between Sepal.Width and Sepal.Length.

```{r}
result <- cor.test(iris$Sepal.Width, iris$Sepal.Length)
result
```

```{r}
iris %>% ggplot( aes( y = Sepal.Width, x = Sepal.Length)) +
  geom_point()
```

Using the BayesFactor R package.

```{r}
library(BayesFactor)
result <- correlationBF(iris$Sepal.Width, iris$Sepal.Length)
```

```{r}
describe_posterior(result)
```

```{r}
bayesfactor(result)
```

```{r}
library(see)

plot(bayesfactor(result)) +
  scale_fill_pizza()
```

**Question:** Is there evidence that the correlation is different from 0?  Recall that if the Bayes Factor is bigger than 3 or less than 1/3 then there is evidence that the correlation is different from zero.

**Answer:** <<< Type your answer here. >>>


# T-tests

**Question 3:**  Run the following code for each Species in the iris dataset.

**Question:** Is there a difference in Sepal.Width for each pair of Species.

**Answer:** <<< Type your answer here. >>>

```{r}
iris %>% distinct(Species)
```


```{r}
library(dplyr)
library(ggplot2)

# Select only two relevant species
data <- iris %>% 
  filter(Species != "setosa") %>% 
  droplevels()

# Visualise distributions and observations
data %>% 
  ggplot(aes(x = Species, y = Sepal.Width, fill = Species)) +
  geom_violindot(fill_dots = "black", size_dots = 1) +
  scale_fill_material() +
  theme_modern()

```

```{r}
result <- BayesFactor::ttestBF(formula = Sepal.Width ~ Species, data = data)
describe_posterior(result)
```


```{r}
library(see)

plot(bayesfactor(result)) +
  scale_fill_pizza()
```