Simple Bayesian methods of linear regression and testing for significant differences between regression line slopes

bayes-regression-slopes.R

# Investigating some simple Bayesian methods of linear regression
# and testing for significant differences between regression line slopes

n <- 100
x <- seq(n)
y1 = 5 * x + 150
y2 = 1.5 * x + 50

d1 <- data.frame(x1 = x  + rnorm(n, sd = n/10),
                 y1 = y1)

d2 <- data.frame(x2 = x  + rnorm(n, sd = n/10),
                 y2 = y2)


library(ggplot2)
ggplot() +
  geom_point(aes(x1, y1), d1, col = "blue") +
  geom_point(aes(x2, y2), d2, col = "red") +
  theme_minimal()

library(MCMCpack)
posterior1 <- (MCMCregress(y1 ~ x1, d1))
# windows()
plot(posterior1)
raftery.diag(posterior1)
summary(posterior1)

# intercept
hist(posterior1[,1])
# x1 (ie. slope)
hist(posterior1[,2])
# get distribution of slope
x1_post <- as.numeric(posterior1[,2])

posterior2 <- (MCMCregress(y2 ~ x2, d2))
# windows()
plot(posterior2)
raftery.diag(posterior2)
summary(posterior2)

# intercept
hist(posterior2[,1])
# x1 (ie. slope)
hist(posterior2[,2])
# get distribution of slope
x2_post <- as.numeric(posterior2[,2])

# now compare the two distributions of slopes


# BEST http://www.indiana.edu/~kruschke/BEST/
# devtools::install_github('mikemeredith/BEST')
library(BEST)
slopes_test <- BESTmcmc(x1_post, x2_post, verbose = TRUE)
summary(slopes_test)
plot(slopes_test)
windows()
plotAll(slopes_test, credMass=0.8, ROPEm=c(-0.1,0.1), 
        ROPEeff=c(-0.2,0.2), compValm=0.5) 
pairs(slopes_test)

# Bayes Factors http://bayesfactorpcl.r-forge.r-project.org/
#devtools::install_github( username = "richarddmorey", repo = "BayesFactor", subdir='pkg/BayesFactor', dependencies = TRUE)
library(BayesFactor)
bf = ttestBF(x1_post, x2_post)
chains = posterior(bf, iterations = 10000)
windows()
plot(chains[, 1:4])
plot(chains[,"mu"])

# https://github.com/rasmusab/bayesian_first_aid
# devtools::install_github("rasmusab/bayesian_first_aid")
library(BayesianFirstAid)
bayes.t.test(x1_post, x2_post)


# http://madere.biol.mcgill.ca/cchivers/intro_bayesian_chivers.R
library(MHadaptive)

## A LINEAR REGRESSION EXAMPLE ####
## Define a Bayesian linear regression model
li_reg<-function(pars,data)
{
  a<-pars[1]      #intercept
  b<-pars[2]      #slope
  sd_e<-pars[3]   #error (residuals)
  if(sd_e<=0){return(NaN)}
  pred <- a + b * data[,1]
  log_likelihood<-sum( dnorm(data[,2],pred,sd_e, log=TRUE) )
  prior<- prior_reg(pars)
  return(log_likelihood + prior)
}

## Define the Prior distributions
prior_reg<-function(pars)
{
  a<-pars[1]          #intercept
  b<-pars[2]          #slope  
  epsilon<-pars[3]    #error
  
  prior_a<-dnorm(a,0,100,log=TRUE)     ## non-informative (flat) priors on all 
  prior_b<-dnorm(b,0,100,log=TRUE)     ## parameters.  
  prior_epsilon<-dgamma(epsilon,1,1/100,log=TRUE)      
  
  return(prior_a + prior_b + prior_epsilon)
}

mcmc_r<-Metro_Hastings(li_func=li_reg,pars=c(0,1,1),
                       par_names=c('a','b','epsilon'),data=d1)

mcmc_r<-Metro_Hastings(li_func=li_reg,pars=c(0,1,1),
                       prop_sigma=mcmc_r$prop_sigma,par_names=c('a','b','epsilon'),data=d1)
mcmc_r<-mcmc_thin(mcmc_r)
plotMH(mcmc_r)