rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Smiles and Leniency
# Greg Francis
# PSY 646
# 14 September 2018

# Bayesian t-tests

# It can take a few minutes for the code to get moving. During this time, the program is rewriting the model into another programming language and compiling it. Once it starts running, you will see output on the console window.

# load data file
PWdata<-read.csv(file="PhysiciansWeight.csv",header=TRUE,stringsAsFactors=TRUE)

# load the brms library
library(brms)

# By default brms provides priors that are adapted to fit within the range of your data. Generally, this produces results similar to standard linear regression. Later, we will introduce informative priors

model1 = brm(Time ~ Weight, data = PWdata, iter = 2000, warmup = 200, chains = 3, thin = 2 )

# print out summary of model
print(summary(model1))

# Check on model convergence and posteriors
plot(model1)

# Plot model fit to data
dev.new()
plot(conditional_effects(model1), points = TRUE)

# Get posterior distributions

post<-posterior_samples(model1)

# check probability that difference is less than zero
dev.new()
plot(density(post$b_Weightoverweight))
cat("Probability that group mean difference is greater than zero is: ", sum(post$b_Weightoverweight > 0) / length(post$b_Weightoverweight))

# Compare to traditional t-test
traditional <- t.test(Time ~ Weight, data=PWdata, var.equal=TRUE)

print(traditional)

# Compare to traditional linear regression
print(summary(lm(Time~Weight, data=PWdata)))

# Compare to a model with only one mean
model2 = brm(Time ~ 1, data = PWdata, iter = 2000, warmup = 200, chains = 3, thin = 2 )

print(LOO(model1, model2) )

# --------------------------------

# Bayesian Welch's test (unequal variances test)

# Define a formula that includes variations in sigma by group
uneq_var_frm <- brmsformula(Time ~ Weight, sigma ~ Weight)

model3 = brm(uneq_var_frm, data = PWdata, iter = 2000, warmup = 200, chains = 3, thin = 2 )

# print out summary of model
print(summary(model3))

# Check on model convergence and posteriors
dev.new()
plot(model3)

# Plot model fit to data
dev.new()
plot(conditional_effects(model3), points = TRUE)

# Get posteriors

post<-posterior_samples(model3)

# check probability that difference is less than zero
cat("Probability that group mean difference is greater than zero is: ", sum(post$b_Weightoverweight > 0) / length(post$b_Weightoverweight), "\n")

# Compare to traditional Welch's test
traditional <- t.test(Time ~ Weight, data=PWdata, var.equal=FALSE)
print(traditional)

# -----------------
# Compare one sigma model against two sigmas model

print("Leave One Out Cross Validation (smaller is better)")
print(LOO(model1,model3))

# -----------------
# Compare models against each other

print("Leave One Out Cross Validation (smaller is better)")
print(LOO(model1, model2, model3))


#------------------
# Set non-default priors on means

model5 = brm(Time ~ Weight, data = PWdata, iter = 2000, warmup = 200, chains = 3, thin = 2, prior = c(prior(normal(0, 10), class = "Intercept"), prior(normal(0, 10), class = "b"))  )
print(summary(model5))


# Set a kind of crazy prior on the first mean

model6 = brm(Time ~ Weight, data = PWdata, iter = 2000, warmup = 200, chains = 3, thin = 2, prior = c(prior(normal(10, 1), class = "Intercept"), prior(normal(0, 10), class = "b"))  )
print(summary(model6))