install.packages("lme4")
install.packages("sjstats")
install.packages("haven")
install.packages("lmerTest")
install.packages("MuMIn")
install.packages("ggplot2")
install.packages("etaSquared")
install.packages("effsize")
install.packages("prettyR")
installed.packages("psych")

'Load libraries'
library(haven)
library(sjstats)
library(lmerTest)
library(lme4)
library(MuMIn)
library(ggplot2)
library(etaSquared)
library(emmeans)
library(effsize)
library(MASS) 
library(prettyR)
library(psych)

'Load data set'
descriptives <- read_sav("D:/data_aggregated.sav")

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

descriptives$gender.f <- factor(descriptives$gender)
is.factor(descriptives$gender.f)

'Creating Table 1 (externalization: 1 = high, 0 = low; gender: 1 = female, 0 = male; responder: 1 = responder, 0 = nonresponder)'
Table_1_age <- aggregate(descriptives$age, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_age <- do.call(data.frame, Table_1_age) 
print(Table_1_age)

testosterone_pooled <- subset(descriptives,descriptives$testosterone_pooled > 0) 
Table_1_testosterone_pooled <- aggregate(testosterone_pooled$testosterone_pooled, by = list(Gender = testosterone_pooled$gender.f, Externalization = testosterone_pooled$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_testosterone_pooled <- do.call(data.frame, Table_1_testosterone_pooled)
print(Table_1_testosterone_pooled)

Table_1_TriPM_disinhibition <- aggregate(descriptives$TriPM_Disinhibition, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_TriPM_disinhibition <- do.call(data.frame, Table_1_TriPM_disinhibition)
print(Table_1_TriPM_disinhibition)

Table_1_TriPM_Meanness <- aggregate(descriptives$TriPM_Meanness, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_TriPM_Meanness <- do.call(data.frame, Table_1_TriPM_Meanness)
print(Table_1_TriPM_Meanness)

Table_1_TriPM_Boldness <- aggregate(descriptives$TriPM_Boldness, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_TriPM_Boldness <- do.call(data.frame, Table_1_TriPM_Boldness)
print(Table_1_TriPM_Boldness)

Table_1_K_FAF_spontanous <- aggregate(descriptives$K_FAF_spontanous, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_K_FAF_spontanous <- do.call(data.frame, Table_1_K_FAF_spontanous)
print(Table_1_K_FAF_spontanous)

Table_1_K_FAF_reactive <- aggregate(descriptives$K_FAF_reactive, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_K_FAF_reactive <- do.call(data.frame, Table_1_K_FAF_reactive)
print(Table_1_K_FAF_reactive)

Table_1_K_FAF_excit <- aggregate(descriptives$K_FAF_excit, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_K_FAF_excit <- do.call(data.frame, Table_1_K_FAF_excit)
print(Table_1_K_FAF_excit)

Table_1_K_FAF_self <- aggregate(descriptives$K_FAF_self, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_K_FAF_self <- do.call(data.frame, Table_1_K_FAF_self)
print(Table_1_K_FAF_self)

Table_1_K_FAF_inhibition <- aggregate(descriptives$K_FAF_inhibition, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_K_FAF_inhibition <- do.call(data.frame, Table_1_K_FAF_inhibition)
print(Table_1_K_FAF_inhibition)

Table_1_BUSS_anger <- aggregate(descriptives$BUSS_anger, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_BUSS_anger <- do.call(data.frame, Table_1_BUSS_anger)
print(Table_1_BUSS_anger)

Table_1_BUSS_phys_agg <- aggregate(descriptives$BUSS_phys_agg, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_BUSS_phys_agg <- do.call(data.frame, Table_1_BUSS_phys_agg)
print(Table_1_BUSS_phys_agg)

Table_1_BUSS_hostility <- aggregate(descriptives$BUSS_hostility, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_BUSS_hostility <- do.call(data.frame, Table_1_BUSS_hostility)
print(Table_1_BUSS_hostility)

Table_1_BUSS_verb_Agg <- aggregate(descriptives$BUSS_verb_Agg, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_BUSS_verb_Agg <- do.call(data.frame, Table_1_BUSS_verb_Agg)
print(Table_1_BUSS_verb_Agg)

Table_1_BUSS_SUM <- aggregate(descriptives$BUSS_SUM, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_BUSS_SUM <- do.call(data.frame, Table_1_BUSS_SUM)
print(Table_1_BUSS_SUM)

Table_1_BGHA_Child <- aggregate(descriptives$BGHA_Child, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_BGHA_Child <- do.call(data.frame, Table_1_BGHA_Child)
print(Table_1_BGHA_Child)

Table_1_BGHA_Adules <- aggregate(descriptives$BGHA_Adules, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_BGHA_Adules <- do.call(data.frame, Table_1_BGHA_Adules)
print(Table_1_BGHA_Adules)

Table_1_BGHA_Adult <- aggregate(descriptives$BGHA_Adult, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_BGHA_Adult <- do.call(data.frame, Table_1_BGHA_Adult)
print(Table_1_BGHA_Adult)

Table_1_RPQ_proactive <- aggregate(descriptives$RPQ_proactive, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_RPQ_proactive <- do.call(data.frame, Table_1_RPQ_proactive)
print(Table_1_RPQ_proactive)


Table_1_RPQ_reactive <- aggregate(descriptives$RPQ_reactive, by = list(Gender = descriptives$gender.f, Externalization = descriptives$externalization.f), FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
Table_1_RPQ_reactive <- do.call(data.frame, Table_1_RPQ_reactive)
print(Table_1_RPQ_reactive)

'welch tests between the two externalization groups'
t.test(descriptives$age  ~descriptives$externalization.f)
cohen.d(descriptives$age, descriptives$externalization.f)
p.adjust(0.2986, method = "bonferroni", n = 20)

t.test(descriptives$testosterone_pooled  ~descriptives$externalization.f)
cohen.d(descriptives$testosterone_pooled ,descriptives$externalization.f)
p.adjust(0.4977, method = "bonferroni", n = 20)

t.test(descriptives$TriPM_Disinhibition  ~descriptives$externalization.f)
cohen.d(descriptives$TriPM_Disinhibition ,descriptives$externalization.f)
p.adjust(2.2e-16, method = "bonferroni", n = 20)

t.test(descriptives$TriPM_Meanness  ~descriptives$externalization.f)
cohen.d(descriptives$TriPM_Meanness  ,descriptives$externalization.f)
p.adjust(2.2e-16, method = "bonferroni", n = 20)

t.test(descriptives$TriPM_Boldness  ~descriptives$externalization.f)
cohen.d(descriptives$TriPM_Boldness  ,descriptives$externalization.f)
p.adjust(0.0191, method = "bonferroni", n = 20)

t.test(descriptives$K_FAF_spontanous  ~descriptives$externalization.f)
cohen.d(descriptives$K_FAF_spontanous ,descriptives$externalization.f)
p.adjust(3.95e-10, method = "bonferroni", n = 20)

t.test(descriptives$K_FAF_reactive  ~descriptives$externalization.f)
cohen.d(descriptives$K_FAF_reactive  ,descriptives$externalization.f)
p.adjust(3.971e-10, method = "bonferroni", n = 20)

t.test(descriptives$K_FAF_excit  ~descriptives$externalization.f)
cohen.d(descriptives$K_FAF_excit  ,descriptives$externalization.f)
p.adjust(1.887e-06, method = "bonferroni", n = 20)

t.test(descriptives$K_FAF_self  ~descriptives$externalization.f)
cohen.d(descriptives$K_FAF_self  ,descriptives$externalization.f)
p.adjust(1.508e-07, method = "bonferroni", n = 20)

t.test(descriptives$K_FAF_inhibition  ~descriptives$externalization.f)
cohen.d(descriptives$K_FAF_inhibition ,descriptives$externalization.f)
p.adjust(0.002172, method = "bonferroni", n = 20)

t.test(descriptives$BUSS_anger  ~descriptives$externalization.f)
cohen.d(descriptives$BUSS_anger ,descriptives$externalization.f)
p.adjust(2.659e-06, method = "bonferroni", n = 20)

t.test(descriptives$BUSS_phys_agg  ~descriptives$externalization.f)
cohen.d(descriptives$BUSS_phys_agg ,descriptives$externalization.f)
p.adjust(2.812e-06, method = "bonferroni", n = 20)

t.test(descriptives$BUSS_hostility  ~descriptives$externalization.f)
cohen.d(descriptives$BUSS_hostility ,descriptives$externalization.f)
p.adjust(3.196e-06, method = "bonferroni", n = 20)

t.test(descriptives$BUSS_verb_Agg  ~descriptives$externalization.f)
cohen.d(descriptives$BUSS_verb_Agg  ,descriptives$externalization.f)
p.adjust(0.0003076, method = "bonferroni", n = 20)

t.test(descriptives$BUSS_SUM  ~descriptives$externalization.f)
cohen.d(descriptives$BUSS_SUM  ,descriptives$externalization.f)
p.adjust(2.712e-08, method = "bonferroni", n = 20)

t.test(descriptives$BGHA_Child  ~descriptives$externalization.f)
cohen.d(descriptives$BGHA_Child  ,descriptives$externalization.f)
p.adjust(0.01308, method = "bonferroni", n = 20)

t.test(descriptives$BGHA_Adules  ~descriptives$externalization.f)
cohen.d(descriptives$BGHA_Adules  ,descriptives$externalization.f)
p.adjust(1.469e-05, method = "bonferroni", n = 20)

t.test(descriptives$BGHA_Adult  ~descriptives$externalization.f)
cohen.d(descriptives$BGHA_Adult  ,descriptives$externalization.f)
p.adjust(0.009582, method = "bonferroni", n = 20)

t.test(descriptives$RPQ_proactive  ~descriptives$externalization.f)
cohen.d(descriptives$RPQ_proactive  ,descriptives$externalization.f)
p.adjust(0.0009686, method = "bonferroni", n = 20)

t.test(descriptives$RPQ_reactive  ~descriptives$externalization.f)
cohen.d(descriptives$RPQ_reactive  ,descriptives$externalization.f)
p.adjust(0.0001561, method = "bonferroni", n = 20)