---
title: "Study1"
author: "Double Blind"
date: "2023-01-14"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

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(report)
library(jtools)
library(rstatix)
library(sjPlot)

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

## Loading data and preparing data

```{r load data}


data <- import(here("data", "study1_complete_151222.xlsx")) 
data <- data %>% 
  mutate(condition = as.factor(ifelse( is.na(c_flight_1), 'Hipocrisy', 'Control')))%>%
  mutate(attention_check = as.factor(attention_check))

data_long <- gather(data, topic, judgment, c_flight_1:h_attention_1, factor_key=TRUE)%>%
  filter(judgment!= 'NA')

data_long2 <- gather(data_long, env_topic, env_values, ev_needs_1:ev_attention_1, factor_key=TRUE) %>%
  filter(judgment!= 'NA')



data_ok1 <- data_long %>% 
  mutate(topic = gsub("c_", "", as.factor(topic))) %>%
  mutate(topic = gsub("h_", "", as.factor(topic))) %>% 
  mutate(topic = gsub("_1", "", as.factor(topic))) %>%
  filter(attention_check != 'NA')%>%
  filter(attention_check == 'Hamburguer')%>%
  filter(topic != 'attention')%>%
  filter(workerId != 'A1VAZ4O0J5QD7S')%>% #removing participants that did not complete the survey 
  filter(workerId != 'A1EDLMO8ACUTBP')%>%
  filter(workerId != 'A1DC5O5XEIMM1M')%>%
    rowwise()%>%
  mutate(pol_mean = mean(c(pol_electricity_1, pol_gasoline_1, pol_school_1, pol_tax_1)))%>%
  mutate_at('judgment', as.numeric)


data_okwide <- data_ok1%>%
  pivot_wider(names_from = topic, values_from = judgment)%>%
  rowwise()%>%
  mutate(judgment_m = (flight + recycle + beef+ car)/4)%>%
  rowwise()%>%
  mutate(pol_mean = mean(c(pol_electricity_1, pol_gasoline_1, pol_school_1, pol_tax_1)))%>%
  mutate_at('judgment_m', as.numeric) %>%
  mutate(pol_or_num = case_when(pol_orientation == "Other" ~ 0, pol_orientation == "Very Liberal/Left wing" ~ 1, pol_orientation == "Slightly Liberal/Left wing" ~ 2, pol_orientation == "Moderate" ~ 3, pol_orientation == "Slightly Conservative/Right wing" ~ 4, pol_orientation == "Very Conservative/Right wing" ~ 5 ))%>%
  mutate_at("pol_or_num", as.numeric)



data_mean <- data %>% 
  rowwise() %>% 
  dplyr::mutate(pol_mean = mean(c(pol_electricity_1, pol_gasoline_1, pol_school_1, pol_tax_1))) %>%
  dplyr::mutate(control_mean = mean(c(c_flight_1, c_recycle_1, c_beef_1, c_car_1))) %>%
  dplyr::mutate(hipocrisy_mean = mean(c(h_flight_1, h_recycle_1, h_beef_1, h_car_1))) %>%
  dplyr::mutate(all_mean = mean(c(na.omit(c_flight_1), na.omit(c_recycle_1), na.omit(c_beef_1), na.omit(c_car_1), na.omit(h_flight_1), na.omit(h_recycle_1), na.omit(h_beef_1), na.omit(h_car_1)))) 


```

## Demographics

```{r dems}

data_okwide %<>% mutate_at('age', as.numeric)%>%
  filter(age != 'NA')

age_mean <- mean(data_okwide$age)

age_mean

age_sd <- sd(data_okwide$age)

age_sd

data_pol<- data_okwide %>%
  select(pol_orientation)%>%
  group_by(pol_orientation)%>%
  dplyr::summarise(avg = (n()/256)*100)

data_pol

data_gen<- data_okwide %>%
  select(gender)%>%
  group_by(gender)%>%
  dplyr::summarise(avg = (n()/256)*100)

data_gen

data_edu<- data_okwide %>%
  select(education)%>%
  group_by(education)%>%
  dplyr::summarise(avg = (n()/256)*100)

data_edu


# political orientation distribution - quick look 
hist(data_okwide$pol_or_num)
```



## Hypothesis 1: Testing Hypocrisy penalty 


```{r meanjudg}

## First, lets check if we can average Judgment across topics 

data_c <- data %>%
    filter(attention_check != 'NA')%>%
  filter(attention_check == 'Hamburguer')%>%
  filter(workerId != 'A1VAZ4O0J5QD7S')%>% #removing participants that did not complete the survey 
  filter(workerId != 'A1EDLMO8ACUTBP')%>%
  filter(workerId != 'A1DC5O5XEIMM1M')%>%
  select(starts_with('c_'), -starts_with('c_attention'))%>%
  na.omit()

alpha(data_c)

data_h <- data %>%
     filter(attention_check != 'NA')%>%
  filter(attention_check == 'Hamburguer')%>%
  filter(workerId != 'A1VAZ4O0J5QD7S')%>% #removing participants that did not complete the survey 
  filter(workerId != 'A1EDLMO8ACUTBP')%>%
  filter(workerId != 'A1DC5O5XEIMM1M')%>%
  select(starts_with('h_'), -starts_with('h_attention'))%>%
  na.omit()

alpha(data_h)

#alphas are .78 and .85 (0.82 on average), which is acceptable for combining through topics/vignettes

## Pre-registered analysis: 

mh1_ok <- lm(judgment_m ~ condition, data= data_okwide)

summary(mh1_ok)

qqnorm(residuals(mh1_ok))

# the model does not follow normality of residuals, so since it is just a mean comparison, we use an appropriate test for our data: 


# double check there is no missing data before running wilcox test 

data_okwide2 <- data_okwide %>% 
  select(condition, judgment_m)%>%
  filter(!is.na(judgment_m) & !is.na(condition))


## wilcoxon test 

test_result <- wilcox.test(judgment_m ~ condition, data= data_okwide2)

# obtaining CIs 

wilcox.test(judgment_m ~ condition, data= data_okwide2, 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 <- sum(data_okwide2$condition == "Control")  # Number of observations in the Control group
n2 <- sum(data_okwide2$condition == "Hipocrisy")  # 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


## let's get the means and sd of Mena judgment by condition

means_sd_hc <- data_okwide2 %>%
  group_by(condition) %>%
  dplyr::summarise(mean = mean(judgment_m), sd = sd(judgment_m))


```

## Hypothesis 2: interaciton with CC mitigation policy support 

```{r h2re}

# let's check alpha for policy support across items

ccp_supp <- data_okwide %>%
  select(pol_electricity_1, pol_gasoline_1, pol_school_1, pol_tax_1)

alpha(ccp_supp) #alpha of 0.85, very good 

# quick look at distribution: 
  
  hist(data_okwide$pol_mean)

# model as pre-registered (no random effects)

mh2_ok<- lm(judgment_m~condition*pol_mean, data= data_okwide)

summary(mh2_ok)

eta_2 <- effectsize::eta_squared(mh2_ok)

m_eta_2 <- eta_2$Eta2_partial[3]

d_2 <- sqrt((4*m_eta_2)/(1 - m_eta_2))

qqnorm(residuals(mh2_ok)) # normality of residuals does not look great

## Let's try same model but random effects for vignettes' topics

mh2_pol <- lmer(judgment ~ condition*pol_mean + (1|topic), data = data_ok1)

summary(mh2_pol)

qqnorm(residuals(mh2_pol)) #normality of residuals improves in the mixed effects model, we keep this for reporting in the paper

```


### H2: effect size and tables

```{r tables}


m2_eta2 <- effectsize::eta_squared(mh2_pol)
m2_c_eta2 <- m2_eta2$Eta2_partial[1]
m2_p_eta2 <- m2_eta2$Eta2_partial[2]
m2_int_eta2 <- m2_eta2$Eta2_partial[3]

d_c_m2 <- sqrt((4*m2_c_eta2)/(1 - m2_c_eta2))
d_p_m2 <- sqrt((4*m2_p_eta2)/(1 - m2_p_eta2))
d_int_m2 <- sqrt((4*m2_int_eta2)/(1 - m2_int_eta2))

# interaction term effect size


d_int_m2 

# tables

### random effect model 

tab_model(mh2_pol,   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.")



### model without random effects

tab_model(mh2_ok,   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.")

```

### Figure


```{r figureh2}

okabe_ito_colors <- c( "#E69F00", "maroon3")
okabe_ito_colors_1 <- c("springgreen4","lightgreen")
okabe_ito_colors_2 <- c("maroon3", "pink2")

# Define CMYK-safe colors
okabe_ito_colors_cmyk <- c(rgb(0.9, 0.62, 0.0), rgb(0.8, 0.2, 0.47))

data_okwide%>%
  mutate_at('condition', as.factor)%>%
ggplot( aes(x = pol_mean, y = judgment_m, color = condition)) +
  stat_summary_bin(fun = "mean", bins = 7, geom = "point", size = 3) +  # Binned points
  geom_smooth(method = "loess", span = 0.5, se = TRUE) + # Loess smooth line with CI
  geom_jitter(alpha= 0.3)+
  scale_color_manual(values = okabe_ito_colors_cmyk, labels = c("Control", "Hypocrisy") ) +  # Color-blind friendly 
  xlim(1,7)+
  theme_minimal() +
  labs(title = "CC Hypocrisy Penalty", subtitle = "Judgment of Climate Hypocrisy vs Climate Harm, by levels of CC mitigation support",
       x = "CC Mitigation support (Low to High)",
       y = "Moral Judgment (Wrong to Right)", color= 'Condition') 


```