---
title: "Study 2"
author: "Double Blind"
date: "2024-08-02"
output: html_document
---

```{r setup, include=FALSE}


library(aod)
#install.packages("aod")
library(ggplot2)
library(here)
library(tidyverse)
library(rio)
library(magrittr)
library(skimr)
library(janitor)
library(fastDummies)
#install.packages("fastDummies")
library(reshape2)
library(wesanderson)
library(ggridges)
#library(sundry)
#library(colorblindr)
library(purrr)
#library(lme4)
library(lmerTest)
library(vcd)
library(vcdExtra)
library(car)
library(pander) 
library(xtable)
library(apaTables)
library(dplyr)
library(tidyr)
library(svglite)
library(psych)
library(pwr)
library(olsrr)
library(effectsize)
library(simpleboot)
library(modelr)
library(ggthemes)  # For colorblind-friendly palettes
library(ggpubr)
library(rstatix)
library(jtools)
library(sjPlot)
library(ggpubr)
library(boot)

knitr::opts_chunk$set(include = TRUE, echo = TRUE, warning= FALSE, message = FALSE, error = FALSE)
```

# Loading and preparing data for Analysis

```{r datawork}
data <- import(here("data2", "batch2_150324.csv")) 

data %<>% slice(-1:-2) #getting rid of unnecessary 

#formatting date format to do more filters 
data$EndDate <- as.Date(data$EndDate , format = "%m/%d/%y")

#now I will filter for it to only show data since the day data collection started, taking out pilot and test runs

data_dc <- data%>%
  filter(EndDate > '2024-02-21')%>%
  filter(EndDate < '2024-03-11')%>%
  filter(Status != 'Survey Preview')%>%
  filter(Finished != 'FALSE')%>%
  filter(att_ch_noun == 'Library')%>% # filter att check fails
  filter(att_ch_fruit == 'Orange') %>%
  filter( cc_inf_qy %in% c('Meat Industry,Fossil Fuels Use (planes, cars, etc)', 'Meat Industry', 'Fossil Fuels Use (planes, cars, etc)') |  cc_inf_qn %in% c('Recycling', 'Recycling,None of the above'))%>%
  filter( f_inf_qy %in% c('Keep kitchen clean,Store raw meat separately',  'Keep kitchen clean', 'Store raw meat separately') |  f_inf_qn %in% c('Store all foods at room temperature,Don’t clean kitchen utensils', 'Store all foods at room temperature', 'Don’t clean kitchen utensils'))%>%
    filter( mh_inf_qy %in% c('Have a good sleep,Do exercise', 'Have a good sleep', 'Do exercise') |  mh_inf_qn %in% c('None of the above'))

```

### Now lets prepare variables for analysis


```{r datawork2}


## conditions: Hypocrisy/Control ; Left/Right 

data_a <- data_dc %>%
  mutate(condition_h = case_when(judgment_hlf_1  == "" & judgment_hlm_1 == "" & judgment_hlj_1 == "" & judgment_hrf_1  == "" & judgment_hrm_1 == "" & judgment_hrj_1 == "" ~ 'Control'))%>%
  mutate(condition_hyp = case_when(is.na(condition_h) ~ 'Hypocrisy'))%>%
  mutate(condition_hc = as.factor(case_when(condition_hyp== 'Hypocrisy' ~ 'Hypocrisy', condition_h == 'Control' ~ 'Control')))%>%
  mutate(condition_l = case_when(judgment_hrf_1  == "" & judgment_hrm_1 == "" & judgment_hrj_1 == "" & judgment_crf_1  == "" & judgment_crm_1 == "" & judgment_crj_1 == "" ~ 'Left'))%>%
  mutate(condition_r = case_when(is.na(condition_l) ~ 'Right'))%>%
  mutate(condition_lr = as.factor(case_when(condition_l== 'Left' ~ 'Left', condition_r == 'Right' ~ 'Right')))%>%
  dplyr::rename(id = 'Random ID')%>%
  dplyr::rename(id_w = 'workerId')

long_hc <- data_a%>%
  select(judgment_hrf_1, judgment_hrm_1, judgment_hrj_1, judgment_crf_1, judgment_crm_1, judgment_crj_1, judgment_hlf_1, judgment_hlm_1, judgment_hlj_1, judgment_clf_1, judgment_clm_1, judgment_clj_1, id, id_w, condition_hc, condition_lr, starts_with('pol_'), -starts_with('pol_orientation_6'), starts_with('ccp_'), political_party, starts_with('emo_'), starts_with('donation_'))%>%
  pivot_longer( cols = starts_with("judgment_"),
    names_to = "topic",
    names_prefix = "judgment_",
    values_to = "judgment")%>%
  mutate(topic_ok= case_when(topic == 'hrf_1' | topic == 'hlf_1' | topic == 'crf_1' | topic == 'clf_1' ~ 'Fuel', topic == 'hrm_1' | topic == 'hlm_1' | topic == 'crm_1' | topic == 'clm_1' ~ 'Meat', topic == 'hrj_1' | topic == 'hlj_1' | topic == 'crj_1' | topic == 'clj_1' ~ 'Jet'))%>%
  dplyr::rename(topic_old = topic)%>%
  dplyr::rename(topic = topic_ok)%>%
  select(-topic_old)%>%
  filter(judgment != "")%>% 
  mutate_at('topic', factor) %>%
  dplyr::rename(clean_air = 'donation_1', co2_coalition = 'donation_4', goodwill = 'donation_5' )%>%
  mutate(donation = as.numeric(clean_air) + (-1)*as.numeric(co2_coalition) + 0*as.numeric(goodwill)) %>%
  mutate(pol_group = case_when(political_party == 'Democratic party' ~ 'Left', political_party == 'Republican party'  ~ 'Right' , political_party == 'None' | political_party == 'Other:' | political_party == 'Prefer not to say'~ 'Other')) %>%
  mutate(orientation_group = case_when(pol_orientation == 'Very Liberal/Left wing' | pol_orientation == 'Liberal/Left wing' ~ 'Left', pol_orientation == 'Very Conservative/Right wing' | pol_orientation == 'Conservative/Right wing' ~ 'Right', pol_orientation == 'Moderate' | pol_orientation == 'Prefer not to say' | pol_orientation == 'These options don’t describe my beliefs' ~ 'Other' )) %>% 
  mutate(group = case_when(pol_group == 'Left' & orientation_group == 'Left' ~ 'Left', pol_group == 'Right' & orientation_group == 'Right' ~ 'Right',  pol_group == 'Other' | orientation_group == 'Other'  ~ 'Other', pol_group == 'Right' & orientation_group == 'Left'  ~ 'Inconsistent', pol_group == 'Left' & orientation_group == 'Right'  ~ 'Inconsistent'))%>%
  filter(group != 'Other') %>%
  filter(group !=  'Inconsistent') %>%
  mutate(ingroup = case_when(group == condition_lr ~ 1, group != condition_lr ~ 0))%>%
  mutate(ingroup_f = as.factor(case_when(ingroup== 1 ~ 'ingroup', ingroup == 0 ~ 'outgroup')))%>%
   mutate(donation_p = donation/100) %>%
  mutate(across(c(clean_air, co2_coalition, goodwill), as.numeric)) %>%
  mutate(clean_air_p = clean_air/100, co2_coalition_p = co2_coalition/100, goodwill_p = goodwill/100)%>%
  mutate(across(starts_with('ccp_'), as.numeric))%>%
  mutate(cc_policy = (ccp_gasoline_1 + ccp_electric_1 + ccp_meal_1 + ccp_flights_1)/4)



long_likely <- data_a%>%
  select(starts_with('likely_'), id, id_w, condition_hc, condition_lr, starts_with('pol_'), -starts_with('pol_orientation_6'), starts_with('ccp_'), political_party, starts_with('emo_'), starts_with('donation_'))%>%
  pivot_longer( cols = starts_with("likely_"),
    names_to = "topic",
    names_prefix = "likely_",
    values_to = "likely")%>%
  mutate(topic_ok= case_when(topic == 'hrf_1' | topic == 'hlf_1' | topic == 'crf_1' | topic == 'clf_1' ~ 'Fuel', topic == 'hrm_1' | topic == 'hlm_1' | topic == 'crm_1' | topic == 'clm_1' ~ 'Meat', topic == 'hrj_1' | topic == 'hlj_1' | topic == 'crj_1' | topic == 'clj_1' ~ 'Jet'))%>%
  dplyr::rename(topic_old = topic)%>%
  dplyr::rename(topic = topic_ok)%>%
  select(-topic_old)%>%
  filter(likely != "")%>% 
  mutate_at('topic', factor)%>%
  mutate_at('likely', as.numeric)

likely_id_long <- long_likely %>%
  select(likely, id_w) 

long_full2 <- merge(long_likely, long_hc, by = 'id')

likely_id_wide <- likely_id_long%>% distinct(id_w, .keep_all = TRUE)

# IMPORTANT: we can change these filters, as of now, they only KEEP people who answered either both Right in Orientation AND Party, or both Left in Orientation AND Party, and filtering out moderates, other, prefer not to say, none, or if they were inconsistent (eg, left in party but right on orientation), the justification for this is that it does not match the filters we used on mTurk to have a balanced sample between left and right wings, meaning that we cannot entirely trust their data to test the ingroup out group effects especially - when appropriate we will urn analysis including moderates to double check consistency of results

## Below we DON'T get rid of moderates  so we can do other analyses 

long_hc2 <- data_a%>%
  select(judgment_hrf_1, judgment_hrm_1, judgment_hrj_1, judgment_crf_1, judgment_crm_1, judgment_crj_1, judgment_hlf_1, judgment_hlm_1, judgment_hlj_1, judgment_clf_1, judgment_clm_1, judgment_clj_1, id, id_w, condition_hc, condition_lr, starts_with('pol_'), -starts_with('pol_orientation_6'), starts_with('ccp_'), political_party, starts_with('emo_'), starts_with('donation_'))%>%
  pivot_longer( cols = starts_with("judgment_"),
    names_to = "topic",
    names_prefix = "judgment_",
    values_to = "judgment")%>%
  mutate(topic_ok= case_when(topic == 'hrf_1' | topic == 'hlf_1' | topic == 'crf_1' | topic == 'clf_1' ~ 'Fuel', topic == 'hrm_1' | topic == 'hlm_1' | topic == 'crm_1' | topic == 'clm_1' ~ 'Meat', topic == 'hrj_1' | topic == 'hlj_1' | topic == 'crj_1' | topic == 'clj_1' ~ 'Jet'))%>%
  dplyr::rename(topic_old = topic)%>%
  dplyr::rename(topic = topic_ok)%>%
  select(-topic_old)%>%
  filter(judgment != "")%>% 
  mutate_at('topic', factor) %>%
  dplyr::rename(clean_air = 'donation_1', co2_coalition = 'donation_4', goodwill = 'donation_5' )%>%
  mutate(donation = as.numeric(clean_air) + (-1)*as.numeric(co2_coalition) + 0*as.numeric(goodwill)) %>%
  mutate(pol_group = case_when(political_party == 'Democratic party' ~ 'Left', political_party == 'Republican party'  ~ 'Right' , political_party == 'None' | political_party == 'Other:' | political_party == 'Prefer not to say'~ 'Other')) %>%
  mutate(orientation_group = case_when(pol_orientation == 'Very Liberal/Left wing'~ 'Very Left', pol_orientation == 'Liberal/Left wing' ~ 'Left', pol_orientation == 'Very Conservative/Right wing' ~ 'Very Right', pol_orientation == 'Conservative/Right wing' ~ 'Right', pol_orientation == 'Moderate'  ~ 'Moderate', pol_orientation == 'Prefer not to say' | pol_orientation == 'These options don’t describe my beliefs' ~ 'Other' )) %>% 
  mutate(group = case_when(pol_group == 'Left' & orientation_group == 'Left' ~ 'Left', pol_group == 'Right' & orientation_group == 'Right' ~ 'Right',  pol_group == 'Other' | orientation_group == 'Other'  ~ 'Other', pol_group == 'Right' & orientation_group == 'Left'  ~ 'Inconsistent', pol_group == 'Left' & orientation_group == 'Right'  ~ 'Inconsistent'))%>%
   mutate(donation_p = donation/100) %>%
  mutate(across(c(clean_air, co2_coalition, goodwill), as.numeric)) %>%
  mutate(clean_air_p = clean_air/100, co2_coalition_p = co2_coalition/100, goodwill_p = goodwill/100)%>%
  mutate(across(starts_with('ccp_'), as.numeric))%>%
  mutate(cc_policy = (ccp_gasoline_1 + ccp_electric_1 + ccp_meal_1 + ccp_flights_1)/4)

# Now wide format

wide_hc <- long_hc %>%
   pivot_wider(names_from = topic, values_from = judgment)%>%
  mutate(judgment_m = (as.numeric(Fuel) + as.numeric(Jet) + as.numeric(Meat))/3)

wide_hc_lik <- merge(likely_id_wide, wide_hc, by = 'id_w')


# Including moderates wide format 

wide_hc2 <- long_hc2 %>%
   pivot_wider(names_from = topic, values_from = judgment) %>%
   mutate(judgment_m = (as.numeric(Fuel) + as.numeric(Jet) + as.numeric(Meat))/3) %>%
   filter(orientation_group != "Other") %>%
   filter(!is.na(orientation_group)) %>%
   mutate(orientation_group = factor(orientation_group)) %>%
   mutate(orientation_num = case_when(
      orientation_group == 'Very Left' ~ 1,
      orientation_group == 'Left' ~ 2,
      orientation_group == 'Moderate' ~ 3,
      orientation_group == 'Right' ~ 4,
      orientation_group == 'Very Right' ~ 5
   )) %>%
   mutate(orientation_num = as.numeric(orientation_num))



## Let's prep the emotion ratings in a data set for exploratory analysis

emo <- long_hc %>%
  select(starts_with("emo_"), topic, id)%>%
    pivot_longer( cols = starts_with("emo_"),
    names_to = "emo_type",
    names_prefix = "emo_",
    values_to = "emo_rating")%>%
  separate(emo_type, c("topic", "emotion"), sep = "_")%>%
   mutate(topic= case_when(topic == 'hrf' | topic == 'hlf' | topic == 'crf' | topic == 'clf' ~ 'Fuel', topic == 'hrm' | topic == 'hlm' | topic == 'crm' | topic == 'clm' ~ 'Meat', topic == 'hrj' | topic == 'hlj' | topic == 'crj' | topic == 'clj' ~ 'Jet'))%>%
  mutate(emotion = case_when(emotion == '1' ~ 'enthusiastic', emotion == '2' ~ "worried", emotion == '3' ~ "proud", emotion == '4' ~ "angry", emotion == '5' ~ "hopeful", emotion == '6' ~  "frustrated"))%>%
  filter(emo_rating != "")%>%
  mutate_at('emo_rating', as.numeric)%>%
  mutate_at('emotion', as.factor) %>%
  group_by(id, emotion, topic)%>%
  dplyr::summarise(emo_m = mean(emo_rating))%>%
  pivot_wider(names_from = emotion, values_from = emo_m)
  
# long data fram with Likely (ie. plausibility) variable (not including moderates)

long_full <- merge(long_hc, likely_id_wide, by= "id_w")

```

### Quick checks on data

```{r alphas}
#let's check alpha for judgment across topics

judgment <- wide_hc%>%
  select(Fuel, Jet, Meat)%>%
  mutate(across(c(Fuel, Jet, Meat), as.numeric))

alpha(judgment) #alpha of 0.69, acceptable

# let's check alpha for policy support across items

ccp_supp <- wide_hc %>%
  select(ccp_gasoline_1, ccp_electric_1, ccp_meal_1, ccp_flights_1)

alpha(ccp_supp) #alpha of 0.83, very good 
```




## Checking the likely variable (Plausibility) - graphs and summaries (how likely participants think the social media posts would be true)

```{r likely}

sum_likely <- long_likely %>%
  group_by(topic, condition_hc)%>%
  #dplyr::summarise(avg = mean(likely))%>%
   dplyr::summarise(sd = sd(likely))




# Define colorblind-friendly colors
cb_colors <- c("#E69F00", "#56B4E9")  # Orange and Blue from ColorBrewer

# Create the histogram
ggplot(long_full, aes(x = likely, fill = condition_hc)) +
  geom_histogram(binwidth = 1, position = "dodge", color = "black", alpha = 0.8) +
  facet_wrap(~condition_hc) +
  scale_fill_manual(values = cb_colors) +
  theme_minimal(base_size = 14) +
  theme(
    strip.text = element_text(size = 16, face = "bold"),
    legend.position = "none",  # Remove legend since facets indicate groups
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14, face = "bold")
  ) +
  labs(
    x = "Plausibility",
    y = "Count",
    title = "Distribution of Plausibility Scores by Condition"
  )

```





## Obtaining demographics

```{r dems}
data_wide_dem <- data_a %>%
  select(political_party, pol_orientation, cc_important, id, age, gender, race, religion, education, education_parents, )%>%
  mutate(pol_group = case_when(political_party == 'Democratic party' ~ 'Left', political_party == 'Republican party'  ~ 'Right' , political_party == 'None' | political_party == 'Other:' | political_party == 'Prefer not to say'~ 'Other')) %>%
  mutate(orientation_group = case_when(pol_orientation == 'Very Liberal/Left wing' | pol_orientation == 'Liberal/Left wing' ~ 'Left', pol_orientation == 'Very Conservative/Right wing' | pol_orientation == 'Conservative/Right wing' ~ 'Right', pol_orientation == 'Moderate' | pol_orientation == 'Prefer not to say' | pol_orientation == 'These options don’t describe my beliefs' ~ 'Other' )) %>% 
  mutate(group = case_when(pol_group == 'Left' & orientation_group == 'Left' ~ 'Left', pol_group == 'Right' & orientation_group == 'Right' ~ 'Right',  pol_group == 'Other' | orientation_group == 'Other'  ~ 'Other', pol_group == 'Right' & orientation_group == 'Left'  ~ 'Inconsistent', pol_group == 'Left' & orientation_group == 'Right'  ~ 'Inconsistent'))%>%
  filter(group != 'Other') %>%
  filter(group !=  'Inconsistent')
```


## Analyses

### Hipothesis 1.2 

```{r}

## quick plots to observe the data 

hist(wide_hc$judgment_m)

hist(scale(wide_hc$judgment_m))



ggplot(wide_hc, aes(x = condition_hc, y = judgment_m)) +
  geom_violin(draw_quantiles =  0.5) +
  labs(title = "Violin Plot of Judgment Hypocrisy vs Control",
       x = "Condition",
       y = "Mean Jdgment") +
  stat_compare_means(method = "wilcox.test", label.x = 1.3, label.y = 4) 

# judgment is not normally distributed so for a simple menat ocmaprison test we will use wilcoxon (and check with a t test as well)

judgment_m_sd <- wide_hc%>%
  group_by(condition_hc)%>%
  dplyr::summarise(mean = mean(judgment_m), sd = sd(judgment_m))

wilcox_hc <- wide_hc%>%
  pivot_wider(names_from = condition_hc, values_from = judgment_m)%>%
  select(Control, Hypocrisy)

control_j <- wilcox_hc%>%
  select(Control)%>%
  na.omit()

hypocrisy_j <- wilcox_hc%>%
  select(Hypocrisy)%>%
  na.omit()

wilcox.test(control_j$Control, hypocrisy_j$Hypocrisy)

t.test(control_j$Control, hypocrisy_j$Hypocrisy)

### effect size

## prepping the needed info 

wide_hc %<>% mutate_at('condition_hc', as.factor)
test_result <- wilcox.test(judgment_m ~ condition_hc, data= wide_hc)

wilcox.test(judgment_m ~ condition_hc, data= wide_hc, conf.int = T, conf.level = .95)

# Extract W statistic from the test result
W <- test_result$statistic

# Get the sample sizes for each group (assuming 'condition' is a factor with two levels)
n1 <- as.numeric(nrow(hypocrisy_j)) # Number of observations in the Control group
n2 <- as.numeric(nrow(control_j)) # Number of observations in the Hipocrisy group

# Calculate the mean (mu_W) and standard deviation (sigma_W) of W
mu_W <- (n1 * n2) / 2
sigma_W <- sqrt((n1 * n2 * (n1 + n2 + 1)) / 12)

# Calculate the Z statistic (standardized test statistic)
Z <- (W - mu_W) / sigma_W

# Calculate the rank-biserial correlation (effect size)
r <- Z / sqrt(n1 + n2)

# Print the Z statistic and the effect size
Z
r #effect size

```


# Hypothesis 2: Now the interaction

```{r inth2}

#quick data observation 

hist(wide_hc$cc_policy)

# let's run models (judgment averaged across topics)
m_int <- lm(judgment_m ~ condition_hc*cc_policy, data=wide_hc)

m_int_all <- lm(judgment_m ~ condition_hc*cc_policy, data=wide_hc2) #including moderates


# let's run models (topics as random effects)

m_int_2 <- lmer(as.numeric(judgment) ~ condition_hc*cc_policy+ (1|topic), data=long_hc)

m_int_2_all <- lmer(as.numeric(judgment) ~ condition_hc*cc_policy+ (1|topic), data=long_hc2) #including moderates

# results

summary(m_int_2)

summary(m_int_all)

summary(m_int_2_all)

qqnorm(residuals(m_int_2))

qqnorm(residuals(m_int_2_all))

hist(residuals(m_int))

summary(m_int_2)

hist(residuals(m_int_2))


## effect size



int_eta <- eta_squared(m_int)


m_int_eta <- int_eta[3]

d_int <- sqrt((4*m_int_eta)/(1 - m_int_eta))


export_summs(m_int_2, confint = TRUE, robust = "HC0")

# tables

tab_model(m_int_2_all,   pred.labels = c("Intercept", "Condition (Hipocrisy)", "Policy Support", "Hypocrisy*PolicySupp"),
  dv.labels = c("Judgment"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")

tab_model(m_int,   pred.labels = c("Intercept", "Condition (Hipocrisy)", "Policy Support", "Hypocrisy*PolicySupp"),
  dv.labels = c("Judgment"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")

tab_model(m_int_2,   pred.labels = c("Intercept", "Condition (Hipocrisy)", "Policy Support", "Hypocrisy*PolicySupp"),
  dv.labels = c("Judgment"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")
```

### H2: runninG the analysis controlling for Likely (ie. plausibility)

```{r H2_LIKELY}
# let's run models (topics as random effects)

m_2_likely <- lmer(as.numeric(judgment) ~ condition_hc*cc_policy + likely + (1|topic), data=long_full)

summary(m_2_likely)

tab_model(m_2_likely,   pred.labels = c("Intercept", "Condition (Hipocrisy)", "Policy Support", "Plausibility", "Hypocrisy*PolicySupp"),
  dv.labels = c("Judgment"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")
```


## Hypotheses 3 & 4: ingroup-outgroup effect

```{r inoutgroup}

# Hypothesis 3

m_inout <- lmer(as.numeric(judgment) ~ condition_hc + ingroup_f + (1|topic), data = long_hc)
summary(m_inout)

tab_model(m_inout,
  dv.labels = c("Judgment"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")

qqnorm(residuals(m_inout))

hist(residuals(m_inout))

# effect size 

inout_eta <- effectsize::eta_squared(m_inout)

m_inout_eta <- inout_eta$Eta2_partial[2]

d_inout <- sqrt((4*m_inout_eta)/(1 - m_inout_eta))

d_inout

# Hypothesis 4

m_inout2 <- lmer(as.numeric(judgment) ~ condition_hc + ingroup_f *cc_policy + (1|topic), data = long_hc)
summary(m_inout2)

tab_model(m_inout2,
  dv.labels = c("Judgment"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")

qqnorm(residuals(m_inout2))



inout2_eta <- effectsize::eta_squared(m_inout2)

m_inout2_eta <- inout2_eta$Eta2_partial[4]

d_inout2 <- sqrt((4*m_inout2_eta)/(1 - m_inout2_eta))


```

### H3 & 4: running the analysis controlling for Likely (ie. plausibility)

```{r H3_4_LIKELY}
# let's run models (topics as random effects)

m_3_likely <- lmer(as.numeric(judgment) ~ condition_hc + ingroup_f + likely +  (1|topic), data = long_full)

summary(m_3_likely)

tab_model(m_3_likely,
  dv.labels = c("Judgment"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")

m_4_likely <- lmer(as.numeric(judgment) ~ condition_hc + ingroup_f *cc_policy  + likely + (1|topic), data = long_full)

summary(m_4_likely)

tab_model(m_4_likely,
  dv.labels = c("Judgment"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")
```


## Plotting ingroup-outgroup comparison 

```{r plotinout}

# Custom labels for the x-axis
custom_labels <- c(
  "ingroup.Hypocrisy" = "Hypocrisy",
  "outgroup.Hypocrisy" = "Hypocrisy",
  "ingroup.Control" = "Control",
  "outgroup.Control" = "Control"
)

# Plot using ggplot
p <- ggplot(wide_hc, aes(x = interaction(ingroup_f, condition_hc), y = judgment_m, fill = ingroup_f)) +
  geom_violin(outlier.shape = NA) + # Boxplot without outliers
  geom_point(position = position_jitter(width = 0.2), alpha = 0.4)+
   stat_summary(fun = mean, fun.min = mean, fun.max = mean, 
               geom = "crossbar", width = 0.3, color = "darkgray") + # Mean line
  #stat_summary(fun = mean, geom = "point", shape = 20, size = 3, color = "blue") + # Median line
  #stat_boxplot(geom = "errorbar", width = 0.25) + # Error lines
  scale_fill_manual(values = c("#CC79A7", "#009E73")) + # Customize colors
  labs(title = "Leniency to the Political Ingroup", x = "Condition", y = "Judgment (Wrong to Right)", fill = "Group") +
    scale_x_discrete(labels = custom_labels) + # Use custom labels for the x-axis
  theme_minimal()



# Add significance annotation using ggsignif
p + geom_signif(comparisons = list(c("ingroup.Hypocrisy", "ingroup.Control")),
                map_signif_level = TRUE,
                y_position = max(wide_hc$judgment_m) + 0.5) +
    geom_signif(comparisons = list(c("outgroup.Hypocrisy", "outgroup.Control")),
                map_signif_level = TRUE,
                y_position = max(wide_hc$judgment_m) + 1) +
    geom_signif(comparisons = list(c("ingroup.Hypocrisy", "outgroup.Hypocrisy")),
                map_signif_level = TRUE,
                y_position = max(wide_hc$judgment_m) + 1.5) +
    geom_signif(comparisons = list(c("ingroup.Control", "outgroup.Control")), 
                map_signif_level = TRUE,
                y_position = max(wide_hc$judgment_m) + 2) 




```

## Improved version of the plot: 

```{r improved plot}

# Define custom labels
custom_labels <- c(
  "ingroup.Hypocrisy" = "Hypocrisy",
  "outgroup.Hypocrisy" = "Hypocrisy",
  "ingroup.Control" = "Control",
  "outgroup.Control" = "Control"
)

# Compute group means
means_df <- wide_hc %>%
  group_by(group = interaction(ingroup_f, condition_hc)) %>%
  dplyr::summarise(mu_hat = mean(judgment_m), .groups = "drop")

# Base plot
p <- ggplot(wide_hc, aes(x = interaction(ingroup_f, condition_hc), y = judgment_m)) +
  geom_violin(color = "black", fill = NA) +  # No fill, just contour
  geom_jitter(aes(color = ingroup_f), width = 0.2, alpha = 0.6, size = 1.5) +  # Points colored by group
  geom_boxplot(width = 0.15, outlier.shape = NA, alpha = 0.3) +  # Boxplot for spread (pale)
  geom_point(data = means_df, aes(x = group, y = mu_hat), color = "black", size = 3) +  # Mean dot
  geom_text(data = means_df, aes(x = as.numeric(group) + 0.25,, y = mu_hat + 0.15, 
                                 label = paste0("μ̂ = ", round(mu_hat, 2))),
            size = 3, vjust = 0, fontface = "italic") +  # Mean label
  scale_color_manual(values = c(
  "ingroup" = rgb(0.8, 0.2, 0.6),   # CMYK-safe magenta
  "outgroup" = rgb(0.2, 0.6, 0.4)   # CMYK-safe green
)) +
  scale_x_discrete(labels = custom_labels) +
  labs(title = "Leniency to the Political Ingroup",
       x = "Condition", y = "Judgment (Wrong to Right)", color = "Group") +
  theme_minimal() +
  theme(legend.position = "top")

# Add significance bars
p_final <- p +
  geom_signif(comparisons = list(c("ingroup.Hypocrisy", "ingroup.Control")),
              map_signif_level = TRUE,
              y_position = max(wide_hc$judgment_m) + 0.5) +
  geom_signif(comparisons = list(c("outgroup.Hypocrisy", "outgroup.Control")),
              map_signif_level = TRUE,
              y_position = max(wide_hc$judgment_m) + 1) +
  geom_signif(comparisons = list(c("ingroup.Hypocrisy", "outgroup.Hypocrisy")),
              map_signif_level = TRUE,
              y_position = max(wide_hc$judgment_m) + 1.5) +
  geom_signif(comparisons = list(c("ingroup.Control", "outgroup.Control")),
              map_signif_level = TRUE,
              y_position = max(wide_hc$judgment_m) + 2)

# Show the final plot
print(p_final)

# Save the plot in high definition
ggsave("figure2_high_quality.png", p_final, width = 10, height = 6, dpi = 600, bg = "white")


tiff("/Users/tamar/Desktop/R/Hipocrisy/figure_2_cmyk.tiff", width = 7, height = 5, units = "in", res = 300, type = "cairo")
print(p_final)
dev.off()

```



# Hypothesis 1.1 & Hypothesis 4

## Donation analysis

```{r boots}

# Quick mean comparison 
wilcox_d <- wide_hc%>%
  pivot_wider(names_from = condition_hc, values_from = donation_p)%>%
  select(Control, Hypocrisy)

control_d <- wilcox_d%>%
  select(Control)%>%
  na.omit()

hypocrisy_d <- wilcox_d%>%
  select(Hypocrisy)%>%
  na.omit()

wilcox.test(control_d$Control, hypocrisy_d$Hypocrisy)

#model

m_d <- lm(donation_p ~ condition_hc*cc_policy + ingroup_f*cc_policy + condition_hc*ingroup_f, data= wide_hc)

summary(m_d)

tab_model(m_d,
  dv.labels = c("Donation"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")

qqnorm(residuals(m_d))

#quick bootstrap 

library(broom)
bm_d <- lm.boot(m_d, 1000, rows = FALSE) #LOAD LIBRARY

summary(bm_d)

## checcking cc policy effect

m_d0 <- lm(donation_p ~ cc_policy, data= wide_hc)
summary(m_d0)

tab_model(m_d0,
  dv.labels = c("Donation"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")

## effect size for cc policy 

d_eta <- effectsize::eta_squared(m_d)

m_d_eta <- d_eta$Eta2_partial[2]

d_d <- sqrt((4*m_d_eta)/(1 - m_d_eta))


bm_d0 <- lm.boot(m_d0, 1000, rows = FALSE)
summary(bm_d0)

```


# Donation explained by plausibility 

```{r don_plaus}
#aggregating plausibility to see the overall effect of plausibility on donation  

agg <- aggregate(likely ~ id_w, data = long_full, FUN = mean) 

agg$donation_p <- tapply(long_full$donation_p, long_full$id_w, unique)

# linear regression 

model_dp <- lm(donation_p ~ likely, data = agg)

summary(model_dp)

# table 

tab_model(model_dp,
  dv.labels = c("Donation"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")
```

 
 
# Full bootstrap for donation analysis 

```{r bottstrapping}

# Bootstrap 95% CI for regression coefficients

# function to obtain regression weights

bs <- function(formula, data, indices)
{
  d <- data[indices,] # allows boot to select sample
  fit <- lm(formula, data=d)
  return(coef(fit))
}
# bootstrapping with 1000 replications
results <- boot(data=wide_hc, statistic=bs,
   R=1000, formula=donation_p ~ condition_hc*cc_policy + ingroup_f*cc_policy + condition_hc*ingroup_f)



# view results
results
plot(results, index=1) # intercept
plot(results, index=2) # 
plot(results, index=3) # 
plot(results, index=4) # 
plot(results, index=5) #
plot(results, index=6) #
plot(results, index=7) #

# get 95% confidence intervals
boot.ci(results, type="bca", index=1) # 
boot.ci(results, type="bca", index=2) # 
boot.ci(results, type="bca", index=3) # 
boot.ci(results, type="bca", index=4) # 
boot.ci(results, type="bca", index=5) # 
boot.ci(results, type="bca", index=6) # 
boot.ci(results, type="bca", index=7) # 
```

# Exploratory analysis

### Emotions analysis

```{r expemo}

##grouping by positive and negative emotions

negative<- as.data.frame(emo) %>%
  select(angry, frustrated, worried)

alpha(negative)


positive<- as.data.frame(emo) %>%
  select(enthusiastic, hopeful, proud)

alpha(positive)


emo_g <- as.data.frame(emo)%>%
  rowwise()%>%
  mutate(positive = mean(enthusiastic, hopeful, proud), negative = mean(angry, frustrated, worried))

long_emo <- merge(long_hc, emo_g)

#negative 

neg_full <- lmer(negative~condition_hc + ingroup_f + (1|topic), data = long_emo)

summary(neg_full)

tab_model(neg_full,
  dv.labels = c("Negative Emotions"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")


pos_full <- lmer(positive~condition_hc + ingroup_f + (1|topic), data = long_emo)

summary(pos_full)

tab_model(pos_full,
  dv.labels = c("Positive Emotions"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")


mean(long_emo$positive)
sd(long_emo$positive)


##effect sizes

neg_full_eta <- effectsize::eta_squared(neg_full)

neg_full_eta_hc <- neg_full_eta$Eta2_partial[1]

d_neg_full_hc <- sqrt((4*neg_full_eta_hc)/(1 - neg_full_eta_hc))

neg_full_eta_out <- neg_full_eta$Eta2_partial[2]

d_neg_full_out <- sqrt((4*neg_full_eta_out)/(1 - neg_full_eta_out))


### positive

pos_full_eta <- effectsize::eta_squared(pos_full)

pos_full_eta_hc <- pos_full_eta$Eta2_partial[1]

d_pos_full_hc <- sqrt((4*pos_full_eta_hc)/(1 - pos_full_eta_hc))

pos_full_eta_out <- pos_full_eta$Eta2_partial[2]

d_pos_full_out <- sqrt((4*pos_full_eta_out)/(1 - pos_full_eta_out))
```

### Plotting exploratory analysis 

```{r expplot}
#negative emotions

p_neg <- ggplot(long_emo, aes(x = interaction(ingroup_f, condition_hc), y = negative, fill = ingroup_f)) +
  geom_violin(outlier.shape = NA) + # Boxplot without outliers
  geom_point(position = position_jitter(width = 0.2), alpha = 0.4)+
   stat_summary(fun = mean, fun.min = mean, fun.max = mean, 
               geom = "crossbar", width = 0.3, color = "darkgray") + # Mean line
  #stat_summary(fun = mean, geom = "point", shape = 20, size = 3, color = "blue") + # Median line
  #stat_boxplot(geom = "errorbar", width = 0.25) + # Error lines
  scale_fill_manual(values = c("#CC57A7", "#001E73")) + # Customize colors
  labs(title = "Effects on negative emotions", x = "Condition", y = "Negative emotions (low to high)", fill = "Group") +
    scale_x_discrete(labels = custom_labels) + # Use custom labels for the x-axis
  theme_minimal()



# Add significance annotation using ggsignif
p_neg + geom_signif(comparisons = list(c("ingroup.Hypocrisy", "ingroup.Control")),
                map_signif_level = TRUE,
                y_position = max(long_emo$negative) + 1) +
    geom_signif(comparisons = list(c("outgroup.Hypocrisy", "outgroup.Control")),
                map_signif_level = TRUE,
                y_position = max(long_emo$negative) + 4 ) +
    geom_signif(comparisons = list(c("ingroup.Hypocrisy", "outgroup.Hypocrisy")),
                map_signif_level = TRUE,
                y_position = max(long_emo$negative) + 7) +
    geom_signif(comparisons = list(c("ingroup.Control", "outgroup.Control")), 
                map_signif_level = TRUE,
                y_position = max(long_emo$negative) + 10) 

# positive emotions 

p_pos <- ggplot(long_emo, aes(x = interaction(ingroup_f, condition_hc), y = positive, fill = ingroup_f)) +
  geom_violin(outlier.shape = NA) + # Boxplot without outliers
  geom_point(position = position_jitter(width = 0.2), alpha = 0.4)+
   stat_summary(fun = mean, fun.min = mean, fun.max = mean, 
               geom = "crossbar", width = 0.3, color = "darkgray") + # Mean line
  #stat_summary(fun = mean, geom = "point", shape = 20, size = 3, color = "blue") + # Median line
  #stat_boxplot(geom = "errorbar", width = 0.25) + # Error lines
  scale_fill_manual(values = c("#AA57B7", "#281E73")) + # Customize colors
  labs(title = "Effects on positive emotions", x = "Condition", y = "Positive emotions (low to high)", fill = "Group") +
    scale_x_discrete(labels = custom_labels) + # Use custom labels for the x-axis
  theme_minimal()



# Add significance annotation using ggsignif
p_pos + geom_signif(comparisons = list(c("ingroup.Hypocrisy", "ingroup.Control")),
                map_signif_level = TRUE,
                y_position = max(long_emo$positive) + 1) +
    geom_signif(comparisons = list(c("outgroup.Hypocrisy", "outgroup.Control")),
                map_signif_level = TRUE,
                y_position = max(long_emo$positive) + 4 ) +
    geom_signif(comparisons = list(c("ingroup.Hypocrisy", "outgroup.Hypocrisy")),
                map_signif_level = TRUE,
                y_position = max(long_emo$positive) + 7) +
    geom_signif(comparisons = list(c("ingroup.Control", "outgroup.Control")), 
                map_signif_level = TRUE,
                y_position = max(long_emo$positive) + 10) 

```

### Likely (plausibility) analysis (likely ~ policy supp)

```{r likpol}
mlikely <- lmer(likely ~ cc_policy + (1|topic) , data = long_full)

mlikely2 <- lm(likely ~ cc_policy , data = long_full)

summary(mlikely)

ggplot(long_full, aes(x = cc_policy, y = likely)) +
  geom_smooth(method = "lm", se = TRUE, color = "blue", size = 1) +
  labs(x = "Climate Change Policy Support (cc_policy)",
       y = "Plausibility",
       title = "Plausibility by Climate Policy Support") +
  theme_minimal()


tab_model(mlikely,
  dv.labels = c("Plausibility"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")



## aggregated so that i see likely overall

#aggregating plausibility to see the overall effect of plausibility on donation  

agg2 <- aggregate(likely ~ id_w, data = long_full, FUN = mean) 

agg2$cc_policy <- tapply(long_full$cc_policy, long_full$id_w, unique)

mlikely2 <- lm(likely ~ cc_policy , data = agg2)

summary(mlikely2)

tab_model(mlikely2,
  dv.labels = c("Plausibility"),
  string.pred = "Coefficient",
  string.ci = "CI (95%)",
  string.p = "p",
  show.std = TRUE, 
  string.std_ci = "Std. CI",
  string.std = "Std. Beta",
  string.est = "Est.")



```