#########################################
#   Bayesian Regression Using blavaan   #
#              N Schalken               #
#           09 February 2017            #
#########################################

#------------------------------------------------------------------------------------------------------------------------------------------------------------------------
##############
# Exercise 1 #
##############

#1a. Download and load the package:
install.packages("blavaan")
require(blavaan)

#1b. Load data (make sure you set the working directory to the folder where your data is):

# To set your working directory to the right folder use:
# setwd("C:/your/folder/here")

#OR go to the Menu --> Session --> Set Working Directory --> Choose Directory

data <- read.table("regression.dat")

#1c. Give names to the variables in the dataset.
colnames(data)<- c("y", "x1", "x2")

#1d. Specify the regression model.
model.regression <- '#regression
                    y ~ x1 + x2

                    #variances
                    y ~~ y

                    #intercepts
                    y ~ 1'
#---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
##############
# Exercise 2 #
##############

#2a. Fitting the Bayesian regression model in blavaan. First specify the model object, and then specify the dataset.
fit.regression1 <- blavaan(model.regression, data=data)

#2b. Obtaining summary estimates of the parameters.
summary(fit.regression1)

#2c. Obtaining different plots of the parameters of the fitted model.
plot(fit.regression1, pars= 1:4, plot.type="trace")     #Traceplots 
plot(fit.regression1, pars= 1:4, plot.type="autocorr")  #Autocorrelation plots
plot(fit.regression1, pars= 1:4, plot.type="ecdf")      #ECDF plots
plot(fit.regression1, pars= 1:4, plot.type="density")   #Density plots

#---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
##############
# Exercise 3 #
##############

#3a. Fitting the same model in blavaan, but now with specified priors. In the previous model, default priors were automatically used.
fit.regression2 <- blavaan(model.regression, data=data, dp = dpriors(itheta = "dgamma(0.01, 0.01)", nu= "dnorm(0.0, 0.001)", beta = "dnorm(0.0, 0.001)"))

#3b. Summary estimates of the parameters in the new model.
summary(fit.regression2)

#Compare model 1 and model 2, are the results comparable? 

#2c. Obtaining different plots of the parameters of the second fitted model.
plot(fit.regression2, pars= 1:4, plot.type="trace")     #Traceplots 
plot(fit.regression2, pars= 1:4, plot.type="autocorr")  #Autocorrelation plots
plot(fit.regression2, pars= 1:4, plot.type="ecdf")      #ECDF plots
plot(fit.regression2, pars= 1:4, plot.type="density")   #Density plots


############################################################################################################################################################################