'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")

'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)

'Seczion: Decision phase.'
'Load data set'
ROI_Decision <- read_sav("D:/ROI_Decision_phase.sav")

'Factors'
ROI_Decision$Provocation_previous.f <- factor(ROI_Decision$Provocation_previous)
is.factor(ROI_Decision$Provocation_previous.f)

ROI_Decision$Externalization.f <- factor(ROI_Decision$Externalization)
is.factor(ROI_Decision$Externalization.f)



'ROI: Provocation x externalization interaction'
ACC <- aov_ez("id_subject","ROI_ACC",ROI_Decision,between=c("Externalization.f"),within=c("Provocation_previous.f"), anova_table = list(es = "none"))

nice(ACC, es = "pes", observed = NULL,
     correction = c("GG", "HF", "none"), MSE = TRUE,
     sig_symbols = c(" +", " *", " **", " ***"))

p.adjust( .014, method = "bonferroni", n = 2)
p.adjust(  .008, method = "bonferroni", n = 2)

OFC <- aov_ez("id_subject","ROI_OFC",ROI_Decision,between=c("Externalization.f"),within=c("Provocation_previous.f"), anova_table = list(es = "none"))

nice(OFC, es = "pes", observed = NULL,
     correction = c("GG", "HF", "none"), MSE = TRUE,
     sig_symbols = c(" +", " *", " **", " ***"))
p.adjust(  .100, method = "bonferroni", n = 2)
p.adjust( .036, method = "bonferroni", n = 2)

'Section: Parametric modulation of aggression response-related activation'
'Load data set'
ROI_Parametric <- read_sav("D:/ROI_Parametric_modulation.sav")

'Factors'
ROI_Parametric$gender.f <- factor(ROI_Parametric$gender)
is.factor(ROI_Parametric$gender.f)

ROI_Parametric$Externalization.f <- factor(ROI_Parametric$Externalization)
is.factor(ROI_Parametric$Externalization.f)

'ROI: Provocation x externalization interaction'
L_Precentral <- lm(ROI_L_4a ~ Ztestosterone_pooled*gender.f, data=ROI_Parametric)
anova(L_Precentral)
eta_sq(L_Precentral)

'Without variational case'
no_variational_case <- subset(ROI_Parametric,ROI_Parametric$id_subject!=213)

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

no_variational_case$Externalization.f <- factor(no_variational_case$Externalization)
is.factor(no_variational_case$Externalization.f)

L_Precentral_no_variational_case <- lm(ROI_L_4a ~ Ztestosterone_pooled*gender.f, data=no_variational_case)
anova(L_Precentral_no_variational_case)
eta_sq(L_Precentral_no_variational_case)

'Section: Feedback phase.'
'Load data set'
ROI_Feedback <- read_sav("D:/ROI_Feedback_phase.sav")

'Factors'
ROI_Feedback$Provocation_previous.f <- factor(ROI_Feedback$Provocation_previous)
is.factor(ROI_Feedback$Provocation_previous.f)

ROI_Feedback$Externalization.f <- factor(ROI_Feedback$Externalization)
is.factor(ROI_Feedback$Externalization.f)

'ROI: Provocation x externalization interaction'
ACC <- aov_ez("id_subject","ROI_ACC",ROI_Feedback,between=c("Externalization.f"),within=c("Provocation_previous.f"), anova_table = list(es = "none"))

nice(ACC, es = "pes", observed = NULL,
     correction = c("GG", "HF", "none"), MSE = TRUE,
     sig_symbols = c(" +", " *", " **", " ***"))


Insula <- aov_ez("id_subject","ROI_insula",ROI_Feedback,between=c("Externalization.f"),within=c("Provocation_previous.f"), anova_table = list(es = "none"))

nice(Insula, es = "pes", observed = NULL,
     correction = c("GG", "HF", "none"), MSE = TRUE,
     sig_symbols = c(" +", " *", " **", " ***"))

p.adjust( .034 , method = "bonferroni", n = 2)

Amygdala <- aov_ez("id_subject","ROI_amygdala",ROI_Feedback,between=c("Externalization.f"),within=c("Provocation_previous.f"), anova_table = list(es = "none"))

nice(Amygdala, es = "pes", observed = NULL,
     correction = c("GG", "HF", "none"), MSE = TRUE,
     sig_symbols = c(" +", " *", " **", " ***"))