'Install packages'
install.packages("lme4")
install.packages("sjstats")
install.packages("texreg")
install.packages("foreign")
install.packages("MuMIn")
install.packages("lmerTest")
install.packages("simr")
install.packages("multcomp")
install.packages("haven")
install.packages("emmeans")
install.packages("cowplot")
install.packages("simr", dependencies=TRUE)
install.packages("ggplot2")
install.packages("ggsignif")
install.packages("ggplot")
install.packages("gridExtra")
install.packages("sjPlot")

'Load libraries'
library(emmeans)
library(haven)
library(lme4)
library(sjstats)
library(texreg)
library(foreign)
library(MuMIn)
library(multcomp)
library(simr)
library(ggplot2)
library(lmerTest)
library(car)
library(afex)
library(psych)
library(effsize)
library(lsmeans)
library(cowplot)
library(xtable)
library(psychometric)
library(afex)
library(ggsignif)
library(sjPlot)
library(gridExtra)

'Load data set'
data_all <- read_sav("D:/data_all_without_win.sav")
data_all <- subset(data_all,data_all$participant_id!=246)


'Factors'
data_all$provocation_category_in_the_previous_trial.f <- factor(data_all$provocation_category_in_the_previous_trial)
is.factor(data_all$provocation_category_in_the_previous_trial.f)

data_all$externalization.f <- factor(data_all$externalization)
is.factor(data_all$externalization.f)
data_all$deception_check_rater_2.f <- factor(data_all$deception_check_rater_2)
is.factor(data_all$deception_check_rater_2.f)
data_all$gender.f <- factor(data_all$gender)
is.factor(data_all$gender.f)


'Random intercept model'
model_1 <- lmer(aggression_level ~ 1 + (1| participant_id) , data_all)
icc(model_1)
r.squaredGLMM(model_1)

'Random intercept model with fixed effect for provocation'
model_2 <- lmer(aggression_level ~ (1| participant_id) + provocation_category_in_the_previous_trial.f, data_all)
r.squaredGLMM(model_2)

'Random slope model with fixed effect for provocation'
model_3 <- lmer(aggression_level ~ provocation_category_in_the_previous_trial.f + (1 + provocation_category_in_the_previous_trial.f| participant_id) , data_all)
print(model_3)
r.squaredGLMM(model_3)

'Random slope model with all fixed effects (provocation, experimentel condition, gender, K-FAF reactive aggression scale, deception check)'
model_4 <- lmer(aggression_level ~ provocation_category_in_the_previous_trial.f*externalization.f   + ZCortisol_AUCg2*Ztestosterone_pooled*gender.f +  ZK_FAF_reactive + deception_check_rater_2.f + (1+ provocation_category_in_the_previous_trial.f|participant_id), data_all)
anova(model_4)
r.squaredGLMM(model_4)
eta_sq(anova(model_4))

print(model_4)
summary(model_4)

'Likelihood  Ratio (LR) test model_1 vs model_2'
anova(model_1,model_2)
0.5*(1 - pchisq(( 16198.32 -16130.42)*2,1)) + 0.5*(1 - pchisq((16198.32 -16130.42)*2,2))

'Likelihood  Ratio (LR) test model_2 vs model_3 (random provocation effect)'
anova(model_2,model_3)
0.5*(1 - pchisq(( -16039.83+16130.42 )*2,1)) + 0.5*(1 - pchisq((-16039.83+16130.42)*2,2))

'Likelihood  Ratio (LR) test model_3 vs model_4 (fixed effect structure)'
anova(model_3,model_4)
0.5*(1 - pchisq((16039.83 - 16002.26  )*2,1)) + 0.5*(1 - pchisq((16039.83 - 16002.26 )*2,2))

texreg(list(model_1,model_2,model_3,model_4), dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE, label = "tab:4", caption = "Table 4 Four linear mixed models with reactive aggression as dependent variable in monetary mTAP")
'Notes. Model 1 only consists of a random intercept for participant; in model 3 and 4 a random slope for provocation by participant is added; model 4 contains also fixed effects for experimental condition, gender, deception check and the trait reactive aggression scale of the K-FAF; *** p < .001,** p<.01,* p < .05.'

'Without variational case'
no_variational_case <- subset(data_all,data_all$participant_id!=213)

'Factors'
no_variational_case$provocation_category_in_the_previous_trial.f <- factor(no_variational_case$provocation_category_in_the_previous_trial)
is.factor(no_variational_case$provocation_category_in_the_previous_trial.f)

no_variational_case$externalization.f <- factor(no_variational_case$externalization)
is.factor(no_variational_case$externalization.f)
no_variational_case$deception_check_rater_1.f <- factor(no_variational_case$deception_check_rater_1)
is.factor(no_variational_case$deception_check_rater_1.f)
no_variational_case$deception_check_rater_2.f <- factor(no_variational_case$deception_check_rater_2)
is.factor(no_variational_case$deception_check_rater_2.f)
no_variational_case$gender.f <- factor(no_variational_case$gender)
is.factor(no_variational_case$gender.f)

'Random slope model with all fixed effects (provocation, experimentel condition, gender, K-FAF reactive aggression scale, deception check)'
model_4_no_variational_case <- lmer(aggression_level ~ provocation_category_in_the_previous_trial.f*externalization.f   + ZCortisol_AUCg2*Ztestosterone_pooled*gender.f +  ZK_FAF_reactive + deception_check_rater_2.f + (1+ provocation_category_in_the_previous_trial.f|participant_id), no_variational_case)
anova(model_4_no_variational_case)
r.squaredGLMM(model_4_no_variational_case)
eta_sq(anova(model_4_no_variational_case))

'affective response to the stressor: main effect time and externalization for positive and negative affect'
'Load data set'
Affective_response <- read_sav("D:/Affective_response.sav")

Affective_response <- subset(Affective_response,Affective_response$participant_id!=249)
Affective_response <- subset(Affective_response,Affective_response$participant_id!=235)
Affective_response <- subset(Affective_response,Affective_response$participant_id!=242)
Affective_response <- subset(Affective_response,Affective_response$participant_id!=221)
Affective_response <- subset(Affective_response,Affective_response$participant_id!=219)
Affective_response <- subset(Affective_response,Affective_response$participant_id!=227)
Affective_response <- subset(Affective_response,Affective_response$participant_id!=239)
Affective_response <- subset(Affective_response,Affective_response$participant_id!=243)

Affective_response$externalization.f <- factor(Affective_response$externalization)
is.factor(Affective_response$externalization.f)

Affective_response$time.f <- factor(Affective_response$time)
is.factor(Affective_response$time.f)

Pos_Affect <- aov_ez("participant_id","PANAS_PA",Affective_response,between=c("externalization.f"),within=c("time.f"))

nice(Pos_Affect, es = "pes", observed = NULL,
     correction = c("GG", "HF", "none"), MSE = TRUE,
     sig.symbols = c(" +", " *", " **", " ***"))
sink("output.txt")
print(Pos_Affect)
sink()

Neg_Affect <- aov_ez("participant_id","PANAS_NA",Affective_response,between=c("externalization.f"),within=c("time.f"))

nice(Neg_Affect, es = "pes", observed = NULL,
     correction = c("GG", "HF", "none"), MSE = TRUE,
     sig.symbols = c(" +", " *", " **", " ***"))
sink("output.txt")
print(Neg_Affect)
sink()