rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Visual Search
# Greg Francis
# PSY 646
# 11 September 2020

# 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 rethinking library
library(rethinking)

# fit a linear model that predicts response time as a function of the number of distractors
print("LINEAR MODEL")


VSmodel <- map(
			alist( RT_ms ~ dnorm(mu, sigma), 
			mu <- a + b*NumberDistractors,
			a ~ dnorm(500, 10),
			b ~ dnorm(0, 100),
			sigma ~ dunif(0, 500)
			), data=VSdata2, control=list(maxit=1000)  ) 
			
summaryTable <- summary(VSmodel)
print(VSmodel)		

# Another summary of the model estimates (corr does nothing here because the parameters are not correlated)
print(precis(VSmodel, corr=TRUE))	

plot(RT_ms ~ NumberDistractors, data=VSdata2)
abline(a=coef(VSmodel)["a"], b=coef(VSmodel)["b"], col=col.alpha("red",1.0))

numVariableLines=20
post<-extract.samples(VSmodel, n= numVariableLines)
for(i in 1: numVariableLines){
	abline(a=post$a[i], b=post$b[i], col=col.alpha("black",0.3))
}


# Frequentist linear regression

print("Frequentist linear regression")
lmModel <- lm(RT_ms ~ NumberDistractors, data=VSdata2)

print(summary(lmModel))




#plot posterior of a
numVariableLines=2000
post<-extract.samples(VSmodel, n= numVariableLines)

dev.new()
dens(post$a, col=rangi2, lwd=2, xlab="posterior for a")

# Probability that a > 900 ?
length(post$a[post$a >900])/length(post$a)
 
#plot posterior of b

dev.new()
dens(post$b, col=rangi2, lwd=2, xlab="posterior for b")

# Probability b is between 30 and 45
length(post$b[post$b > 30 & post$b<45])/length(post$b)


# posterior distribution for mean RT with 35 distractors
mu_at_35 <- post$a +post$b *35

dev.new()
dens(mu_at_35, col=rangi2, lwd=2, xlab="mu|NumDistract=35")

# Probability that RT_ms <2200 for 35 distractors?

 length(mu_at_35[mu_at_35 <= 2200])/length(mu_at_35)