rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Facial Feedback
# Greg Francis
# PSY 646
# 14 February 2024

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

# load the rethinking library
library(rethinking)

# Dummy variable to indicate condition
FFdata$ConditionIndex <-0*FFdata$Trial
FFdata$ConditionIndex[FFdata$Condition =="PenInTeeth"] =1
FFdata$ConditionIndex[FFdata$Condition =="NoPen"] =2
FFdata$ConditionIndex[FFdata$Condition =="PenInLips"] =3

# Other dummy variable approach
FFdata$PenInTeeth <- ifelse(FFdata$Condition =="PenInTeeth", 1, 0)  
FFdata$NoPen <- ifelse(FFdata$Condition =="NoPen", 1, 0)  
FFdata$PenInLips <- ifelse(FFdata$Condition =="PenInLips", 1, 0)  

FFdataClean<- data.frame(HappinessRating= FFdata$HappinessRating, Participant= FFdata$Participant, PenInTeeth= FFdata$PenInTeeth, NoPen=FFdata$NoPen, PenInLips=FFdata$PenInLips, ConditionIndex = FFdata$ConditionIndex)

# Stan data for null model
FFdataClean2<- data.frame(HappinessRating= FFdata$HappinessRating, Participant= FFdata$Participant)


# Null model: one mean for all conditions
FFmodel0 <- map(
			alist( HappinessRating ~ dnorm(mu, sigma), 
			mu <- a,
			a ~ dnorm(50, 50),
			sigma ~ dunif(0, 100)
			), data= FFdata ) 
			
# Alternative model 1: different mean for each condition (dumy variable coding)	
FFmodel1 <- map(
			alist( HappinessRating ~ dnorm(mu, sigma), 
			mu <- a1*PenInTeeth + a2*NoPen + a3*PenInLips,
			a1 ~ dnorm(50, 50),
			a2 ~ dnorm(50, 50),
			a3 ~ dnorm(50, 50),
			sigma ~ dunif(0, 100)
			), data= FFdata ) 
			
# Alternative model 1: different mean for each condition	
FFmodel2 <- map(
			alist( HappinessRating ~ dnorm(mu, sigma), 
			mu <- a[ConditionIndex],
			a[ConditionIndex] ~ dnorm(50, 50),
			sigma ~ dunif(0, 100)
			), data= FFdata ) 

# Alternative model 2: different means and sds for each condition			
FFmodel3 <- map(
			alist( HappinessRating ~ dnorm(mu, sigma), 
			mu <- a[ConditionIndex],
			a[ConditionIndex] ~ dnorm(50, 50),
			sigma[ConditionIndex] ~ dunif(0, 100)
			), data= FFdata ) 		
	
# Multilevel: different means and sds for each condition			
FFMultilevel <- ulam(
			alist( HappinessRating ~ dnorm(mu, sigma), 
			mu <- b1[Participant]*PenInTeeth + b2[Participant]*NoPen + b3[Participant]*PenInLips,
			c(b1, b2, b3)[Participant] ~ multi_normal(c(mean_teeth, mean_none, mean_lips), Rho, sigma_part),
			mean_teeth ~dnorm(50, 50),
			mean_none ~dnorm(50, 50),
			mean_lips ~dnorm(50, 50),
			sigma_part ~ dunif(0, 100),
			Rho ~ dlkjcorr(2),
			sigma <- c1*PenInTeeth + c2*NoPen + c3*PenInLips,
			c1 ~ dunif(0, 100),
			c2 ~ dunif(0, 100),
			c3 ~ dunif(0, 100)
			), data= FFdataClean, log_lik=TRUE, chains=4 ) 
			
save(FFMultilevel, file="Multilevel")			
			
print(compare(FFmodel0, FFmodel2, FFmodel3, FFMultilevel))

			
# Multilevel: equal means and sds for each condition			
FFMultilevelNull <- ulam(
			alist( HappinessRating ~ dnorm(mu, sigma), 
			mu <- a[Participant],
			a[Participant] ~ dnorm(grandmean, sigma_part),
			grandmean ~dnorm(50, 50),
			sigma_part ~ dunif(0, 100),
			sigma ~ dunif(0, 100)
			), data= FFdataClean2, log_lik=TRUE, chains=4  ) 
			

# Multilevel: teeth and lips deviate from none in predicted directions and sds for each condition			
FFMultilevelOrder <- ulam(
			alist( HappinessRating ~ dnorm(mu, sigma), 
			mu <- (a2[Participant]+a1[Participant])*PenInTeeth + a2[Participant]*NoPen + (a2[Participant]-a3[Participant])*PenInLips,
			c(a1, a2, a3)[Participant] ~ multi_normal(c(deviation_teeth, mean_none, deviation_lips), Rho, sigma_part),
			deviation_teeth ~dunif(0, 10),
			mean_none ~dnorm(50, 50),
			deviation_lips ~dunif(0, 10),
			sigma_part ~ dunif(0, 100),
			Rho ~ dlkjcorr(2),
			sigma <- c1*PenInTeeth + c2*NoPen + c3*PenInLips,
			c1 ~ dunif(0, 100),
			c2 ~ dunif(0, 100),
			c3 ~ dunif(0, 100)
			), data= FFdataClean, log_lik=TRUE, chains=4  ) 
			
save(FFMultilevelOrder, file="MultilevelOrder")	


# Multilevel: teeth and lips deviate from each other in predicted directions and sds for each condition			
FFMultilevelOrder2 <- ulam(
			alist( HappinessRating ~ dnorm(mu, sigma), 
			mu <- (a3[Participant]+a1[Participant])*PenInTeeth + a2[Participant]*NoPen + a3[Participant]*PenInLips,
			c(a1, a2, a3)[Participant] ~ multi_normal(c(deviation_teeth, mean_none, mean_lips), Rho, sigma_part),
			deviation_teeth ~dunif(0, 10),
			mean_none ~dnorm(50, 50),
			mean_lips ~dnorm(50, 50),
			sigma_part ~ dunif(0, 100),
			Rho ~ dlkjcorr(2),
			sigma <- c1*PenInTeeth + c2*NoPen + c3*PenInLips,
			c1 ~ dunif(0, 100),
			c2 ~ dunif(0, 100),
			c3 ~ dunif(0, 100)
			), data= FFdataClean, log_lik=TRUE, chains=4 ) 
			
save(FFMultilevelOrder2, file="MultilevelOrder2")