rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Masking Laws
# Greg Francis
# PSY 646
# 16 February 2024

# load full data file
MLdata<-read.csv(file="MaskingLaws.csv",header=TRUE,stringsAsFactors=FALSE)

# load the rethinking library
library(rethinking)
			

# STA law
MLcurveSTA <- ulam(
			alist( PercentCorrect ~ dnorm(mu, sigma), 
			mu <- a[TMDurationsCondition] +b[TMDurationsCondition]*(SOA-(bottom-MaskDuration+TargetDuration))*(SOA-(bottom-MaskDuration+TargetDuration)), 
			a[TMDurationsCondition] ~  dnorm(50, 50),
			b[TMDurationsCondition] ~ dnorm(0, 1),
			bottom ~ dnorm(125, 500),
			sigma ~ dunif(0, 500)
			), data= MLdata, iter=10000, warmup=4000, log_lik=TRUE, chains=4 ) 	

# FM ISI law
MLcurveFMISI <- ulam(
			alist( PercentCorrect ~ dnorm(mu, sigma), 
			mu <- a[TMDurationsCondition] +b[TMDurationsCondition]*(SOA-(bottom[Mcondition]+TargetDuration))*(SOA-(bottom[Mcondition]+TargetDuration)), 
			a[TMDurationsCondition] ~  dnorm(50, 50),
			b[TMDurationsCondition] ~ dnorm(0, 1),
			bottom[Mcondition] ~ dnorm(125, 255),
			sigma ~ dunif(0, 500)
			), data= MLdata, log_lik=TRUE, chains=4 ) 	

# ED SOA law			
MLcurveEDSOA <- ulam(
			alist( PercentCorrect ~ dnorm(mu, sigma), 
			mu <- a[TMDurationsCondition] +b[TMDurationsCondition]*(SOA-(bottom[EDSOAConditions]))*(SOA-(bottom[EDSOAConditions])), 
			a[TMDurationsCondition] ~  dnorm(50, 50),
			b[TMDurationsCondition] ~ dnorm(0, 1),
			bottom[EDSOAConditions] ~ dnorm(125, 255),
			sigma ~ dunif(0, 500)
			), data= MLdata, iter=5000, warmup=2000, log_lik=TRUE, chains=4 ) 				

# No law
MLcurveNoLaw <- ulam(
			alist( PercentCorrect ~ dnorm(mu, sigma), 
			mu <- a[TMDurationsCondition] +b[TMDurationsCondition]*(SOA-bottom[TMDurationsCondition])*(SOA-bottom[TMDurationsCondition]), 
			a[TMDurationsCondition] ~  dnorm(50, 50),
			b[TMDurationsCondition] ~ dnorm(0, 1),
			bottom[TMDurationsCondition] ~ dnorm(125, 255),
			sigma ~ dunif(0, 500)
			), data= MLdata, iter=5000, warmup=2000, log_lik=TRUE, chains=4 ) 	

print(compare(MLcurveSTA, MLcurveFMISI, MLcurveEDSOA, MLcurveNoLaw))