#---
#title: "Artistic Activism Study 1"
#author: "blinded for review"
#version: 2024-05-17
#---

rm(list=ls())

####Preparations#####

setwd("/…/study1")

library(readr)
library(rstatix)
library(nlme)
library(sjPlot)
library(emmeans)
library(tidyr)
library(tidyverse)
library(ggpubr)
library(ggplot2)
library(dplyr)
library(effectsize)
library(psych)
library(sjmisc)
library(descr)
library(car)

#####import dataset
ds.activism.clean <- read_delim("data_activism_study1.csv","\t", escape_double = FALSE, trim_ws = TRUE, show_col_types = FALSE)

#####Rename Variables#####
ds.activism.analysis <- ds.activism.clean %>% rename(
  CASE	=	CASE,
  t.started	=	STARTED,
  group	=	RD01,
  rotation	=	RD02,
  gender	=	SD01,
  yearofbirth	=	SD02_01,
  education	=	SD03,
  ideology	=	SD04,
  cc.attitude	=	SD05,
  art.interest	=	SD06,
  mc.disrupt1	=	S101,
  mc.art1	=	S102,
  extreme1	=	S103,
  immoral1	=	S104,
  emotion1	=	S105,
  identification1	=	S106,
  supportprotest1	=	S107,
  supportcause1	=	S108,
  mc.disrupt2	=	S201,
  mc.art2	=	S202,
  extreme2	=	S203,
  immoral2	=	S204,
  emotion2	=	S205,
  identification2	=	S206,
  supportprotest2	=	S207,
  supportcause2	=	S208,
  mc.disrupt3	=	S301,
  mc.art3	=	S302,
  extreme3	=	S303,
  immoral3	=	S304,
  emotion3	=	S305,
  identification3	=	S306,
  supportprotest3	=	S307,
  supportcause3	=	S308,
  t.finished	=	LASTDATA,
)

d1 <- subset(ds.activism.analysis, select = c(CASE, t.started, group, rotation, gender, gender.cat, age.cat, yearofbirth, education, ideology, cc.attitude, art.interest, mc.disrupt1, mc.art1, extreme1, immoral1, emotion1, identification1, supportprotest1, supportcause1, mc.disrupt2, mc.art2, extreme2, immoral2, emotion2, identification2, supportprotest2, supportcause2, mc.disrupt3, mc.art3, extreme3, immoral3, emotion3, identification3, supportprotest3, supportcause3, t.finished))
rm(ds.activism.analysis, ds.activism.clean)

#####Recode EG-Groups to factor#####
d1$group <- as.factor(d1$group)
d1$group.disrupt <- as.factor(dplyr::recode(d1$group, '1'='nd', '2'='md', '3'='nd', '4'='md', '5'='hd', '6'='hd'))
d1$group.art <- as.factor(dplyr::recode(d1$group, '1'='noa', '2'='noa', '3'='art', '4'='art', '5'='noa', '6'='art'))
d1$group <- interaction(d1$group.disrupt, d1$group.art)
d1$CASE <- as.factor(d1$CASE)
d1$gender <- as.factor(d1$gender)
d1$rotation <- as.factor(d1$rotation)

#####Recode Age & Gender & education#####
d1$age <- 2021-d1$yearofbirth
d1$gender <- ifelse(d1$gender==4, 3, d1$gender)
d1$edu.cat <- as.factor(ifelse((d1$education<=4 | d1$education == 7), 'else', 'upp/sec'))
d1$gender.cat <- as.factor(ifelse((d1$gender.cat==1), 'female', 'male'))

#####add aggregated variables for DVs######
d1 <- d1 %>%
  mutate(mean.mc.disrupt = rowMeans(select(., c(mc.disrupt1, mc.disrupt2, mc.disrupt3))),
         mean.mc.art = rowMeans(select(., c(mc.art1, mc.art2, mc.art3))),
         mean.extreme = rowMeans(select(., c(extreme1, extreme2, extreme3))),
         mean.immoral = rowMeans(select(., c(immoral1, immoral2, immoral3))),
         mean.emotion = rowMeans(select(., c(emotion1, emotion2, emotion3))),
         mean.identification = rowMeans(select(., c(identification1, identification2, identification3))),
         mean.supportprotest = rowMeans(select(., c(supportprotest1, supportprotest2, supportprotest3))),
         mean.supportcause = rowMeans(select(., c(supportcause1, supportcause2, supportcause3)))
  )

####Descriptives####
## group sizes
table(d1$group)
# NOTE: each group > 30

## Gender
table(d1$gender)
prop.table(table(d1$gender))
round(freq(d1$gender), 1 )

table(d1$gender.cat)
prop.table(table(d1$gender.cat))
round(freq(d1$gender.cat), 1 )

## Age
describe(d1$age)

##Gender x age
table(d1$age.cat, d1$gender.cat)

## Education
table(d1$education)
prop.table(table(d1$education))
round(freq(d1$education), 1 )

table(d1$edu.cat)
prop.table(table(d1$edu.cat))
round(freq(d1$edu.cat), 1 )

###DVs
##perceived as disruptive
describe(d1$mean.mc.disrupt)
hist(d1$mean.mc.disrupt)

## perceived as art
describe(d1$mean.mc.art)
hist(d1$mean.mc.art)

## perceived as extreme
describe(d1$mean.extreme)
hist(d1$mean.extreme)

## perceived as immoral
describe(d1$mean.immoral)
hist(d1$mean.immoral)

## perceived as emotion
describe(d1$mean.emotion)
hist(d1$mean.emotion)

## perceived as identification
describe(d1$mean.identification)
hist(d1$mean.identification)

## support protest
describe(d1$mean.supportprotest)
hist(d1$mean.supportprotest)

## support cause
describe(d1$mean.supportcause)
hist(d1$mean.supportcause)

###other Variables###

## Art interest
describe(d1$art.interest)
hist(d1$art.interest)

## Ideology/political attitude
describe(d1$ideology)
hist(d1$ideology)

## Climate change attitude
describe(d1$cc.attitude)
hist(d1$cc.attitude)

##Interest in art
describe(d1$art.interest)
hist(d1$art.interest)

###intercorrelations (wide-dataset)
d1 %>%
  select(mean.mc.disrupt, mean.mc.art, mean.extreme, mean.immoral, mean.emotion, mean.identification, mean.supportprotest, mean.supportcause, ideology, cc.attitude, art.interest) %>%
  corr.test()

# Correlation repeated measures
d1 %>%
  select(extreme1, extreme2, extreme3) %>%
  corr.test()

d1 %>%
  select(immoral1, immoral2, immoral3) %>%
  corr.test()

d1 %>%
  select(emotion1, emotion2, emotion3) %>%
  corr.test()

d1 %>%
  select(identification1, identification2, identification3) %>%
  corr.test()

d1 %>%
  select(supportprotest1, supportprotest2, supportprotest3) %>%
  corr.test()

d1 %>%
  select(supportcause1, supportcause2, supportcause3) %>%
  corr.test()


### Equal distribution of descriptives/variables in experimental groups?
## Gender
table(d1$group, d1$gender)
round(prop.table(table(d1$group, d1$gender), margin = 1), 2)
chisq.test(d1$gender, d1$group)

table(d1$group, d1$gender.cat)
round(prop.table(table(d1$group, d1$gender.cat), margin = 1), 2)
chi2_genderxgroup <- chisq.test(d1$gender.cat, d1$group)
chi2_genderxgroup
sqrt(chi2_genderxgroup$statistic/length(d1$group))

## Education
#full education
table(d1$group, d1$education)
round(prop.table(table(d1$group, d1$education), margin = 1), 2)
chi2_edu_full <- chisq.test(d1$group, d1$education)
chi2_edu_full
sqrt(chi2_edu_full$statistic/length(d1$group))

#education simple
table(d1$group, d1$edu.cat)
round(prop.table(table(d1$group, d1$edu.cat), margin = 1), 2)
chi2_edu_simple <- chisq.test(d1$group, d1$edu.cat)
chi2_edu_simple
sqrt(chi2_edu_simple$statistic/length(d1$group))

## Age
describeBy(d1$age, d1$group)
distr.group.age <- aov(d1$age ~ d1$group)
summary(distr.group.age) 
eta_squared(distr.group.age)

## Ideology/political orientation
describeBy(d1$ideology, d1$group)
distr.group.ideology <- aov(d1$ideology ~ d1$group)
summary(distr.group.ideology) 

## Climate change attitude
describeBy(d1$cc.attitude, d1$group)
distr.group.cc.attitude <- aov(d1$cc.attitude ~ d1$group)
summary(distr.group.cc.attitude)

## interest in art
describeBy(d1$art.interest, d1$group)
distr.group.art.interest <- aov(d1$art.interest ~ d1$group)
summary(distr.group.art.interest)

####Transform dataset to long####

#####wide to long#####
d1.long <- d1 %>% 
  pivot_longer(
    !c(CASE, t.started, group, rotation, gender, gender.cat, age.cat, yearofbirth, age, education, ideology, cc.attitude, art.interest, t.finished, group.disrupt, group.art, edu.cat, mean.mc.disrupt, mean.mc.art, mean.extreme, mean.immoral, mean.emotion, mean.identification, mean.supportprotest, mean.supportcause), 
    names_to = c(".value", "reportID"), 
    names_pattern = "(.*)(.)"
  )

### recode report names
d1.long$reportID <- as.factor(dplyr::recode(d1.long$reportID, '1'='report1', '2'='report2', '3'='report3'))

#####Manipulation Check: (MC_disruption)##### 
##Mean values (MC_disruption)
mc.disrupt.meansby.artcondition <- d1.long %>%
  group_by(group.art) %>%
  get_summary_stats(mc.disrupt, type = "mean_sd")

mc.disrupt.meansby.disrupttcondition <- d1.long %>%
  group_by(group.disrupt) %>%
  get_summary_stats(mc.disrupt, type = "mean_sd")

mc.disrupt.meansby.eggroups <- d1.long %>%
  group_by(group.art, group.disrupt) %>%
  get_summary_stats(mc.disrupt, type = "mean_sd")

mc.disrupt.meansby.disrupttcondition
mc.disrupt.meansby.artcondition
mc.disrupt.meansby.eggroups

##Test (MC_disruption) 
#Main analytical model: Predictors: Disrupt x Art
mc.disrupt.int.model <- stats::aov(mc.disrupt ~ group.disrupt*group.art + Error(CASE), data = d1.long)
base::summary(mc.disrupt.int.model)
effectsize::eta_squared(mc.disrupt.int.model,partial = TRUE, alternative = "two.sided",ci = 0.95)

##post-hoc tukey test (MC_disruption)
emmeans(mc.disrupt.int.model, list(pairwise ~ group.disrupt), adjust = "tukey")
emmeans(mc.disrupt.int.model, list(pairwise ~ group.art), adjust = "tukey")
emmeans(mc.disrupt.int.model, list(pairwise ~ group.disrupt*group.art), adjust = "tukey")

#explore differences between reports: Disrupt x Art x ReportID (MC_disruption)
mc.disrupt.int.model.report <- stats::aov(mc.disrupt ~ group.disrupt*group.art*reportID + Error(CASE), data = d1.long)
base::summary(mc.disrupt.int.model.report)
effectsize::eta_squared(mc.disrupt.int.model.report,partial = TRUE, alternative = "two.sided",ci = 0.95)

#explore covariates: Disrupt x Art x Covariates (MC_disruption)
mc.disrupt.by.controls <- stats::aov(mc.disrupt ~ group.disrupt*group.art*reportID +
                                       age +
                                       gender.cat + 
                                       edu.cat + 
                                       group.art +
                                       ideology + 
                                       cc.attitude + 
                                       art.interest +
                                       Error(CASE), data = d1.long)

base::summary(mc.disrupt.by.controls)
effectsize::eta_squared(mc.disrupt.by.controls,partial = TRUE, alternative = "two.sided",ci = 0.95)

#####Manipulation Check: (MC_artwork)#####
##Mean values (MC_artwork)
mc.art.meansby.artcondition <- d1.long %>%
  group_by(group.art) %>%
  get_summary_stats(mc.art, type = "mean_sd")

mc.art.meansby.disrupttcondition <- d1.long %>%
  group_by(group.disrupt) %>%
  get_summary_stats(mc.art, type = "mean_sd")

mc.art.meansby.eggroups <- d1.long %>%
  group_by(group.art, group.disrupt) %>%
  get_summary_stats(mc.art, type = "mean_sd")

mc.art.meansby.disrupttcondition
mc.art.meansby.artcondition
mc.art.meansby.eggroups

##Test (MC_artwork)
#Main analytical model: Predictors: Disrupt x Art
mc.art.int.model <- stats::aov(mc.art ~ group.disrupt*group.art + Error(CASE), data = d1.long)
base::summary(mc.art.int.model)
effectsize::eta_squared(mc.art.int.model,partial = TRUE, alternative = "two.sided",ci = 0.95)

##post-hoc tukey test (MC_artwork)
emmeans(mc.art.int.model, list(pairwise ~ group.disrupt), adjust = "tukey")
emmeans(mc.art.int.model, list(pairwise ~ group.art), adjust = "tukey")
emmeans(mc.art.int.model, list(pairwise ~ group.disrupt*group.art), adjust = "tukey")

#explore differences between reports: Disrupt x Art x ReportID (MC_artwork)
mc.art.int.model.report <- stats::aov(mc.art ~ group.disrupt*group.art*reportID + Error(CASE), data = d1.long)
base::summary(mc.art.int.model.report)
effectsize::eta_squared(mc.art.int.model.report, partial = TRUE, alternative = "two.sided", ci = 0.95)

#explore covariates: Disrupt x Art x Covariates (MC_artwork)
mc.art.by.controls <- stats::aov(mc.art ~ group.disrupt*group.art*reportID +
                                       age +
                                       gender.cat + 
                                       edu.cat + 
                                       group.art +
                                       ideology + 
                                       cc.attitude + 
                                       art.interest +
                                       Error(CASE), data = d1.long)

base::summary(mc.art.by.controls)
effectsize::eta_squared(mc.art.by.controls, partial = TRUE, alternative = "two.sided", ci = 0.95)

#####Analysis (DV_extreme)#####
##Mean values (DV_extreme)
extreme.meansby.artcondition <- d1.long %>%
  group_by(group.art) %>%
  get_summary_stats(extreme, type = "mean_sd")

extreme.meansby.disrupttcondition <- d1.long %>%
  group_by(group.disrupt) %>%
  get_summary_stats(extreme, type = "mean_sd")

extreme.meansby.eggroups <- d1.long %>%
  group_by(group.art, group.disrupt) %>%
  get_summary_stats(extreme, type = "mean_sd")

extreme.meansby.disrupttcondition
extreme.meansby.artcondition
extreme.meansby.eggroups

##Test (DV_extreme)
#Main analytical model: Predictors: Disrupt x Art
extreme.int.model <- stats::aov(extreme ~ group.disrupt*group.art + Error(CASE), data = d1.long)
base::summary(extreme.int.model)
effectsize::eta_squared(extreme.int.model, partial = TRUE, alternative = "two.sided", ci = 0.95)

##post-hoc tukey test (DV_extreme)
emmeans(extreme.int.model, list(pairwise ~ group.disrupt), adjust = "tukey")
emmeans(extreme.int.model, list(pairwise ~ group.art), adjust = "tukey")
emmeans(extreme.int.model, list(pairwise ~ group.disrupt*group.art), adjust = "tukey")

#explore differences between reports: Disrupt x Art x ReportID (DV_extreme)
extreme.int.model.report <- stats::aov(extreme ~ group.disrupt*group.art*reportID + Error(CASE), data = d1.long)
base::summary(extreme.int.model.report)
effectsize::eta_squared(extreme.int.model.report, partial = TRUE, alternative = "two.sided", ci = 0.95)

#explore covariates: Disrupt x Art x Covariates (DV_extreme)
extreme.by.controls <- stats::aov(extreme ~ group.disrupt*group.art*reportID +
                                   age +
                                   gender.cat + 
                                   edu.cat + 
                                   group.art +
                                   ideology + 
                                   cc.attitude + 
                                   art.interest +
                                   Error(CASE), data = d1.long)

base::summary(extreme.by.controls)
effectsize::eta_squared(extreme.by.controls, partial = TRUE, alternative = "two.sided", ci = 0.95)

#####Analysis (DV_immoral)#####
##Mean values (DV_immoral)
immoral.meansby.artcondition <- d1.long %>%
  group_by(group.art) %>%
  get_summary_stats(immoral, type = "mean_sd")

immoral.meansby.disrupttcondition <- d1.long %>%
  group_by(group.disrupt) %>%
  get_summary_stats(immoral, type = "mean_sd")

immoral.meansby.eggroups <- d1.long %>%
  group_by(group.art, group.disrupt) %>%
  get_summary_stats(immoral, type = "mean_sd")

immoral.meansby.disrupttcondition
immoral.meansby.artcondition
immoral.meansby.eggroups

##Test (DV_immoral)
#Main analytical model: Predictors: Disrupt x Art
immoral.int.model <- stats::aov(immoral ~ group.disrupt*group.art + Error(CASE), data = d1.long)
base::summary(immoral.int.model)
effectsize::eta_squared(immoral.int.model, partial = TRUE, alternative = "two.sided", ci = 0.95)

##post-hoc tukey test (DV_immoral)
emmeans(immoral.int.model, list(pairwise ~ group.disrupt), adjust = "tukey")
emmeans(immoral.int.model, list(pairwise ~ group.art), adjust = "tukey")
emmeans(immoral.int.model, list(pairwise ~ group.disrupt*group.art), adjust = "tukey")

#explore differences between reports: Disrupt x Art x ReportID (DV_immoral)
immoral.int.model.report <- stats::aov(immoral ~ group.disrupt*group.art*reportID + Error(CASE), data = d1.long)
base::summary(immoral.int.model.report)
effectsize::eta_squared(immoral.int.model.report, partial = TRUE, alternative = "two.sided", ci = 0.95)

#explore covariates: Disrupt x Art x Covariates (DV_immoral)
immoral.by.controls <- stats::aov(immoral ~ group.disrupt*group.art*reportID +
                                    age +
                                    gender.cat + 
                                    edu.cat + 
                                    group.art +
                                    ideology + 
                                    cc.attitude + 
                                    art.interest +
                                    Error(CASE), data = d1.long)

base::summary(immoral.by.controls)
effectsize::eta_squared(immoral.by.controls, partial = TRUE, alternative = "two.sided", ci = 0.95)

#####Analysis (DV_emotion)#####
##Mean values (DV_emotion)
emotion.meansby.artcondition <- d1.long %>%
  group_by(group.art) %>%
  get_summary_stats(emotion, type = "mean_sd")

emotion.meansby.disrupttcondition <- d1.long %>%
  group_by(group.disrupt) %>%
  get_summary_stats(emotion, type = "mean_sd")

emotion.meansby.eggroups <- d1.long %>%
  group_by(group.art, group.disrupt) %>%
  get_summary_stats(emotion, type = "mean_sd")

emotion.meansby.disrupttcondition
emotion.meansby.artcondition
emotion.meansby.eggroups


##Test (DV_emotion)
#Main analytical model: Predictors: Disrupt x Art
emotion.int.model <- stats::aov(emotion ~ group.disrupt*group.art + Error(CASE), data = d1.long)
base::summary(emotion.int.model)
effectsize::eta_squared(emotion.int.model, partial = TRUE, alternative = "two.sided", ci = 0.95)

##post-hoc tukey test (DV_emotion)
emmeans(emotion.int.model, list(pairwise ~ group.disrupt), adjust = "tukey")
emmeans(emotion.int.model, list(pairwise ~ group.art), adjust = "tukey")
emmeans(emotion.int.model, list(pairwise ~ group.disrupt*group.art), adjust = "tukey")

#explore differences between reports: Disrupt x Art x ReportID (DV_emotion)
emotion.int.model.report <- stats::aov(emotion ~ group.disrupt*group.art*reportID + Error(CASE), data = d1.long)
base::summary(emotion.int.model.report)
effectsize::eta_squared(emotion.int.model.report, partial = TRUE, alternative = "two.sided", ci = 0.95)

#explore covariates: Disrupt x Art x Covariates (DV_emotion)
emotion.by.controls <- stats::aov(emotion ~ group.disrupt*group.art*reportID +
                                    age +
                                    gender.cat + 
                                    edu.cat + 
                                    group.art +
                                    ideology + 
                                    cc.attitude + 
                                    art.interest +
                                    Error(CASE), data = d1.long)

base::summary(emotion.by.controls)
effectsize::eta_squared(emotion.by.controls, partial = TRUE, alternative = "two.sided", ci = 0.95)

#####Analysis (DV_identification)#####
##Mean values (DV_identification)
identification.meansby.artcondition <- d1.long %>%
  group_by(group.art) %>%
  get_summary_stats(identification, type = "mean_sd")

identification.meansby.disrupttcondition <- d1.long %>%
  group_by(group.disrupt) %>%
  get_summary_stats(identification, type = "mean_sd")

identification.meansby.eggroups <- d1.long %>%
  group_by(group.art, group.disrupt) %>%
  get_summary_stats(identification, type = "mean_sd")

identification.meansby.disrupttcondition
identification.meansby.artcondition
identification.meansby.eggroups

##Test (DV_identification)
#Main analytical model: Predictors: Disrupt x Art
identification.int.model <- stats::aov(identification ~ group.disrupt*group.art + Error(CASE), data = d1.long)
base::summary(identification.int.model)
effectsize::eta_squared(identification.int.model, partial = TRUE, alternative = "two.sided", ci = 0.95)

##post-hoc tukey test (DV_identification)
emmeans(identification.int.model, list(pairwise ~ group.disrupt), adjust = "tukey")
emmeans(identification.int.model, list(pairwise ~ group.art), adjust = "tukey")
emmeans(identification.int.model, list(pairwise ~ group.disrupt*group.art), adjust = "tukey")

#explore differences between reports: Disrupt x Art x ReportID (DV_identification)
identification.int.model.report <- stats::aov(identification ~ group.disrupt*group.art*reportID + Error(CASE), data = d1.long)
base::summary(identification.int.model.report)
effectsize::eta_squared(identification.int.model.report, partial = TRUE, alternative = "two.sided", ci = 0.95)

#explore covariates: Disrupt x Art x Covariates (DV_identification)
identification.by.controls <- stats::aov(identification ~ group.disrupt*group.art*reportID +
                                    age +
                                    gender.cat + 
                                    edu.cat + 
                                    group.art +
                                    ideology + 
                                    cc.attitude + 
                                    art.interest +
                                    Error(CASE), data = d1.long)

base::summary(identification.by.controls)
effectsize::eta_squared(identification.by.controls, partial = TRUE, alternative = "two.sided", ci = 0.95)

#####Analysis (DV_supportprotest)#####
##Mean values (DV_supportprotest)
supportprotest.meansby.artcondition <- d1.long %>%
  group_by(group.art) %>%
  get_summary_stats(supportprotest, type = "mean_sd")

supportprotest.meansby.disrupttcondition <- d1.long %>%
  group_by(group.disrupt) %>%
  get_summary_stats(supportprotest, type = "mean_sd")

supportprotest.meansby.eggroups <- d1.long %>%
  group_by(group.art, group.disrupt) %>%
  get_summary_stats(supportprotest, type = "mean_sd")

supportprotest.meansby.disrupttcondition
supportprotest.meansby.artcondition
supportprotest.meansby.eggroups

##Test (DV_supportprotest)
#Main analytical model: Predictors: Disrupt x Art
supportprotest.int.model <- stats::aov(supportprotest ~ group.disrupt*group.art + Error(CASE), data = d1.long)
base::summary(supportprotest.int.model)
effectsize::eta_squared(supportprotest.int.model, partial = TRUE, alternative = "two.sided", ci = 0.95)

##post-hoc tukey test (DV_supportprotest)
emmeans(supportprotest.int.model, list(pairwise ~ group.disrupt), adjust = "tukey")
emmeans(supportprotest.int.model, list(pairwise ~ group.art), adjust = "tukey")
emmeans(supportprotest.int.model, list(pairwise ~ group.disrupt*group.art), adjust = "tukey")

#explore differences between reports: Disrupt x Art x ReportID (DV_supportprotest)
supportprotest.int.model.report <- stats::aov(supportprotest ~ group.disrupt*group.art*reportID + Error(CASE), data = d1.long)
base::summary(supportprotest.int.model.report)
effectsize::eta_squared(supportprotest.int.model.report, partial = TRUE, alternative = "two.sided", ci = 0.95)

#explore covariates: Disrupt x Art x Covariates (DV_supportprotest)
supportprotest.by.controls <- stats::aov(supportprotest ~ group.disrupt*group.art*reportID +
                                           age +
                                           gender.cat + 
                                           edu.cat + 
                                           group.art +
                                           ideology + 
                                           cc.attitude + 
                                           art.interest +
                                           Error(CASE), data = d1.long)

base::summary(supportprotest.by.controls)
effectsize::eta_squared(supportprotest.by.controls, partial = TRUE, alternative = "two.sided", ci = 0.95)

#####Analysis (DV_supportcause)#####
##Mean values (DV_supportcause)
supportcause.meansby.artcondition <- d1.long %>%
  group_by(group.art) %>%
  get_summary_stats(supportcause, type = "mean_sd")

supportcause.meansby.disrupttcondition <- d1.long %>%
  group_by(group.disrupt) %>%
  get_summary_stats(supportcause, type = "mean_sd")

supportcause.meansby.eggroups <- d1.long %>%
  group_by(group.art, group.disrupt) %>%
  get_summary_stats(supportcause, type = "mean_sd")

supportcause.meansby.disrupttcondition
supportcause.meansby.artcondition
supportcause.meansby.eggroups

##Test (DV_supportcause)
#Main analytical model: Predictors: Disrupt x Art
supportcause.int.model <- stats::aov(supportcause ~ group.disrupt*group.art + Error(CASE), data = d1.long)
base::summary(supportcause.int.model)
effectsize::eta_squared(supportcause.int.model, partial = TRUE, alternative = "two.sided", ci = 0.95)

##post-hoc tukey test (DV_supportcause)
emmeans(supportcause.int.model, list(pairwise ~ group.disrupt), adjust = "tukey")
emmeans(supportcause.int.model, list(pairwise ~ group.art), adjust = "tukey")
emmeans(supportcause.int.model, list(pairwise ~ group.disrupt*group.art), adjust = "tukey")

#explore differences between reports: Disrupt x Art x ReportID (DV_supportcause)
supportcause.int.model.report <- stats::aov(supportcause ~ group.disrupt*group.art*reportID + Error(CASE), data = d1.long)
base::summary(supportcause.int.model.report)
effectsize::eta_squared(supportcause.int.model.report, partial = TRUE, alternative = "two.sided", ci = 0.95)

#explore covariates: Disrupt x Art x Covariates (DV_supportprotest)
supportcause.by.controls <- stats::aov(supportcause ~ group.disrupt*group.art*reportID +
                                           age +
                                           gender.cat + 
                                           edu.cat + 
                                           group.art +
                                           ideology + 
                                           cc.attitude + 
                                           art.interest +
                                           Error(CASE), data = d1.long)

base::summary(supportcause.by.controls)
effectsize::eta_squared(supportcause.by.controls, partial = TRUE, alternative = "two.sided", ci = 0.95)

####Figure####

#####violin plots#####

library(ggplot2)
library(ggpubr) 

# List of variables
variables <- c("mc.disrupt", "mc.art", "extreme", "immoral", "emotion", "identification", "supportprotest", "supportcause")

# Create an empty list to store the plots
plots_list <- list()

# Vector of corresponding labels for the Y-axis
y_axis_labels <- c("Perc. as disr.", "Perc. as art", "Extreme", "Immoral", "Emotion", 
                   "Identification", "Support Protest", "Support Cause")

# Loop through each variable
for (i in seq_along(variables)) {
  variable <- variables[i]
  label <- y_axis_labels[i]
  
  # Create the plot for the current variable
  current_plot <- ggplot(d1.long, aes(x = group.disrupt, y = .data[[variable]], fill = group.art)) +
    geom_violin(trim = TRUE, scale = "width") +
    stat_summary(
      fun.data = "mean_cl_normal",
      geom = "pointrange",
      colour = "white",
      position = position_dodge(width = 0.9), 
      size = 0.1,
    ) +
    scale_x_discrete(labels = c("Peaceful", "Disobedient", "Violent")) +
    scale_fill_manual(values = c("noa" = "grey55", "art" = "grey30")) +
    scale_y_continuous(breaks = 1:6, limits = c(1, 6)) +
    theme_light(base_size = 12) + 
    theme(legend.position = "none") + 
    coord_cartesian(ylim = c(1, 6)) +
    labs(y = variable, x = " ") +
    ylab(label)  # Set Y-axis label for the current variable

  # Store  plot in the list
  plots_list <- c(plots_list, list(current_plot))
}

# Arrange the plots using ggarrange
figure <- ggarrange(plotlist = plots_list, ncol = 2, nrow = 4)

# Print the final plot
print(figure)

# Export the final plot
ggsave("figure1.pdf", figure, height = 19.5, width = 13.5, units = "cm")

#####export data for alternative plots#####
export.meansby.eggroups <- rbind(mc.disrupt.meansby.eggroups, 
                                 mc.art.meansby.eggroups,
                                 extreme.meansby.eggroups,
                                 immoral.meansby.eggroups,
                                 emotion.meansby.eggroups,
                                 identification.meansby.eggroups,
                                 supportprotest.meansby.eggroups,
                                 supportcause.meansby.eggroups)

export.meansby.eggroups

#write.csv(export.meansby.eggroups,"export_meansby_eggroups.csv", row.names=FALSE)