rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Weber's Law
# Greg Francis
# PSY 646
# 01 April 2024

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

# filter parts we do not care about
WLdata2 <- subset(WLdata, WLdata$TrialType=="Experiment")

# Make a data frame with means and standard deviation for each combination of orientaton and wing type
# Compute the means and standard deviations for each condition (wings, orientation, length, participant)
WLmeans <- aggregate(WLdata2$MatchSize_pixels,
						list(WLdata2$SONA_ID, WLdata2$Wings, WLdata2$MatchOrientation, WLdata2$TargetLength), mean)
WLsds <- aggregate(WLdata2$MatchSize_pixels,
						list( WLdata2$SONA_ID, WLdata2$Wings, WLdata2$MatchOrientation, WLdata2$TargetLength), sd)
# Following computes means across all subjects, to test whether code is valid (it is)
WLmeans2 <- aggregate(WLdata2$MatchSize_pixels,
						list(WLdata2$Wings, WLdata2$MatchOrientation, WLdata2$TargetLength), mean)
WLsds2 <- aggregate(WLdata2$MatchSize_pixels,
						list( WLdata2$Wings, WLdata2$MatchOrientation, WLdata2$TargetLength), sd)						

plot(WLmeans2$x, WLsds2$x, main="Physical model")

dev.new()
plot(WLmeans2$Group.3, WLsds2$x, main="Perceptual model")


# Reorganize into a data frame with one row for each participant, wing, orientation, length, mean, sd
WLdata3 <- data.frame(WLmeans)
WLdata3$SD = WLsds$x		
colnames(WLdata3) <-c("Participant", "Wings", "MatchOrientation", "TargetLength", "MeanMatch", "SD")

# Set up indices for participants
SONAIDs<-unique(WLdata3$Participant)
WLdata3$ParticipantIndex = WLdata3$Participant
count=0
for(i in SONAIDs){
	count= count+1
	WLdata3$ParticipantIndex[WLdata3$Participant ==i] = count  
}

# Convert strings to numbers as appropriate
WLdata3$WingIsNone <- ifelse(WLdata3$Wings =="None", 1, 0)  		
WLdata3$WingIsOut <- ifelse(WLdata3$Wings =="Outward", 1, 0)  		
WLdata3$WingIsIn <- ifelse(WLdata3$Wings =="Inward", 1, 0)  		
WLdata3$MatchOrientationIndex <- ifelse(WLdata3$MatchOrientation =="Horizontal", 1, 0)  

																					
# load the rethinking library
library(rethinking)
	
# Full model (Wings, Orientation, Physical length) predicts SD 9(Not implemented)

WLdataFullPhysical <- data.frame(Participant=WLdata3$ParticipantIndex, WingIsNone =WLdata3$WingIsNone, WingIsOut =WLdata3$WingIsOut, WingIsIn =WLdata3$WingIsIn, MatchOrientationIndex =WLdata3$MatchOrientationIndex, MeanMatch =WLdata3$MeanMatch, SD =WLdata3$SD )		


# Physical length model
WLdataPhysical <- data.frame(Participant=WLdata3$ParticipantIndex,  MeanMatch =WLdata3$MeanMatch, SD =WLdata3$SD )		
											
PhysicalLengthLinearModel <- ulam(
			alist( SD ~ dnorm(mu, sigma), 
			mu <- a[Participant] + b[Participant]*MeanMatch,
			a[Participant] ~ dnorm(grand_mua, grand_sa),
			grand_mua ~ dnorm(5, 5),
			grand_sa ~ dcauchy(0, 5),
			b[Participant] ~ dnorm(grand_mub, grand_s),
			grand_mub ~ dnorm(5, 5),
			grand_s ~ dcauchy(0, 5),
			sigma ~ dcauchy(0, 5)
			), data= WLdataPhysical, iter=10000, warmup=4000, log_lik=TRUE, chains=4, cores=4 ) 
			
print("Finished physical model")					

										
# Perceived length model
WLdataPerceptual <- data.frame(Participant=WLdata3$ParticipantIndex, TargetLength =WLdata3$TargetLength, SD =WLdata3$SD )		
PerceptualMatchLinearModel <- ulam(
			alist( SD ~ dnorm(mu, sigma), 
			mu <- a[Participant] + b[Participant]* TargetLength,
			a[Participant] ~ dnorm(grand_mua, grand_sa),
			grand_mua ~ dnorm(5, 5),
			grand_sa ~ dcauchy(0, 5),
			b[Participant] ~ dnorm(grand_mub, grand_s),
			grand_mub ~ dnorm(5, 5),
			grand_s ~ dcauchy(0, 5),
			sigma ~ dcauchy(0, 5)
			), data= WLdataPerceptual, iter=10000, warmup=4000, log_lik=TRUE, chains=4, cores=4 ) 
			
print("Finished perceptual model")					
			
# Model comparison 
print(compare(PerceptualMatchLinearModel, PhysicalLengthLinearModel))