#.______          ___   .___________. __  .__   __.   _______      _______.
#|   _  \        /   \  |           ||  | |  \ |  |  /  _____|    /       |
#|  |_)  |      /  ^  \ `---|  |----`|  | |   \|  | |  |  __     |   (----`
#|      /      /  /_\  \    |  |     |  | |  . `  | |  | |_ |     \   \    
#|  |\  \----./  _____  \   |  |     |  | |  |\   | |  |__| | .----)   |   
#| _| `._____/__/     \__\  |__|     |__| |__| \__|  \______| |_______/    

#Analysis script Rating data "Fear/Discomfort/Physical Arousal"-Ratings (day1 and day 2) and "danger/safe" ratings (day 2)
#written by Madeleine Mueller, January 2024

#libraries
library("lme4")
library("car")
library("effects")
library("emmeans")
library("dplyr")

#   _______       ___   ____    ____     __                                  
#  |       \     /   \  \   \  /   /    /_ |                                 
#  |  .--.  |   /  ^  \  \   \/   /      | |                                 
#  |  |  |  |  /  /_\  \  \_    _/       | |                                 
#  |  '--'  | /  _____  \   |  |         | |                                 
#  |_______/ /__/     \__\  |__|         |_|                                 
##########################
#FEAR/DISCOMFORT/PHYSICAL AROUSAL
##########################
#read in data
O_D1<-read.csv(file='/Users/mmueller/Documents/socialLearn/ana/R/d1_FearRating.csv',head=TRUE,sep=",")

O_D1$stim<-as.factor(O_D1$stim)
O_D1$time<-as.factor(O_D1$prepost)
O_D1$sub<-O_D1$sub

#create models
d1.discomfort<-(lmer(disc_rat~(1|sub)+stim*time, data=O_D1,control = lmerControl(optimizer = "bobyqa")))
d1.fear<-(lmer(fear_rat~(1|sub)+stim*time, data=O_D1,control = lmerControl(optimizer = "bobyqa")))
d1.phys<-(lmer(phys_rat~(1|sub)+stim*time, data=O_D1,control = lmerControl(optimizer = "bobyqa")))

#ANOVAs based on models
Anova(d1.discomfort, type="3", test="F" )
Anova(d1.fear, type="3", test="F" )
Anova(d1.phys, type="3", test="F" )

#post-hoc tests
emmeans(d1.discomfort, list(pairwise ~ stim:time), adjust = "none")
emmeans(d1.fear, list(pairwise ~ stim:time), adjust = "none")
emmeans(d1.phys, list(pairwise ~ stim:time), adjust = "none")

#plot interaction effect
plot(effect("stim*time",d1.discomfort))
plot(effect("stim*time",d1.fear))
plot(effect("stim*time",d1.phys))

##########################
#FEAR/DISCOMFORT/PHYSICAL AROUSAL --> CS DISCRIMINATION CS+ - CS-
##########################
# Create a new data frame with the differences
csdiff <- O_D1 %>%
  # Group by subject and prepost
  group_by(sub, prepost) %>%
  # Calculate the differences for each variable
  summarise(diff_disc_rat = diff(disc_rat),
            diff_fear_rat = diff(fear_rat),
            diff_phys_rat = diff(phys_rat))%>%
  # Adjust the differences to be CSP - CSM
  mutate(across(starts_with("diff"), ~ -1 * .))

csdiff$time<-as.factor(csdiff$prepost)

#create model, now with CS differences as outcome measure, therefore only time as fixed effect
d1.ddiff<-(lmer(diff_disc_rat~(1|sub)+time, data=csdiff,control = lmerControl(optimizer = "bobyqa")))
d1.fdiff<-(lmer(diff_fear_rat~(1|sub)+time, data=csdiff,control = lmerControl(optimizer = "bobyqa")))
d1.pdiff<-(lmer(diff_phys_rat~(1|sub)+time, data=csdiff,control = lmerControl(optimizer = "bobyqa")))

#ANOVAs based on models
Anova(d1.ddiff, type="3", test="F" )
Anova(d1.fdiff, type="3", test="F" )
Anova(d1.pdiff, type="3", test="F" )

#post-hoc tests
emmeans(d1.ddiff, list(pairwise ~ time), adjust = "none")
emmeans(d1.fdiff, list(pairwise ~ time), adjust = "none")
emmeans(d1.pdiff, list(pairwise ~ time), adjust = "none")

#plot effect
plot(effect("time",d1.ddiff))
plot(effect("time",d1.fdiff))
plot(effect("time",d1.pdiff))

#    _______       ___   ____    ____     ___   
#   |       \     /   \  \   \  /   /    |__ \  
#   |  .--.  |   /  ^  \  \   \/   /        ) | 
#   |  |  |  |  /  /_\  \  \_    _/        / /  
#   |  '--'  | /  _____  \   |  |         / /_  
#   |_______/ /__/     \__\  |__|        |____| 

##########################
#FEAR/DISCOMFORT/PHYSICAL AROUSAL
##########################
#read in data
O_D2<-read.csv(file='/Users/mmueller/Documents/socialLearn/ana/R/d2_FearRating.csv',head=TRUE,sep=",")

O_D2$stim<-as.factor(O_D2$stim)
O_D2$time<-as.factor(O_D2$prepost)
O_D2$sub<-O_D2$sub

#create models
d2.disc<-(lmer(disc_rat~(1|sub)+stim*time, data=O_D2,control = lmerControl(optimizer = "bobyqa")))
d2.fear<-(lmer(fear_rat~(1|sub)+stim*time, data=O_D2,control = lmerControl(optimizer = "bobyqa")))
d2.phys<-(lmer(phys_rat~(1|sub)+stim*time, data=O_D2,control = lmerControl(optimizer = "bobyqa")))

#ANOVAs based on models
Anova(d2.disc, type="3", test="F" )
Anova(d2.fear, type="3", test="F" )
Anova(d2.phys, type="3", test="F" )

#post-hoc tests
emmeans(d2.disc, list(pairwise ~ stim:time), adjust = "none")
emmeans(d2.fear, list(pairwise ~ stim:time), adjust = "none")
emmeans(d2.phys, list(pairwise ~ stim:time), adjust = "none")

#plot interaction effect
plot(effect("stim*time",d2.disc))
plot(effect("stim*time",d2.fear))
plot(effect("stim*time",d2.phys))


################
##SAFE/DANGER Ratings
################
#now we have generalised stimuli
#read in data
O_D2s<-read.csv(file='/Users/mmueller/Documents/socialLearn/ana/R/d2_Safe_Ratings.csv',head=TRUE,sep=",")

O_D2s$stim<-as.factor(O_D2s$stim)
O_D2s$block<-as.factor(O_D2s$block)
O_D2s$sub<-O_D2s$sub

#create model
d2.safe<-(lmer(safe_rat~(1|sub)+stim*block, data=O_D2s,control = lmerControl(optimizer = "bobyqa")))
#ANOVA based on model
Anova(d2.safe, type="3", test="F" )
#post-hoc tests
emmeans(d2.safe, list(pairwise ~ stim*block), adjust = "none")
#options(max.print=2000)

#plot effects
plot(effect("stim",d2.safe))
plot(effect("stim*block",d2.safe))