rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Visual Search
# Greg Francis
# PSY 646
# 30 August 2018

# fit a linear model that predicts response time as a function of the number of distractors

# 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 full data file
VSdata<-read.csv(file="VisualSearch.csv",header=TRUE,stringsAsFactors=FALSE)

# pull out just the trials for the first participant's conjunctive target present condition
VSdata2<-subset(VSdata, VSdata$Participant=="Francis200S16-2" & VSdata$Target=="Absent" & VSdata$DistractorType=="Conjunction")
#  
# 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(RT_ms ~ NumberDistractors, data = VSdata2, iter = 2000, warmup = 200, chains = 3,
    thin = 2 )

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

# Check on model convergence and posteriors
plot(model1)

# Extract posteriors from samples
post<-posterior_samples(model1)

# Compute posterior probability that b_NumDistractors is less than 42
probLess42 <-length(post$b_NumberDistractors[post$b_NumberDistractors < 42])/length(post$b_NumberDistractors) 
cat("Probability b_NumberDistractors is less than 42 = ", probLess42, "\n")

# Compute posterior probability for RT_ms for NumberDistractors = 35
RT_at_35 <- post$b_Intercept + post$b_NumberDistractors * 35
dev.new()
plot(density(RT_at35))

# Compute posterior probability that RT_ms for 35 distractors is greater than 2400
probGT2400 <-length(RT_at_35[RT_at_35 > 2400])/length(RT_at_35) 
cat("Probability RT_ms for 35 distractors is greater than 2400 = ", probGT2400, "\n")

# set up a new graphics window, so as to not write over previous one
dev.new()
# Plot fitting line and HPDI
plot(marginal_effects(model1), points = TRUE)




# Try a different prior on the slope parameter N(30,1). Indicates strong belief that the slope will be around 30, with little deviation.

model2 = brm(RT_ms ~ NumberDistractors, data = VSdata2, iter = 2000, warmup = 200, chains = 3,
    thin = 2, prior = c( prior(normal(30, 1), class = "b")) )

print(summary(model2))
dev.new()
plot(marginal_effects(model2), points = TRUE)
dev.new()
plot(model2)


# Try a different prior on the slope parameter N(55,1). Indicates strong belief that the slope will be around 55, with little deviation.

model3 = brm(RT_ms ~ NumberDistractors, data = VSdata2, iter = 2000, warmup = 200, chains = 3,
    thin = 2, prior = c( prior(normal(55, 1), class = "b")) )

print(summary(model3))
dev.new()
plot(marginal_effects(model3), points = TRUE)
dev.new()
plot(model3)

# Try a different prior on the slope parameter N(55,10). Indicates belief that the slope will be around 55, with some deviation.

model4 = brm(RT_ms ~ NumberDistractors, data = VSdata2, iter = 2000, warmup = 200, chains = 3,
    thin = 2, prior = c( prior(normal(55, 10), class = "b")) )

print(summary(model4))
dev.new()
plot(marginal_effects(model4), points = TRUE)
dev.new()
plot(model4)