rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Smiles and Leniency
# Greg Francis
# PSY 626
# 16 September 2020

# simple ANOVA

# load the rethinking library
library(rethinking)

# load data file
SLdata<-read.csv(file="SmilesLeniency.csv",header=TRUE,stringsAsFactors=TRUE)

# Null model: one mean for all smile types
SLmodel0 <- map(
			alist( Leniency ~ dnorm(mu, sigma), 
			mu <- a,
			a~ dnorm(5, 5),
			sigma ~ dunif(0, 10)
			), data= SLdata ) 
cat("Finished SLmodel0\n")

print(SLmodel0)


# Dummy variable approach
SLdata$False <- ifelse(SLdata$SmileType =="FALSE", 1, 0)  
SLdata$Felt <- ifelse(SLdata$SmileType =="Felt", 1, 0)  
SLdata$Miserable <- ifelse(SLdata$SmileType =="Miserable", 1, 0)  
SLdata$Neutral <- ifelse(SLdata$SmileType =="Neutral", 1, 0)  


# Alternative model 1: different mean for each condition (dumy variable coding)	
SLmodel1 <- map(
			alist( Leniency ~ dnorm(mu, sigma), 
			mu <- a1*False + a2*Felt + a3*Miserable + a4*Neutral,
			a1 ~ dnorm(5, 5),
			a2 ~ dnorm(5, 5),
			a3 ~ dnorm(5, 5),
			a4 ~ dnorm(5, 5),
			sigma ~ dunif(0, 10)
			), data= SLdata  ) 
cat("Finished SLmodel1\n")

plot(coeftab(SLmodel0, SLmodel1))

# Examine posterior distributions
post<-extract.samples(SLmodel1, n= 2000)

# Plot posterior for mean of felt smile type
dev.new()
dens(post$a2, col=rangi2, lwd=2, xlab="Value", main="posterior for a2", show.HPDI=0.89)


# Probability that mean for neutral smile type is less than mean for Felt smile type?

cat("Probability neutral mean is smaller than felt mean:\n")
print( length(post$a4[post$a4 <= post$a2])/length(post$a4) )

cat("Probability neutral mean is smaller than false mean:\n")
print( length(post$a4[post$a4 <= post$a1])/length(post$a4) )

cat("Probability neutral mean is smaller than miserable mean:\n")
print( length(post$a4[post$a4 <= post$a3])/length(post$a4) )

cat("Probability miserable mean is smaller than felt mean:\n")
print( length(post$a3[post$a3 <= post$a2])/length(post$a2) )

cat("Probability felt mean is smaller than false mean:\n")
print( length(post$a2[post$a2 <= post$a1])/length(post$a2) )


# Contrast model: neutral versus any smile type
SLmodel2 <- map(
			alist( Leniency ~ dnorm(mu, sigma), 
			mu <- a1*(1-Neutral) + a4*Neutral,
			a1 ~ dnorm(5, 5),
			a4 ~ dnorm(5, 5),
			sigma ~ dunif(0, 10)
			), data= SLdata  ) 
cat("Finished SLmodel2\n")

print(compare(SLmodel0, SLmodel1, SLmodel2))

print(compare(SLmodel0, SLmodel2))

# Another way to set up the alternative model
# Dummy variable to indicate condition
SLdata$ConditionIndex <-0*SLdata$Leniency # just makes an array of the right length, filled with 0's
SLdata$ConditionIndex[SLdata$SmileType =="FALSE"] =1
SLdata$ConditionIndex[SLdata$SmileType =="Felt"] =2
SLdata$ConditionIndex[SLdata$SmileType =="Miserable"] =3
SLdata$ConditionIndex[SLdata$SmileType =="Neutral"] =4


# Alternative model 1: different mean for each condition (dumy variable coding)	
SLmodel3 <- map(
			alist( Leniency ~ dnorm(mu, sigma), 
			mu <- a[ConditionIndex],
			a[ConditionIndex] ~ dnorm(5, 5),
			sigma ~ dunif(0, 10)
			), data= SLdata  ) 
cat("Finished SLmodel3\n")

print(SLmodel3)

# Show how to extract posterior from this model set up
post<-extract.samples(SLmodel3, n= 2000)

cat("Probability neutral mean is smaller than felt mean:\n")
print( length(post$a[,4][post$a[,4] <= post$a[,2]])/length(post$a[,4]) )