#     ___   ____    ____  ______    __   _______       ___      .__   __.   ______  _______ 
#    /   \  \   \  /   / /  __  \  |  | |       \     /   \     |  \ |  |  /      ||   ____|
#   /  ^  \  \   \/   / |  |  |  | |  | |  .--.  |   /  ^  \    |   \|  | |  ,----'|  |__   
#  /  /_\  \  \      /  |  |  |  | |  | |  |  |  |  /  /_\  \   |  . `  | |  |     |   __|  
# /  _____  \  \    /   |  `--'  | |  | |  '--'  | /  _____  \  |  |\   | |  `----.|  |____ 
#/__/     \__\  \__/     \______/  |__| |_______/ /__/     \__\ |__| \__|  \______||_______|

#Analysis script avoidance data general/stimulus difference analysis, mediation analysis and reactiontime
#written by Madeleine Mueller, January 2024

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

######################################
#AVOIDANCE MAIN ANALYSIS  
#steps as outcome measure
######################################
#read in data
OSLO_grid<-read.csv(file='/Users/mmueller/Documents/socialLearn/ana/R/d2_avoidance.csv',head=TRUE,sep=",")
#exclude fence trials
clearFence<-OSLO_grid[OSLO_grid$CS != "3", ]
#exclude both CS door trials
onlyCS<-clearFence[clearFence$CS != "4", ]

onlyCS$trial<-as.factor(onlyCS$trial)
onlyCS$CS<-as.factor(onlyCS$CS)
onlyCS$pathlength<-as.factor(onlyCS$pathlength)
onlyCS$code<-onlyCS$code

#create model
totalsteps<-(lmer(steps~(1|code)+CS*pathlength+trial, data=onlyCS,control = lmerControl(optimizer = "bobyqa")))
#ANOVA based on model
Anova(totalsteps, type="3", test="F" )
#post-hoc tests
emmeans(totalsteps, list(pairwise ~ CS:pathlength), adjust = "none")
#plot effect
plot(effect("CS*pathlength",totalsteps))
plot(effect("CS",totalsteps))


######################################
#AVOIDANCE STIMULUS DIFFERENCE ANALYSIS 
#steps as outcome measure
######################################
# Calculate the mean steps per participant for each combination of CS and pathlength
mean_steps <- onlyCS %>%
  group_by(code, CS, pathlength) %>%
  summarise(mean_steps = mean(steps))

# Calculate the difference between mean steps for CS=1 and CS=2 for each pathlength
csdiff <- mean_steps %>%
  pivot_wider(names_from = CS, values_from = mean_steps) %>%
  mutate(diff_mean_steps = `1` - `2`)

csdiff$pathlength<-as.factor(csdiff$pathlength)

#create model
csdiff_model<-(lmer(diff_mean_steps~(1|code)+pathlength, data=csdiff,control = lmerControl(optimizer = "bobyqa")))
#ANOVA based on model
Anova(csdiff_model, type="3", test="F" )
#post-hoc test
emmeans(csdiff_model, list(pairwise ~ pathlength), adjust = "none")
#plot effect
plot(effect("pathlength",csdiff_model))


######################################
#AVOIDANCE WALLFACTOR ANALYSIS 
#wallfactor as outcome measure (mean wallfactor over whole trial)
######################################
#create model
wall2<-(lmer(wallfactor~(1|code)+CS*pathlength+trial, data=onlyCS, control = lmerControl(optimizer = "bobyqa")))
#ANOVA based on model
Anova(wall2, type="3", test="F" )
#post-hoc test
emmeans(wall2, list(pairwise ~ trial), adjust = "none")
#plot effects
plot(effect("CS*pathlength",wall2))
plot(effect("pathlength",wall2))
plot(effect("CS",wall2))


######################################
#AVOIDANCE MEDIATION ANALYSIS 
######################################
#libraries
library(mediation)
library(lme4)
set.seed(2014)
data("framing", package = "mediation")

t.test(onlyCS$wallfactor~onlyCS$CS, var.equal = TRUE, alternative = "two.sided")

#example:
#1. MEDIATOR MODEL:
#lmer=mediator model type, as linear model "lm" would also be possible. As GLM, we could use glm/bayesglm/glmer. 
#wallmean is the mediator of this model
#then after~ comes the treatment, which is CS for us, we would need a binary treatment (0=CS+, 1=CS-)
#this translates to "With treatment of CS- or without"
#following that, we can add pre-treatment covariates, e.g. age/gender... 
medoslo.fit<- lmer(wallfactor~(1|code)+CS, data=onlyCS)
summary(medoslo.fit)

#2. OUTCOME MODEL:
#here again lm, because our outcome variable "steps" is not binary, for binary, we would need glmer
#with this model, we expect that CS- decision is increasing wallfactor, which leads to more steps
outoslo.fit <- lmer(steps ~ (1|code)+CS+wallfactor , data = onlyCS)
summary(outoslo.fit)
#
#ACME=average causal mediation effects ("indirect effect")
#ADE=average direct effects
#1000 simulations is default, but it is recommended, but one can choose a higher number, if results vary too much
#it is recommended to use bootstrap, if sims>1000 and computing power is not a problem/results should be very similar without boot
#if boot='FALSE' a quasi-Bayesian approximation is used for confidence intervals; if 'TRUE' nonparametric bootstrap will be used. 

medoslo.out <- mediate(medoslo.fit, outoslo.fit, treat = "CS", mediator = "wallfactor",boot=FALSE, sims = 1000)
summary(medoslo.out)
plot(medoslo.out)

######################################
#AVOIDANCE REACTIONTIME ANALYSIS 
######################################
#read in data
OSLOrt<-read.csv(file='/Users/mmueller/Documents/socialLearn/ana/R/d2_avoidance_reactiontime.csv',head=TRUE,sep=",")

#PRE DECISION
#exclude steps during decision
OSLOrtwodec <- subset(OSLOrt, !grepl("dec", prepost))
#exclude steps after decision
OSLOrtwodec <- subset(OSLOrtwodec, !grepl("post", prepost))
df <- subset(OSLOrtwodec, reactiontime <= 4)
df <- subset(df, reactiontime >0)
rm(OSLOrt)
rm(OSLOrtwodec)

#create model
pre_rt_single<-(lmer(reactiontime~(1|code)+Csdoor*pathlength+steps+trial, data=df,control = lmerControl(optimizer = "bobyqa")))
#ANOVA based on model
Anova(pre_rt_single, type="3", test="F" )
#post-hoc tests
emmeans(pre_rt_single, list(pairwise ~ trial), adjust = "none")
#plot effect
plot(effect("steps",pre_rt_single))

#POST DECISION
#exclude steps during decision
OSLOrtwodec <- subset(OSLOrt, !grepl("dec", prepost))
#exclude steps before decision
OSLOrtwodec <- subset(OSLOrtwodec, !grepl("pre", prepost))
dfpost <- subset(OSLOrtwodec, reactiontime <= 4)
dfpost <- subset(dfpost, reactiontime >0)
rm(OSLOrt)
rm(OSLOrtwodec)

#create model
post_rt_single<-(lmer(reactiontime~(1|code)+Csdoor*pathlength+steps+trial, data=dfpost,control = lmerControl(optimizer = "bobyqa")))
#ANOVA based on model
Anova(post_rt_single, type="3", test="F" )
#plot effect
plot(effect("steps",post_rt_single))