---
title: Analysis of Variance (ANOVA)
subtitle: for Each Zero-Sum Belief Items using Political Party Affiliation 
author:
  - name: Aashia Khan
    corresponding: true
    roles:
      - Conceptualization
      - Investigation
      - Methodology
      - Project administration
      - Resources
      - Supervision
    affiliations: 
      - Binghamton University
  - name: Zihan Hei
    corresponding: false
    orcid: 0009-0007-1617-2666 
    roles:
      - Software
      - Formal analysis
      - Visualization
      - Writing - original draft
    affiliations:
      - Binghamton University
  - name: Jeff John
    corresponding: false
    roles:
      - Methodology
      - Project administration
      - Writing - original draft
    affiliations:
      - Binghamton University
  - name: Shane McCarty
    orcid: 0000-0001-8930-7049
    corresponding: true
    roles:
      - Conceptualization
      - Investigation
      - Methodology
      - Project administration
      - Resources
      - Supervision
      - Writing - review & editing
    affiliations:
      - Binghamton University
      - Promote Care & Prevent Harm
keywords:
  - zero sum beliefs
  - social identities
  - political affiliation
  - racial identity
abstract: |
  *Background:* Zero-sum beliefs—the perception that one group's gains necessarily result in another group's losses—are important predictors of political attitudes. However, the referents for zero-sum beliefs as economic or social identity remain underexplored in relation to political ideology, party affiliation, and voting behavior in contemporary elections.
  
  *Method:* We conducted a comprehensive analysis examining three dimensions of zero-sum beliefs (general, economic, and social identity). Using Kruskal-Wallis tests on eleven zero-sum beliefs, we investigated how political party affiliation and racial/ethnic identity influenced endorsement of zero-sum beliefs across multiple domains. Subsequently, we examined whether these zero-sum belief patterns predicted self-reported voting for Donald Trump versus Kamala Harris in the 2024 presidential election.
  
  *Results:* Political party affiliation was a significant predictor for all eight zero-sum social identity beliefs, but none of the economic or general beliefs. Republican voters and certain racial/ethnic groups demonstrated higher endorsement of zero-sum social identity beliefs. A logistic regression shows that after controlling for political ideology, a composite of zero-sum social identity beliefs explains voting behavior in the 2024 presidential election, with stronger zero-sum social identity thinking associated with Trump support and lower zero-sum social identity beliefs predicting Harris support. Other sociodemographic factors and zero-sum economic thinking were not significant predictors.
  
  *Discussion:* Zero-sum social identity beliefs may represent a competitive core belief underlying contemporary political party affiliation and candidate preference. These findings affirm prior work that zero-sum thinking about economics differ from social identities, with similar levels of agreement on zero-sum economic beliefs across political parties but significantly different levels of agreement on zero-sum social identity beliefs by party affiliation. To the best of our knowledge, this study is the first to show that zero-sum thinking about social identities predicts voter preference in the 2024 election. Ultimately, future work needs to examine how to reduce zero-sum social identity thinking. 

plain-language-summary: |
  Many people hold "zero-sum beliefs"—the idea that when one group succeeds, another group must fail. This study examined how these beliefs relate to politics and voting in the 2024 presidential election. We found that while people across political parties hold similar zero-sum beliefs about economic issues, Republicans were much more likely than Democrats to hold zero-sum beliefs about social identities (such as undocumented vs. citizens, White vs. Black, LGBTQ vs. Religious). These social identity zero-sum beliefs were strong predictors of self-reported voting behavior: people with stronger zero-sum social identity beliefs were more likely to vote for Donald Trump, while those with weaker zero-sum social identity beliefs were more likely to vote for Kamala Harris. Post election discourse suggests the economy -- and not social identity -- was the primary driver of voters. Yet, our results show that zero-sum beliefs about economics didn't predict candidate preference. But, the social identity zero-sum beliefs were the second most important factor after accounting for political ideology as the primary predictor. This suggests that how people think about competition between different social groups is a key psychological factor driving political divisions and voting choices in America today.
date: last-modified
number-sections: true
format:
  html:
    toc: true
    code-fold: true
    code-summary: "Show the code"
    comments:
      hypothesis: true
  pdf:
    number-sections: true
    code-fold: true
    code-summary: "Show the code"
    keep-tex: true
---

# Introduction

## Political Realignment (SM)

## Zero-Sum Beliefs Regarding Game, Economic, and Social Identities (AK, JJ)

## Zero-Sum Beliefs and Political Party (AK)

## Zero-Sum Beliefs and Racial Identity (JJ)

## Zero-Sum Beliefs and Voter Behavior (AK)

# Methods

## Study Participants (AK)

More than fifty-five thousand people who are active on the Prolific platform were eligible to complete a 45-minute health beliefs survey with measures on various beliefs associated with politics and health. Using a quota sample by sex (50% male/ 50% female), political affiliation (33% Republican, 33% Democrat, and 34% Independent), and race/ethnicity (White 40%, Black 20%, Asian 20%, Mixed 10%, and Other 10%), Prolific recruited one hundred and twenty-five people to complete the survey.

## Measures (AK)

This needs updating ...

Respondents were asked to report their gender (selections: Girl or woman, boy or man, nonbinary/genderfluid/genderqueer, I am not sure/questioning), race (selections: American Indian or Alaska Native, Asian, Black or African American, Hispanic or Latine, Middle Eastern or North African, Native Hawaiian/Pacific Islander, White), social status (selections: Lowest 1- Highest 10), and political beliefs (selections: far left/leftist, very liberal, liberal, moderate, conservative, very conservative, alt-right/far-right). Respondents were asked to report their level of agreement using a likert scale (1: Strongly Disbelieve, 2: Disbelieve, 3: Somewhat Disbelieve, 4: Neither, 5: Somewhat Believe, 6: Believe, 7: Strongly Believe) with X number of statements.
- political affiliation/party
- political ideology/beliefs
- gender
- education
- social class
- racial identity
- 2024 voting behavior

## Analysis Plan (ZH)

Two-way Analysis of Variance (ANOVA) were run on 11 zero-sum beliefs ( `ZEROSUM_` ) – including X economic zero-sum and X social identity zero-sum beliefs – using the factors of political affiliation ( `POLITICALPARTY` ) and racial/ethnic identity ( `RACIALIDENTITY.4` ).

An explanatory logistic regression was conducted to classy voter's preference for Donald Trump (1) vs. Kamala Harris (0) ( `TRUMPVOTE` ) using `GENDER_MALE`, `RACE_WHITE`, `EDUCATION_HIGH`, `SOCIALCLASS`, `POLITICALBELIEFS`, `ZEROSUM_IDENTITY`, `ZEROSUM_ECONOMIC`, and `ZEROSUM_1`. 

ZH: Additionally, we used `tidymodels` package in R to classify using a predictive model. Using an integrative modeling approach, we provide an explanatory and predictive model for classifying voter preferences. 

# Results

```{r}
#| label: load library
#| echo: false
#| output: false
#| warning: false
#| error: false
#| results: false

# Load necessary libraries
library(ggplot2)
library(plotly)
library(dplyr)
library(tidyr)
library(ggrain)  
library(rmarkdown)
library(readr)
library(patchwork)
library(codebookr)
library(dplyr, warn.conflicts = FALSE)
library(haven)
library(codebook)
library(rempsyc)
library(car)
library(knitr)
library(broom)
library(ggdist)
library(patchwork)
library(car)
library(dataMaid)
library(devtools)
library(MBESS)
library(apaTables)
library(ggpubr)
library(psych)
library(forcats)
library(corrplot)
```


```{r}
#| label: import select_data
select_data <- read.csv("/cloud/project/data/select_data.csv")
```

```{r}
#| label: table ETHNICITY

table(select_data$ETHNICITY)
## notice that the groups are unequal which is a problem for the ANOVA.
## TO DO: consider whether we should do a non-parametric approach

# SOURCE FOR DUMMY CODE
## https://stats.oarc.ucla.edu/other/mult-pkg/faq/general/faqwhat-is-dummy-coding/
## https://www.nathanwhudson.com/courses/methods/resources/10.%20Dummy%20Coding.pdf
## can you find an R specific example? 
```



```{r}
#| label: party colors for plots

party_colors <- c("Democrat" = "#0015BC", "Republican" = "#E9141D", "Independent" = "#732C7B")
```

## Factors Associated with Zero-Sum Beliefs

### Zero Sum 1

```{r}
#| label: groupstats for ZERO_SUM1 

library(dplyr)

## We group by both political party and racial identity to explore whether attitudes about zero-sum beliefs differ across these two demographic factors.
## Then, we summarize the data using dplyr's `summarise()` to calculate:
## - mean of ZEROSUM_1 (agreement with zero-sum beliefs)
## - standard error (SE) of the mean, which reflects the variability in responses

group_stats_1 <- select_data %>%
  group_by(POLITICALPARTY, RACIALIDENTITY.4) %>%
  summarise(
    mean_ZEROSUM_1 = mean(ZEROSUM_1, na.rm = TRUE),  # average score for each group
    se = sd(ZEROSUM_1, na.rm = TRUE) / sqrt(n()),    # standard error for error bars
    .groups = 'drop'   # ungroup for downstream operations
  )

# Quick check of the output's structure and first few rows
str(group_stats_1)
head(group_stats_1)

```

```{r}
#| label: zerosum1 pointrange

## This plot visualizes how different political and racial groups scored on the ZEROSUM_1 scale
## We use a point-range plot to show:
##  - the group mean (point)
##  - variability (standard error bars above and below the mean)

## The color of each point reflects the racial identity of respondents

zesty_four <- c("#E69F00", "#009E73", "#999999", "#CC79A7")  # custom color palette

plot.zerosum1 <- ggplot(group_stats_1, aes(POLITICALPARTY, mean_ZEROSUM_1)) +
  geom_pointrange(aes(
      color = RACIALIDENTITY.4,
      ymin = mean_ZEROSUM_1 - se,  # lower bound of SE
      ymax = mean_ZEROSUM_1 + se   # upper bound of SE
    ),
    position = position_dodge(0.3)) +
  theme_bw() +
  labs(
    title = "Life is so devised that when somebody gains, others have to lose.",
    x = "Political Affiliation",
    y = "Agreement with Zero Sum Beliefs",
    color = "RACIALIDENTITY.4"
  ) +
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8)
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", "2: Disbelieve", "3: Somewhat Disbelieve",
               "4: Neither", "5: Somewhat Believe", "6: Believe", "7: Strongly Believe"),
    limits = c(1, 7)
  ) +
  scale_color_manual(values = zesty_four)

# Print and save the plot for documentation
print(plot.zerosum1)
ggsave("plot1:gain-lose.png", plot = plot.zerosum1, width = 10, height = 8, dpi = 300)

```

```{r}
#| label: ZEROSUM_1 average score with 95% CI

## This chunk computes the 95% confidence intervals (CIs) for the mean ZEROSUM_1 score by political party.
## Unlike SE, a CI gives a sense of the range in which the true population mean likely falls.

group_stats_1_avg <- select_data %>%
  group_by(POLITICALPARTY) %>%
  summarise(
    mean_ZEROSUM_1 = mean(ZEROSUM_1, na.rm = TRUE),
    se = sd(ZEROSUM_1, na.rm = TRUE) / sqrt(n()),  # standard error
    n = n(),
    .groups = 'drop'
  ) %>%
  mutate(
    # Calculate t-critical value for 95% CI using the t-distribution
    t_value = qt(0.975, df = n - 1),
    ci_lower = mean_ZEROSUM_1 - t_value * se,
    ci_upper = mean_ZEROSUM_1 + t_value * se
  )

## These statistics will be used for overlaying on the raincloud plot later

```



```{r}
#| label: plot.zerosum1_enhanced and Raincloud plot

## This chunk creates a "raincloud plot" to visualize the distribution of ZEROSUM_1 scores by political party. 
## Raincloud plots combine density shapes, raw data points, and summary stats (mean + CI).

plot.zerosum1_raincloud <- select_data %>%
  filter(!is.na(ZEROSUM_1), !is.na(POLITICALPARTY)) %>%
  ggplot(aes(x = POLITICALPARTY, y = ZEROSUM_1, fill = POLITICALPARTY, color = POLITICALPARTY)) +

## Half-eye density shows the distribution shape (like a smoothed histogram).
  stat_halfeye(
    adjust = 0.5,         
    width = 0.6,          
    .width = 0,           
    justification = -0.2, 
    point_colour = NA,
    alpha = 0.7           
  ) +

## Add 95% confidence intervals from earlier calculations
  geom_pointrange(
    data = group_stats_1_avg,
    aes(
      x = POLITICALPARTY,
      y = mean_ZEROSUM_1,
      ymin = ci_lower,
      ymax = ci_upper
    ),
    color = "black",      
    size = 0.9,
    shape = 21,
    fill = "white",
    stroke = 1.2,
    inherit.aes = FALSE   
  ) +

## Add jittered raw data points (each dot = individual response)
  geom_point(
    aes(shape = POLITICALPARTY),
    size = 1.5,
    alpha = 0.5,
    position = position_jitter(
      seed = 1, width = 0.1  # keeps jitter reproducible
    )
  ) +

## Overlay group means as white dots
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3,
    color = "white",
    stroke = 1.5,
    shape = 21
  ) +

## Color palette and theming
  scale_fill_manual(values = party_colors) +
  scale_color_manual(values = party_colors) +
  theme_bw() +
  labs(
    title = "Life is so devised that when somebody gains, others have to lose.",
    x = "Political Affiliation",
    y = "Level of Agreement",
    caption = "White dots = means; colored dots = individual responses"
  ) +
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    plot.caption = element_text(size = 9, color = "gray40", hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    legend.position = "none",
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3)
  ) +
  
## Match axis scale to Likert-style labels for interpretation
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7)
  )

## Show the final plot
plot.zerosum1_raincloud

```





## Inconsistent ANOVA results

aov() function is used for NO interactions, but we have an interaction of POLITICALPARTY \* RACIALIDENTITY.4therefore, we should use the newer package within the linear model lm function using type = 2.

THIS IS THE CODE BLOCK TO USE THROUGHOUT THE REPORT

```{r}
#| label: ZEROSUM_1 ANOVA II

## two way ANOVA using  
library(car)
model.ZEROSUM_1.PE <- lm(ZEROSUM_1 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)
Anova(model.ZEROSUM_1.PE, type = 2)

## There is no statistically significant difference in ZEROSUM_1 scores between political parties, ethnicities, or their interaction.
```

```{r}
#| label: APA format ZEROSUM_1 ANOVA II

# Two-way ANOVA
model.ZEROSUM_1.PE <- lm(ZEROSUM_1 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)

anova_df <- as.data.frame(Anova(model.ZEROSUM_1.PE, type = 2))%>%
  tibble::rownames_to_column(var = "Predictor")
names(anova_df)[names(anova_df) == "Pr(>F)"] <- "p"

# APA format
nice_table(anova_df, 
           title = c("Table 1", "ZeroSum_1 ANOVA Table"), 
           highlight = 0.05, 
           stars = TRUE)

## There is no statistically significant difference in ZEROSUM_1 scores between political parties, ethnicities, or their interaction.
```

```{r}
#| label: ZEROSUM_1 ANOVA plot

plot(model.ZEROSUM_1.PE)  # diagnostic plots
shapiro.test(residuals(model.ZEROSUM_1.PE))  # normality test
```

```{r}
#| label: model.ZEROSUM_1.PE on residuals

# Check how bad the non-normality is
hist(residuals(model.ZEROSUM_1.PE))
# check: for normal distribution
# results: not normal, but it's close (and likely would be with larger sample size)
# interpretation: this is okay

qqnorm(residuals(model.ZEROSUM_1.PE))
# check: points should fall near the diagonal line 
# results: potential outliers at the extremes

qqline(residuals(model.ZEROSUM_1.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation
```

```{r}
#| label: model.ZEROSUM_1.PE boxplot to check for outliers

boxplot(residuals(model.ZEROSUM_1.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

## overall: the visual inspection across the previous plots suggests we are okay to move forward with the ANOVA test, but must report the shapiro.test results here. 
```




### Zero Sum 2

```{r}
#| label: groupstats for ZERO_SUM2 

library(dplyr)
group_stats_2 <- select_data %>%
  group_by(POLITICALPARTY, RACIALIDENTITY.4) %>%
  summarise(
    mean_ZEROSUM_2 = mean(ZEROSUM_2, na.rm = TRUE),
    se = sd(ZEROSUM_2, na.rm = TRUE) / sqrt(n()),
    .groups = 'drop'
  )

# Check the structure
str(group_stats_2)
head(group_stats_2)
```


```{r}
#| label: zerosum2 pointrange

# Define the zesty_four palette
zesty_four <- c("#F5793A", "#A95AA1", "#85C0F9", "#0F2080") 

# Create the plot with pointrange
plot.zerosum2 <- ggplot(group_stats_2, 
                        aes(POLITICALPARTY, mean_ZEROSUM_2)) +
  geom_pointrange(aes(color = RACIALIDENTITY.4,
                      ymin = mean_ZEROSUM_2 - se,
                      # Use standard error instead of sd
                      ymax = mean_ZEROSUM_2 + se),
                      # Use standard error instead of sd
                      position = position_dodge(0.3)) +
  theme_bw() +
  labs(title = "When some people are getting poorer, \n it means that other people are getting richer.",
       x = "Political Affiliation",
       y = "Agreement with Zero Sum Beliefs",
       color = "RACIALIDENTITY.4") + 
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8)
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7)  # Set the y-axis limits to match the scale
  )+
  scale_color_manual(values = zesty_four)

# Print and save the plots
print(plot.zerosum2)
ggsave("plot2:poor-rich.png", 
       plot = plot.zerosum2, 
       width = 10, 
       height = 8, 
       dpi = 300)
```

```{r}
#| label: ZEROSUM_2 average score with 95% CI

# ZEROSUM_1 average score with 95% CI
group_stats_2_avg <- select_data %>%
  group_by(POLITICALPARTY) %>%
  summarise(
    mean_ZEROSUM_2 = mean(ZEROSUM_2, na.rm = TRUE),
    se = sd(ZEROSUM_2, na.rm = TRUE) / sqrt(n()),
    n = n(),
    .groups = 'drop'
  ) %>%
  mutate(
    # Calculate 95% confidence interval using t-distribution
    t_value = qt(0.975, df = n - 1),  # 95% CI, two-tailed
    ci_lower = mean_ZEROSUM_2 - t_value * se,
    ci_upper = mean_ZEROSUM_2 + t_value * se
  )
```



```{r}
#| label: plot.zerosum2_enhanced and Raincloud plot

library(ggdist)
library(patchwork)

plot.zerosum2_raincloud <- select_data %>%
  filter(!is.na(ZEROSUM_2), !is.na(POLITICALPARTY)) %>%
  ggplot(aes(x = POLITICALPARTY, y = ZEROSUM_2, fill = POLITICALPARTY, color = POLITICALPARTY)) +
  stat_halfeye(
    adjust = 0.5,
    width = 0.6,
    .width = 0,
    justification = -0.2,
    point_colour = NA,
    alpha = 0.7
  ) +
  geom_pointrange(
      data = group_stats_2_avg,
      aes(
        x = POLITICALPARTY,
        y = mean_ZEROSUM_2,
        ymin = ci_lower,
        ymax = ci_upper
      ),
      color = "black",
      size = 0.9,
      shape = 21,
      fill = "white",
      stroke = 1.2,
      inherit.aes = FALSE
    ) +
  geom_point(
    aes(shape = POLITICALPARTY),
    size = 1.5,
    alpha = 0.5,
    position = position_jitter(
      seed = 1, width = 0.1
    )
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3,
    color = "white",
    stroke = 1.5,
    shape = 21
  ) +
  scale_fill_manual(values = party_colors) +
  scale_color_manual(values = party_colors) +
  theme_bw() +
  labs(
    title = "When some people are getting poorer, \n it means that other people are getting richer.",
    x = "Political Affiliation",
    y = "Level of Agreement",
    caption = "White dots = means; colored dots = individual responses"
  ) +
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    plot.caption = element_text(size = 9, color = "gray40", hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    legend.position = "none",
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3)  # Fixed
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7) 
  )
plot.zerosum2_raincloud

```


```{r}
#| label: ZEROSUM_2 ANOVA

# Two-way ANOVA
model.ZEROSUM_2.PE <- lm(ZEROSUM_2 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)
Anova(model.ZEROSUM_2.PE)

## There is no statistically significant difference in ZEROSUM_2 scores between political parties, ethnicities, or their interaction.
```

```{r}
#| label: APA format ZEROSUM_2 ANOVA

# Two-way ANOVA
model.ZEROSUM_2.PE <- lm(ZEROSUM_2 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)

anova_df <- as.data.frame(Anova(model.ZEROSUM_2.PE, type = 2))%>%
  tibble::rownames_to_column(var = "Predictor")
names(anova_df)[names(anova_df) == "Pr(>F)"] <- "p"

# APA format
nice_table(anova_df, 
           title = c("Table 2", "ZeroSum_2 ANOVA Table"), 
           highlight = 0.05, 
           stars = TRUE)

## There is no statistically significant difference in ZEROSUM_2 scores between political parties, ethnicities, or their interaction.
```

```{r}
#| label: ZEROSUM_2 ANOVA plot

plot(model.ZEROSUM_2.PE)  # diagnostic plots
shapiro.test(residuals(model.ZEROSUM_2.PE))  # normality test
```

```{r}
#| label: model.ZEROSUM_2.PE on residuals

# Check how bad the non-normality is
hist(residuals(model.ZEROSUM_2.PE))
# check: for normal distribution
# results: approximately normal, but not perfectly normal
# interpretation: this is okay

qqnorm(residuals(model.ZEROSUM_2.PE))
# check: points should fall near the diagonal line 
# results: potential outliers at the extremes

qqline(residuals(model.ZEROSUM_2.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

```

```{r}
#| label: model.ZEROSUM_2.PE boxplot to check for outliers
 
boxplot(residuals(model.ZEROSUM_2.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

## overall: the visual inspection across the previous plots suggests we are okay to move forward with the ANOVA test, but must report the shapiro.test results here. 

```

### \*Zero Sum 3 (recheck the outliers)


```{r}
#| label: groupstats for ZERO_SUM3 

library(dplyr)
group_stats_3 <- select_data %>%
  group_by(POLITICALPARTY, RACIALIDENTITY.4) %>%
  summarise(
    mean_ZEROSUM_3 = mean(ZEROSUM_3, na.rm = TRUE),
    se = sd(ZEROSUM_3, na.rm = TRUE) / sqrt(n()),
    .groups = 'drop'
  )

# Check the structure
str(group_stats_3)
head(group_stats_3)
```

```{r}
#| label: zerosum3 pointrange

# Define the zesty_four palette
zesty_four <- c("#F5793A", "#A95AA1", "#85C0F9", "#0F2080") 

# Create the plot with pointrange
plot.zerosum3 <- ggplot(group_stats_3,
                        aes(POLITICALPARTY, mean_ZEROSUM_3)) +
  geom_pointrange(aes(color = RACIALIDENTITY.4,
                      ymin = mean_ZEROSUM_3 - se,
                      # Use standard error instead of sd
                      ymax = mean_ZEROSUM_3 + se),
                      # Use standard error instead of sd
                      position = position_dodge(0.3)) +
  theme_bw() +
  labs(title = "The wealth of a few is acquired at the expense of many.",
       x = "Political Affiliation",
       y = "Agreement with Zero Sum Beliefs",
       color = "RACIALIDENTITY.4") + 
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8)
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7)  # Set the y-axis limits to match the scale
  )+
  scale_color_manual(values = zesty_four)

# Print and save the plots
print(plot.zerosum3)
ggsave("plot3:wealthfew-many.png", 
       plot = plot.zerosum3, 
       width = 10, 
       height = 8, 
       dpi = 300)

```


```{r}
#| label: ZEROSUM_3 average score with 95% CI

# ZEROSUM_3 average score with 95% CI
group_stats_3_avg <- select_data %>%
  group_by(POLITICALPARTY) %>%
  summarise(
    mean_ZEROSUM_3 = mean(ZEROSUM_3, na.rm = TRUE),
    se = sd(ZEROSUM_3, na.rm = TRUE) / sqrt(n()),
    n = n(),
    .groups = 'drop'
  ) %>%
  mutate(
    # Calculate 95% confidence interval using t-distribution
    t_value = qt(0.975, df = n - 1),  # 95% CI, two-tailed
    ci_lower = mean_ZEROSUM_3 - t_value * se,
    ci_upper = mean_ZEROSUM_3 + t_value * se
  )
```



```{r}
#| label: plot.zerosum3_enhanced and Raincloud plot

plot.zerosum3_raincloud <- select_data %>%
  filter(!is.na(ZEROSUM_3), !is.na(POLITICALPARTY)) %>%
  ggplot(aes(x = POLITICALPARTY, y = ZEROSUM_3, fill = POLITICALPARTY, color = POLITICALPARTY)) +
  stat_halfeye(
    adjust = 0.5,
    width = 0.6,
    .width = 0,
    justification = -0.2,
    point_colour = NA,
    alpha = 0.7
  ) +
  geom_pointrange(
      data = group_stats_3_avg,
      aes(
        x = POLITICALPARTY,
        y = mean_ZEROSUM_3,
        ymin = ci_lower,
        ymax = ci_upper
      ),
      color = "black",
      size = 0.9,
      shape = 21,
      fill = "white",
      stroke = 1.2,
      inherit.aes = FALSE
    ) +
  geom_point(
    aes(shape = POLITICALPARTY),
    size = 1.5,
    alpha = 0.5,
    position = position_jitter(
      seed = 1, width = 0.1
    )
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3,
    color = "white",
    stroke = 1.5,
    shape = 21
  ) +
  scale_fill_manual(values = party_colors) +
  scale_color_manual(values = party_colors) +
  theme_bw() +
  labs(
    title = "The wealth of a few is acquired at the expense of many.",
    x = "Political Affiliation",
    y = "Level of Agreement",
    caption = "White dots = means; colored dots = individual responses"
  ) +
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    plot.caption = element_text(size = 9, color = "gray40", hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    legend.position = "none",
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3)  # Fixed
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7) 
  )
plot.zerosum3_raincloud

```


```{r}
#| label: ZEROSUM_3 ANOVA

# Two-way ANOVA
model.ZEROSUM_3.PE <- lm(ZEROSUM_3 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)

Anova(model.ZEROSUM_3.PE, type = 2)

## There is a statistically significant difference in ZEROSUM_3 scores between ethnicities.

```

```{r}
#| label:  APA format ZEROSUM_3 ANOVA

# Two-way ANOVA
## ERROR here I think 
### I commented it out so that Quarto would allow me to publish this again. 
### You can uncomment and fix whenever now...
#model.ZEROSUM_3.PE <- lm(ZEROSUM_3 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)

#anova_df <- as.data.frame(Anova(model.ZEROSUM_3.PE, type = 2))%>%
  #tibble::rownames_to_column(var = "Predictor")
#names(anova_df)[names(anova_df) == "Pr(>F)"] <- "p"

# APA format
#nice_table(anova_df, 
           #title = c("Table 3", "ZeroSum_3 ANOVA Table"), 
           #highlight = 0.05, 
           #stars = TRUE)

## There is a statistically significant difference in ZEROSUM_3 scores between ethnicities.

```

```{r}
#| label: model.ZEROSUM_3.PE Tukey HSD test
#| include: false
#| eval: false
# Perform Tukey HSD test
#tukey_result <- TukeyHSD(model.ZEROSUM_3.PE)

# Print the Tukey HSD results
#print(tukey_result)

```

```{r}
#| label: ZEROSUM_3 ANOVA plot

plot(model.ZEROSUM_3.PE)  # diagnostic plots
shapiro.test(residuals(model.ZEROSUM_3.PE))  # normality test
```

```{r}
#| label: model.ZEROSUM_3.PE on residuals

# Check how bad the non-normality is
hist(residuals(model.ZEROSUM_3.PE))
# check: for normal distribution
# results: not normal, but it's close (and likely would be with larger sample size)
# interpretation: this is okay

qqnorm(residuals(model.ZEROSUM_3.PE))
# check: points should fall near the diagonal line 
# results: potential outliers at the extremes

qqline(residuals(model.ZEROSUM_3.PE))
# check: for outliers
# results: have outliers
# interpretation: we can drop the outliers to improve interpretation

```

```{r}
#| label: model.ZEROSUM_3.PE boxplot to check for outliers

boxplot(residuals(model.ZEROSUM_3.PE))
# check: for outliers
# results: have outliers
# interpretation: we can drop the outliers to improve interpretation
```

```{r}
#| label: model.ZEROSUM_3.PE.clean drop outliers

resid_vals <- residuals(model.ZEROSUM_3.PE)

# Identify outlier values in the residuals using the IQR method from boxplot.stats
outlier_values <- boxplot.stats(resid_vals)$out
non_outlier_index <- !(resid_vals %in% outlier_values)

# Filter dataset
alldata_no_outliers <- select_data[non_outlier_index, ]

# Re-run ANOVA with cleaned data
model.ZEROSUM_3.PE.clean <- aov(ZEROSUM_3 ~ POLITICALPARTY * RACIALIDENTITY.4, data = alldata_no_outliers)
summary(model.ZEROSUM_3.PE.clean)

```

```{r}
#| label: model.ZEROSUM_3.PE.clean ANOVA plot

plot(model.ZEROSUM_3.PE.clean)  # diagnostic plots
shapiro.test(residuals(model.ZEROSUM_3.PE.clean))  # normality test

```

```{r}
#| label: model.ZEROSUM_3.PE.clean on residuals

# Check how bad the non-normality is
hist(residuals(model.ZEROSUM_3.PE.clean))
# check: for normal distribution
# results: not normal, but it's close (and likely would be with larger sample size)
# interpretation: this is okay

qqnorm(residuals(model.ZEROSUM_3.PE.clean))
# check: points should fall near the diagonal line 
# results: potential outliers at the extremes

qqline(residuals(model.ZEROSUM_3.PE.clean))
# check: for outliers
# results: still have outliers
# interpretation: 
```

```{r}
#| label: model.ZEROSUM_3.PE.clean boxplot to check for outliers

boxplot(residuals(model.ZEROSUM_3.PE.clean))
# check: for outliers
# results: have outliers
# interpretation: 

```

### Zero Sum 4



```{r}
#| label: groupstats for ZERO_SUM4 

group_stats_4 <- select_data %>%
  group_by(POLITICALPARTY, RACIALIDENTITY.4) %>%
  summarise(
    mean_ZEROSUM_4 = mean(ZEROSUM_4, na.rm = TRUE),
    se = sd(ZEROSUM_4, na.rm = TRUE) / sqrt(n()),
    .groups = 'drop'
  )

#check the structure
str(group_stats_4)
head(group_stats_4)
```


```{r}
#| label: zerosum4 pointrange

# Create the plot with pointrange
plot.zerosum4 <- ggplot(group_stats_4, aes(POLITICALPARTY,
                                           mean_ZEROSUM_4)) +
  geom_pointrange(aes(color = RACIALIDENTITY.4, 
                      ymin = mean_ZEROSUM_4 - se,
                      # Use standard error instead of sd
                      ymax = mean_ZEROSUM_4 + se),
                      # Use standard error instead of sd
                      position = position_dodge(0.3)) +
  theme_bw() +
  labs(title = "As women face less sexism, men end up facing more sexism.",
       x = "Political Affiliation",
       y = "Agreement with Zero Sum Beliefs",
       color = "RACIALIDENTITY.4") + 
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8)
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7)  # Set the y-axis limits to match the scale
  )+
  scale_color_manual(values = zesty_four)

# Print and save the plots
print(plot.zerosum4)
ggsave("plot4:women-men.png", 
       plot = plot.zerosum4, 
       width = 10, 
       height = 8, 
       dpi = 300)


```


```{r}
#| label: ZEROSUM_4 average score with 95% CI

# ZEROSUM_4 average score with 95% CI
group_stats_4_avg <- select_data %>%
  group_by(POLITICALPARTY) %>%
  summarise(
    mean_ZEROSUM_4 = mean(ZEROSUM_4, na.rm = TRUE),
    se = sd(ZEROSUM_4, na.rm = TRUE) / sqrt(n()),
    n = n(),
    .groups = 'drop'
  ) %>%
  mutate(
    # Calculate 95% confidence interval using t-distribution
    t_value = qt(0.975, df = n - 1),  # 95% CI, two-tailed
    ci_lower = mean_ZEROSUM_4 - t_value * se,
    ci_upper = mean_ZEROSUM_4 + t_value * se
  )
```



```{r}
#| label: plot.zerosum4_enhanced and Raincloud plot

plot.zerosum4_raincloud <- select_data %>%
  filter(!is.na(ZEROSUM_4), !is.na(POLITICALPARTY)) %>%
  ggplot(aes(x = POLITICALPARTY, y = ZEROSUM_4, fill = POLITICALPARTY, color = POLITICALPARTY)) +
  stat_halfeye(
    adjust = 0.5,
    width = 0.6,
    .width = 0,
    justification = -0.2,
    point_colour = NA,
    alpha = 0.7
  ) +
  geom_pointrange(
      data = group_stats_4_avg,
      aes(
        x = POLITICALPARTY,
        y = mean_ZEROSUM_4,
        ymin = ci_lower,
        ymax = ci_upper
      ),
      color = "black",
      size = 0.9,
      shape = 21,
      fill = "white",
      stroke = 1.2,
      inherit.aes = FALSE
    ) +
  geom_point(
    aes(shape = POLITICALPARTY),
    size = 1.5,
    alpha = 0.5,
    position = position_jitter(
      seed = 1, width = 0.1
    )
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3,
    color = "white",
    stroke = 1.5,
    shape = 21
  ) +
  scale_fill_manual(values = party_colors) +
  scale_color_manual(values = party_colors) +
  theme_bw() +
  labs(
    title = "As women face less sexism, men end up facing more sexism.",
    x = "Political Affiliation",
    y = "Level of Agreement",
    caption = "White dots = means; colored dots = individual responses"
  ) +
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    plot.caption = element_text(size = 9, color = "gray40", hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    legend.position = "none",
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3)  # Fixed
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7) 
  )
plot.zerosum4_raincloud

```


```{r}
#| label: ZEROSUM_4 ANOVA

# Two-way ANOVA
model.ZEROSUM_4.PE <- lm(ZEROSUM_4 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)
Anova(model.ZEROSUM_4.PE)

## There is a statistically significant difference in ZEROSUM_4 scores between political parties.

```

```{r}
#| label: APA format ZEROSUM_4 ANOVA

# Two-way ANOVA
model.ZEROSUM_4.PE <- lm(ZEROSUM_4 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)

anova_df <- as.data.frame(Anova(model.ZEROSUM_4.PE, type = 2))%>%
  tibble::rownames_to_column(var = "Predictor")
names(anova_df)[names(anova_df) == "Pr(>F)"] <- "p"

# APA format
nice_table(anova_df, 
           title = c("Table 4", "ZeroSum_4 ANOVA Table"), 
           highlight = 0.05, 
           stars = TRUE)

## There is a statistically significant difference in ZEROSUM_4 scores between political parties.

```

```{r}
#| label: model.ZEROSUM_4.PE Tukey HSD test
#| include: false
#| eval: false
# Perform Tukey HSD test
tukey_result <- TukeyHSD(model.ZEROSUM_4.PE)

# Print the Tukey HSD results
print(tukey_result)

```

```{r}
#| label: ZEROSUM_4 ANOVA plot

plot(model.ZEROSUM_4.PE)  # diagnostic plots
shapiro.test(residuals(model.ZEROSUM_4.PE))  # normality test

```

```{r}
#| label: model.ZEROSUM_4.PE on residuals

# Check how bad the non-normality is
hist(residuals(model.ZEROSUM_4.PE))
# check: for normal distribution
# results: not normal, but it's close (and likely would be with larger sample size)
# interpretation: this is okay

qqnorm(residuals(model.ZEROSUM_4.PE))
# check: points should fall near the diagonal line 
# results: potential outliers at the extremes

qqline(residuals(model.ZEROSUM_4.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

```

```{r}
#| label: model.ZEROSUM_4.PE boxplot to check for outliers

boxplot(residuals(model.ZEROSUM_4.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

## overall: the visual inspection across the previous plots suggests we are okay to move forward with the ANOVA test, but must report the shapiro.test results here. 

```

### Zero Sum 5

```{r}
#| label: groupstats for ZERO_SUM5 

library(dplyr)
group_stats_5 <- select_data %>%
  group_by(POLITICALPARTY, RACIALIDENTITY.4) %>%
  summarise(
    mean_ZEROSUM_5 = mean(ZEROSUM_5, na.rm = TRUE),
    se = sd(ZEROSUM_5, na.rm = TRUE) / sqrt(n()),
    .groups = 'drop'
  )

# Check the structure
str(group_stats_5)
head(group_stats_5)
```

```{r}
#| label: zerosum5 pointrange

# Define the zesty_four palette
zesty_four <- c("#F5793A", "#A95AA1", "#85C0F9", "#0F2080") 

# Create the plot with pointrange
plot.zerosum5 <- ggplot(group_stats_5, aes(POLITICALPARTY, mean_ZEROSUM_5)) +
  geom_pointrange(aes(color = RACIALIDENTITY.4, 
                      ymin = mean_ZEROSUM_5 - se,  # Use standard error instead of sd
                      ymax = mean_ZEROSUM_5 + se), # Use standard error instead of sd
                  position = position_dodge(0.3)) +
  theme_bw() +
  labs(title = "Less discrimination against minorities means more \n discrimination against whites.",
       x = "Political Affiliation",
       y = "Agreement with Zero Sum Beliefs",
       color = "RACIALIDENTITY.4") + 
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8)
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7)  # Set the y-axis limits to match the scale
  )+
  scale_color_manual(values = zesty_four)

# Print and save the plots
print(plot.zerosum5)
ggsave("plot5:minorities-whites.png", 
       plot = plot.zerosum5, 
       width = 10, 
       height = 8, 
       dpi = 300)

```

```{r}
#| label: ZEROSUM_5 average score with 95% CI

# ZEROSUM_5 average score with 95% CI
group_stats_5_avg <- select_data %>%
  group_by(POLITICALPARTY) %>%
  summarise(
    mean_ZEROSUM_5 = mean(ZEROSUM_5, na.rm = TRUE),
    se = sd(ZEROSUM_5, na.rm = TRUE) / sqrt(n()),
    n = n(),
    .groups = 'drop'
  ) %>%
  mutate(
    # Calculate 95% confidence interval using t-distribution
    t_value = qt(0.975, df = n - 1),  # 95% CI, two-tailed
    ci_lower = mean_ZEROSUM_5 - t_value * se,
    ci_upper = mean_ZEROSUM_5 + t_value * se
  )
```



```{r}
#| label: plot.zerosum5_enhanced and Raincloud plot

plot.zerosum5_raincloud <- select_data %>%
  filter(!is.na(ZEROSUM_5), !is.na(POLITICALPARTY)) %>%
  ggplot(aes(x = POLITICALPARTY, y = ZEROSUM_5, fill = POLITICALPARTY, color = POLITICALPARTY)) +
  stat_halfeye(
    adjust = 0.5,
    width = 0.6,
    .width = 0,
    justification = -0.2,
    point_colour = NA,
    alpha = 0.7
  ) +
  geom_pointrange(
      data = group_stats_5_avg,
      aes(
        x = POLITICALPARTY,
        y = mean_ZEROSUM_5,
        ymin = ci_lower,
        ymax = ci_upper
      ),
      color = "black",
      size = 0.9,
      shape = 21,
      fill = "white",
      stroke = 1.2,
      inherit.aes = FALSE
    ) +
  geom_point(
    aes(shape = POLITICALPARTY),
    size = 1.5,
    alpha = 0.5,
    position = position_jitter(
      seed = 1, width = 0.1
    )
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3,
    color = "white",
    stroke = 1.5,
    shape = 21
  ) +
  scale_fill_manual(values = party_colors) +
  scale_color_manual(values = party_colors) +
  theme_bw() +
  labs(
    title = "Less discrimination against minorities means more \n discrimination against whites.",
    x = "Political Affiliation",
    y = "Level of Agreement",
    caption = "White dots = means; colored dots = individual responses"
  ) +
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    plot.caption = element_text(size = 9, color = "gray40", hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    legend.position = "none",
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3)  # Fixed
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7) 
  )
plot.zerosum5_raincloud

```



```{r}
#| label: ZEROSUM_5 ANOVA

# Two-way ANOVA
model.ZEROSUM_5.PE <- lm(ZEROSUM_5 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)
Anova(model.ZEROSUM_5.PE)

## There is a statistically significant difference in ZEROSUM_5 scores between political parties.

```

```{r}
#| label: APA format ZEROSUM_5 ANOVA

# Two-way ANOVA
model.ZEROSUM_5.PE <- lm(ZEROSUM_5 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)

anova_df <- as.data.frame(Anova(model.ZEROSUM_5.PE, type = 2))%>%
  tibble::rownames_to_column(var = "Predictor")
names(anova_df)[names(anova_df) == "Pr(>F)"] <- "p"

# APA format
nice_table(anova_df, 
           title = c("Table 5", "ZeroSum_5 ANOVA Table"), 
           highlight = 0.05, 
           stars = TRUE)

## There is a statistically significant difference in ZEROSUM_5 scores between political parties.

```

```{r}
#| label: model.ZEROSUM_5.PE Tukey HSD test
#| include: false
#| eval: false
# Perform Tukey HSD test
tukey_result <- TukeyHSD(model.ZEROSUM_5.PE)

# Print the Tukey HSD results
print(tukey_result)

```

```{r}
#| label: ZEROSUM_5 ANOVA plot

plot(model.ZEROSUM_5.PE)  # diagnostic plots
shapiro.test(residuals(model.ZEROSUM_5.PE))  # normality test

```

```{r}
#| label: model.ZEROSUM_5.PE on residuals

# Check how bad the non-normality is
hist(residuals(model.ZEROSUM_5.PE))
# check: for normal distribution
# results: not normal, skew to the right
# interpretation: this is okay

qqnorm(residuals(model.ZEROSUM_5.PE))
# check: points should fall near the diagonal line 
# results: potential outliers at the extremes

qqline(residuals(model.ZEROSUM_5.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

```

```{r}
#| label: model.ZEROSUM_5.PE boxplot to check for outliers

boxplot(residuals(model.ZEROSUM_5.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

## overall: the visual inspection across the previous plots suggests we are okay to move forward with the ANOVA test, but must report the shapiro.test results here. 
```

### Zero Sum 6

```{r}
#| label: groupstats for ZERO_SUM6 

library(dplyr)
group_stats_6 <- select_data %>%
  group_by(POLITICALPARTY, RACIALIDENTITY.4) %>%
  summarise(
    mean_ZEROSUM_6 = mean(ZEROSUM_6, na.rm = TRUE),
    se = sd(ZEROSUM_6, na.rm = TRUE) / sqrt(n()),
    .groups = 'drop'
  )

# Check the structure
str(group_stats_6)
head(group_stats_6)
```

```{r}
#| label: zerosum6 pointrange

# Define the zesty_four palette
zesty_four <- c("#F5793A", "#A95AA1", "#85C0F9", "#0F2080") 

# Create the plot with pointrange
plot.zerosum6 <- ggplot(group_stats_6, aes(POLITICALPARTY, mean_ZEROSUM_6)) +
  geom_pointrange(aes(color = RACIALIDENTITY.4, 
                      ymin = mean_ZEROSUM_6 - se,  # Use standard error instead of sd
                      ymax = mean_ZEROSUM_6 + se), # Use standard error instead of sd
                  position = position_dodge(0.3)) +
  theme_bw() +
  labs(title = "More opportunity for transwomen means less opportunity \n for people who are assigned female at birth.",
       x = "Political Affiliation",
       y = "Agreement with Zero Sum Beliefs",
       color = "RACIALIDENTITY.4") + 
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8)
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7)  # Set the y-axis limits to match the scale
  )+
  scale_color_manual(values = zesty_four)

# Print and save the plots
print(plot.zerosum6)
ggsave("plot6:trans-cis.png", 
       plot = plot.zerosum6, 
       width = 10, 
       height = 8, 
       dpi = 300)
```


```{r}
#| label: ZEROSUM_6 average score with 95% CI

# ZEROSUM_6 average score with 95% CI
group_stats_6_avg <- select_data %>%
  group_by(POLITICALPARTY) %>%
  summarise(
    mean_ZEROSUM_6 = mean(ZEROSUM_6, na.rm = TRUE),
    se = sd(ZEROSUM_6, na.rm = TRUE) / sqrt(n()),
    n = n(),
    .groups = 'drop'
  ) %>%
  mutate(
    # Calculate 95% confidence interval using t-distribution
    t_value = qt(0.975, df = n - 1),  # 95% CI, two-tailed
    ci_lower = mean_ZEROSUM_6 - t_value * se,
    ci_upper = mean_ZEROSUM_6 + t_value * se
  )
```



```{r}
#| label: plot.zerosum6_enhanced and Raincloud plot

plot.zerosum6_raincloud <- select_data %>%
  filter(!is.na(ZEROSUM_6), !is.na(POLITICALPARTY)) %>%
  ggplot(aes(x = POLITICALPARTY, y = ZEROSUM_6, fill = POLITICALPARTY, color = POLITICALPARTY)) +
  stat_halfeye(
    adjust = 0.5,
    width = 0.6,
    .width = 0,
    justification = -0.2,
    point_colour = NA,
    alpha = 0.7
  ) +
  geom_pointrange(
      data = group_stats_6_avg,
      aes(
        x = POLITICALPARTY,
        y = mean_ZEROSUM_6,
        ymin = ci_lower,
        ymax = ci_upper
      ),
      color = "black",
      size = 0.9,
      shape = 21,
      fill = "white",
      stroke = 1.2,
      inherit.aes = FALSE
    ) +
  geom_point(
    aes(shape = POLITICALPARTY),
    size = 1.5,
    alpha = 0.5,
    position = position_jitter(
      seed = 1, width = 0.1
    )
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3,
    color = "white",
    stroke = 1.5,
    shape = 21
  ) +
  scale_fill_manual(values = party_colors) +
  scale_color_manual(values = party_colors) +
  theme_bw() +
  labs(
    title = "More opportunity for transwomen means less opportunity \n for people who are assigned female at birth.",
    x = "Political Affiliation",
    y = "Level of Agreement",
    caption = "White dots = means; colored dots = individual responses"
  ) +
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    plot.caption = element_text(size = 9, color = "gray40", hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    legend.position = "none",
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3)  # Fixed
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7) 
  )
plot.zerosum6_raincloud

```


```{r}
#| label: ZEROSUM_6 ANOVA

# Two-way ANOVA
model.ZEROSUM_6.PE <- lm(ZEROSUM_6 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)
Anova(model.ZEROSUM_6.PE)

## There is a statistically significant difference in ZEROSUM_6 scores between political parties.

```

```{r}
#| label: APA format ZEROSUM_6 ANOVA

# Two-way ANOVA
model.ZEROSUM_6.PE <- lm(ZEROSUM_6 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)

anova_df <- as.data.frame(Anova(model.ZEROSUM_6.PE, type = 2))%>%
  tibble::rownames_to_column(var = "Predictor")
names(anova_df)[names(anova_df) == "Pr(>F)"] <- "p"

# APA format
nice_table(anova_df, 
           title = c("Table 6", "ZeroSum_6 ANOVA Table"), 
           highlight = 0.05, 
           stars = TRUE)

## There is a statistically significant difference in ZEROSUM_6 scores between political parties.

```

```{r}
#| label: model.ZEROSUM_6.PE Tukey HSD test
#| include: false
#| eval: false
# Perform Tukey HSD test
tukey_result <- TukeyHSD(model.ZEROSUM_6.PE)

# Print the Tukey HSD results
print(tukey_result)
```

```{r}
#| label: ZEROSUM_6 ANOVA plot

plot(model.ZEROSUM_6.PE)  # diagnostic plots
shapiro.test(residuals(model.ZEROSUM_6.PE))  # normality test

```

```{r}
#| label: model.ZEROSUM_6.PE on residuals

# Check how bad the non-normality is
hist(residuals(model.ZEROSUM_6.PE))
# check: for normal distribution
# results: not normal
# interpretation: this is okay

qqnorm(residuals(model.ZEROSUM_6.PE))
# check: points should fall near the diagonal line 
# results: potential outliers at the extremes, some of the point in the middle are far away the line

qqline(residuals(model.ZEROSUM_6.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation


```

```{r}
#| label: model.ZEROSUM_6.PE boxplot to check for outliers

boxplot(residuals(model.ZEROSUM_6.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

## overall: the visual inspection across the previous plots suggests we are okay to move forward with the ANOVA test, but must report the shapiro.test results here. 

```

### Zero Sum 7

```{r}
#| label: groupstats for ZERO_SUM7 

library(dplyr)
group_stats_7 <- select_data %>%
  group_by(POLITICALPARTY, RACIALIDENTITY.4) %>%
  summarise(
    mean_ZEROSUM_7 = mean(ZEROSUM_7, na.rm = TRUE),
    se = sd(ZEROSUM_7, na.rm = TRUE) / sqrt(n()),
    .groups = 'drop'
  )

# Check the structure
str(group_stats_7)
head(group_stats_7)
```

```{r}
#| label: zerosum7 pointrange

# Define the zesty_four palette
zesty_four <- c("#F5793A", "#A95AA1", "#85C0F9", "#0F2080") 

# Create the plot with pointrange
plot.zerosum7 <- ggplot(group_stats_7, aes(POLITICALPARTY, mean_ZEROSUM_7)) +
  geom_pointrange(aes(color = RACIALIDENTITY.4, 
                      ymin = mean_ZEROSUM_7 - se,  # Use standard error instead of sd
                      ymax = mean_ZEROSUM_7 + se), # Use standard error instead of sd
                  position = position_dodge(0.3)) +
  theme_bw() +
  labs(title = "More health care access for undocumented immigrants \n means less access for U.S. citizens.",
       x = "Political Affiliation",
       y = "Agreement with Zero Sum Beliefs",
       color = "RACIALIDENTITY.4") + 
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8)
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7)  # Set the y-axis limits to match the scale
  )+
  scale_color_manual(values = zesty_four)

# Print and save the plots
print(plot.zerosum7)
ggsave("plot7:undoc-citizens.png", 
       plot = plot.zerosum7, 
       width = 10, 
       height = 8, 
       dpi = 300)
```


```{r}
#| label: ZEROSUM_7 average score with 95% CI

# ZEROSUM_7 average score with 95% CI
group_stats_7_avg <- select_data %>%
  group_by(POLITICALPARTY) %>%
  summarise(
    mean_ZEROSUM_7 = mean(ZEROSUM_7, na.rm = TRUE),
    se = sd(ZEROSUM_7, na.rm = TRUE) / sqrt(n()),
    n = n(),
    .groups = 'drop'
  ) %>%
  mutate(
    # Calculate 95% confidence interval using t-distribution
    t_value = qt(0.975, df = n - 1),  # 95% CI, two-tailed
    ci_lower = mean_ZEROSUM_7 - t_value * se,
    ci_upper = mean_ZEROSUM_7 + t_value * se
  )
```



```{r}
#| label: plot.zerosum7_enhanced and Raincloud plot

plot.zerosum7_raincloud <- select_data %>%
  filter(!is.na(ZEROSUM_7), !is.na(POLITICALPARTY)) %>%
  ggplot(aes(x = POLITICALPARTY, y = ZEROSUM_7, fill = POLITICALPARTY, color = POLITICALPARTY)) +
  stat_halfeye(
    adjust = 0.5,
    width = 0.6,
    .width = 0,
    justification = -0.2,
    point_colour = NA,
    alpha = 0.7
  ) +
  geom_pointrange(
      data = group_stats_7_avg,
      aes(
        x = POLITICALPARTY,
        y = mean_ZEROSUM_7,
        ymin = ci_lower,
        ymax = ci_upper
      ),
      color = "black",
      size = 0.9,
      shape = 21,
      fill = "white",
      stroke = 1.2,
      inherit.aes = FALSE
    ) +
  geom_point(
    aes(shape = POLITICALPARTY),
    size = 1.5,
    alpha = 0.5,
    position = position_jitter(
      seed = 1, width = 0.1
    )
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3,
    color = "white",
    stroke = 1.5,
    shape = 21
  ) +
  scale_fill_manual(values = party_colors) +
  scale_color_manual(values = party_colors) +
  theme_bw() +
  labs(
    title = "More health care access for undocumented immigrants \n means less access for U.S. citizens.",
    x = "Political Affiliation",
    y = "Level of Agreement",
    caption = "White dots = means; colored dots = individual responses"
  ) +
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    plot.caption = element_text(size = 9, color = "gray40", hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    legend.position = "none",
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3)  # Fixed
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7) 
  )
plot.zerosum7_raincloud

```


```{r}
#| label: ZEROSUM_7 ANOVA

# Two-way ANOVA
model.ZEROSUM_7.PE <- lm(ZEROSUM_7 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)
Anova(model.ZEROSUM_7.PE)

## There is a statistically significant difference in ZEROSUM_7 scores between political parties.

```

```{r}
#| label: APA format ZEROSUM_7 ANOVA

# Two-way ANOVA
model.ZEROSUM_7.PE <- lm(ZEROSUM_7 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)

anova_df <- as.data.frame(Anova(model.ZEROSUM_7.PE, type = 2))%>%
  tibble::rownames_to_column(var = "Predictor")
names(anova_df)[names(anova_df) == "Pr(>F)"] <- "p"

# APA format
nice_table(anova_df, 
           title = c("Table 7", "ZeroSum_7 ANOVA Table"), 
           highlight = 0.05, 
           stars = TRUE)

## There is a statistically significant difference in ZEROSUM_7 scores between political parties.

```

```{r}
#| label: model.ZEROSUM_7.PE Tukey HSD test
#| include: false
#| eval: false
# Perform Tukey HSD test
tukey_result <- TukeyHSD(model.ZEROSUM_7.PE)

# Print the Tukey HSD results
print(tukey_result)

```

```{r}
#| label: ZEROSUM_7 ANOVA plot

plot(model.ZEROSUM_7.PE)  # diagnostic plots
shapiro.test(residuals(model.ZEROSUM_7.PE))  # normality test

```

```{r}
#| label: model.ZEROSUM_7.PE on residuals

# Check how bad the non-normality is
hist(residuals(model.ZEROSUM_7.PE))
# check: for normal distribution
# results: not normal, but it's close (and likely would be with larger sample size)
# interpretation: this is okay

qqnorm(residuals(model.ZEROSUM_7.PE))
# check: points should fall near the diagonal line 
# results: potential outliers at the extremes

qqline(residuals(model.ZEROSUM_7.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

```

```{r}
#| label: model.ZEROSUM_7.PE boxplot to check for outliers

boxplot(residuals(model.ZEROSUM_7.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

## overall: the visual inspection across the previous plots suggests we are okay to move forward with the ANOVA test, but must report the shapiro.test results here. 

```

### Zero Sum 8

```{r}
#| label: groupstats for ZERO_SUM8 

library(dplyr)
group_stats_8 <- select_data %>%
  group_by(POLITICALPARTY, RACIALIDENTITY.4) %>%
  summarise(
    mean_ZEROSUM_8 = mean(ZEROSUM_8, na.rm = TRUE),
    se = sd(ZEROSUM_8, na.rm = TRUE) / sqrt(n()),
    .groups = 'drop'
  )

# Check the structure
str(group_stats_8)
head(group_stats_8)
```

```{r}
#| label: zerosum8 pointrange

# Define the zesty_four palette
zesty_four <- c("#F5793A", "#A95AA1", "#85C0F9", "#0F2080") 

# Create the plot with pointrange
plot.zerosum8 <- ggplot(group_stats_8, aes(POLITICALPARTY, mean_ZEROSUM_8)) +
  geom_pointrange(aes(color = RACIALIDENTITY.4, 
                      ymin = mean_ZEROSUM_8 - se,  # Use standard error instead of sd
                      ymax = mean_ZEROSUM_8 + se), # Use standard error instead of sd
                  position = position_dodge(0.3)) +
  theme_bw() +
  labs(title = "If there is equal pay for women, men will get lower wages.",
       x = "Political Affiliation",
       y = "Agreement with Zero Sum Beliefs",
       color = "RACIALIDENTITY.4") + 
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8)
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7)  # Set the y-axis limits to match the scale
  )+
  scale_color_manual(values = zesty_four)

# Print and save the plots
print(plot.zerosum8)
ggsave("plot8:paywomen-men.png", 
       plot = plot.zerosum8, 
       width = 10, 
       height = 8, 
       dpi = 300)
```

```{r}
#| label: ZEROSUM_8 average score with 95% CI

# ZEROSUM_8 average score with 95% CI
group_stats_8_avg <- select_data %>%
  group_by(POLITICALPARTY) %>%
  summarise(
    mean_ZEROSUM_8 = mean(ZEROSUM_8, na.rm = TRUE),
    se = sd(ZEROSUM_8, na.rm = TRUE) / sqrt(n()),
    n = n(),
    .groups = 'drop'
  ) %>%
  mutate(
    # Calculate 95% confidence interval using t-distribution
    t_value = qt(0.975, df = n - 1),  # 95% CI, two-tailed
    ci_lower = mean_ZEROSUM_8 - t_value * se,
    ci_upper = mean_ZEROSUM_8 + t_value * se
  )
```



```{r}
#| label: plot.zerosum8_enhanced and Raincloud plot

plot.zerosum8_raincloud <- select_data %>%
  filter(!is.na(ZEROSUM_8), !is.na(POLITICALPARTY)) %>%
  ggplot(aes(x = POLITICALPARTY, y = ZEROSUM_8, fill = POLITICALPARTY, color = POLITICALPARTY)) +
  stat_halfeye(
    adjust = 0.5,
    width = 0.6,
    .width = 0,
    justification = -0.2,
    point_colour = NA,
    alpha = 0.7
  ) +
  geom_pointrange(
      data = group_stats_8_avg,
      aes(
        x = POLITICALPARTY,
        y = mean_ZEROSUM_8,
        ymin = ci_lower,
        ymax = ci_upper
      ),
      color = "black",
      size = 0.9,
      shape = 21,
      fill = "white",
      stroke = 1.2,
      inherit.aes = FALSE
    ) +
  geom_point(
    aes(shape = POLITICALPARTY),
    size = 1.5,
    alpha = 0.5,
    position = position_jitter(
      seed = 1, width = 0.1
    )
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3,
    color = "white",
    stroke = 1.5,
    shape = 21
  ) +
  scale_fill_manual(values = party_colors) +
  scale_color_manual(values = party_colors) +
  theme_bw() +
  labs(
    title = "If there is equal pay for women, men will get lower wages.",
    x = "Political Affiliation",
    y = "Level of Agreement",
    caption = "White dots = means; colored dots = individual responses"
  ) +
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    plot.caption = element_text(size = 9, color = "gray40", hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    legend.position = "none",
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3)  # Fixed
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7) 
  )
plot.zerosum8_raincloud

```



```{r}
#| label: ZEROSUM_8 ANOVA

# Two-way ANOVA
model.ZEROSUM_8.PE <- lm(ZEROSUM_8 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)
Anova(model.ZEROSUM_8.PE)

## There is a statistically significant difference in ZEROSUM_8 scores between political parties. 
## There is a statistically significant difference in ZEROSUM_8 scores between ethnicities.

```

```{r}
#| label: APA format ZEROSUM_8 ANOVA

# Two-way ANOVA
model.ZEROSUM_8.PE <- lm(ZEROSUM_8 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)

anova_df <- as.data.frame(Anova(model.ZEROSUM_8.PE, type = 2))%>%
  tibble::rownames_to_column(var = "Predictor")
names(anova_df)[names(anova_df) == "Pr(>F)"] <- "p"

# APA format
nice_table(anova_df, 
           title = c("Table 8", "ZeroSum_8 ANOVA Table"), 
           highlight = 0.05, 
           stars = TRUE)

## There is a statistically significant difference in ZEROSUM_8 scores between political parties.
## There is a statistically significant difference in ZEROSUM_8 scores between ethnicities.
```

```{r}
#| label: model.ZEROSUM_8.PE Tukey HSD test
#| include: false
#| eval: false
# Perform Tukey HSD test
tukey_result <- TukeyHSD(model.ZEROSUM_8.PE)

# Print the Tukey HSD results
print(tukey_result)

```

```{r}
#| label: ZEROSUM_8 ANOVA plot

plot(model.ZEROSUM_8.PE)  # diagnostic plots
shapiro.test(residuals(model.ZEROSUM_8.PE))  # normality test

```

```{r}
#| label: model.ZEROSUM_8.PE on residuals

# Check how bad the non-normality is
hist(residuals(model.ZEROSUM_8.PE))
# check: for normal distribution
# results: approximately normal
# interpretation: this is okay

qqnorm(residuals(model.ZEROSUM_8.PE))
# check: points should fall near the diagonal line 
# results: potential outliers at the extremes

qqline(residuals(model.ZEROSUM_8.PE))
# check: for outliers
# results: no outliers
# interpretation: Since the distribution appears approximately normal, and the Shapiro-Wilk test gives a p-value of 0.1332 (which is greater than 0.05), we do not have evidence to reject normality. Therefore, it is appropriate to proceed with tests for significant differences.

```

```{r}
#| label: model.ZEROSUM_8.PE boxplot to check for outliers

boxplot(residuals(model.ZEROSUM_8.PE))
# check: for outliers
# results: no outliers
# interpretation: Since the distribution appears approximately normal, and the Shapiro-Wilk test gives a p-value of 0.1332 (which is greater than 0.05), we do not have evidence to reject normality. Therefore, it is appropriate to proceed with tests for significant differences.

## overall: the visual inspection across the previous plots suggests we are okay to move forward with the ANOVA test, but must report the shapiro.test results here. 

```

### Zero Sum 9

```{r}
#| label: groupstats for ZERO_SUM9 

library(dplyr)
group_stats_9 <- select_data %>%
  group_by(POLITICALPARTY, RACIALIDENTITY.4) %>%
  summarise(
    mean_ZEROSUM_9 = mean(ZEROSUM_9, na.rm = TRUE),
    se = sd(ZEROSUM_9, na.rm = TRUE) / sqrt(n()),
    .groups = 'drop'
  )

# Check the structure
str(group_stats_9)
head(group_stats_9)
```

```{r}
#| label: zerosum9 pointrange

# Define the zesty_four palette
zesty_four <- c("#F5793A", "#A95AA1", "#85C0F9", "#0F2080") 

# Create the plot with pointrange
plot.zerosum9 <- ggplot(group_stats_9, aes(POLITICALPARTY, mean_ZEROSUM_9)) +
  geom_pointrange(aes(color = RACIALIDENTITY.4, 
                      ymin = mean_ZEROSUM_9 - se,  # Use standard error instead of sd
                      ymax = mean_ZEROSUM_9 + se), # Use standard error instead of sd
                  position = position_dodge(0.3)) +
  theme_bw() +
  labs(title = "LGBTQ+ rights mean less freedom for religious groups.",
       x = "Political Affiliation",
       y = "Agreement with Zero Sum Beliefs",
       color = "RACIALIDENTITY.4") + 
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8)
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7)  # Set the y-axis limits to match the scale
  )+
  scale_color_manual(values = zesty_four)

# Print and save the plots
print(plot.zerosum9)
ggsave("plot9:LQBTQ-religious.png", 
       plot = plot.zerosum9, 
       width = 10, 
       height = 8, 
       dpi = 300)
```


```{r}
#| label: ZEROSUM_9 average score with 95% CI

# ZEROSUM_9 average score with 95% CI
group_stats_9_avg <- select_data %>%
  group_by(POLITICALPARTY) %>%
  summarise(
    mean_ZEROSUM_9 = mean(ZEROSUM_9, na.rm = TRUE),
    se = sd(ZEROSUM_9, na.rm = TRUE) / sqrt(n()),
    n = n(),
    .groups = 'drop'
  ) %>%
  mutate(
    # Calculate 95% confidence interval using t-distribution
    t_value = qt(0.975, df = n - 1),  # 95% CI, two-tailed
    ci_lower = mean_ZEROSUM_9 - t_value * se,
    ci_upper = mean_ZEROSUM_9 + t_value * se
  )
```



```{r}
#| label: plot.zerosum9_enhanced and Raincloud plot

plot.zerosum9_raincloud <- select_data %>%
  filter(!is.na(ZEROSUM_9), !is.na(POLITICALPARTY)) %>%
  ggplot(aes(x = POLITICALPARTY, y = ZEROSUM_9, fill = POLITICALPARTY, color = POLITICALPARTY)) +
  stat_halfeye(
    adjust = 0.5,
    width = 0.6,
    .width = 0,
    justification = -0.2,
    point_colour = NA,
    alpha = 0.7
  ) +
  geom_pointrange(
      data = group_stats_9_avg,
      aes(
        x = POLITICALPARTY,
        y = mean_ZEROSUM_9,
        ymin = ci_lower,
        ymax = ci_upper
      ),
      color = "black",
      size = 0.9,
      shape = 21,
      fill = "white",
      stroke = 1.2,
      inherit.aes = FALSE
    ) +
  geom_point(
    aes(shape = POLITICALPARTY),
    size = 1.5,
    alpha = 0.5,
    position = position_jitter(
      seed = 1, width = 0.1
    )
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3,
    color = "white",
    stroke = 1.5,
    shape = 21
  ) +
  scale_fill_manual(values = party_colors) +
  scale_color_manual(values = party_colors) +
  theme_bw() +
  labs(
    title = "LGBTQ+ rights mean less freedom for religious groups.",
    x = "Political Affiliation",
    y = "Level of Agreement",
    caption = "White dots = means; colored dots = individual responses"
  ) +
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    plot.caption = element_text(size = 9, color = "gray40", hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    legend.position = "none",
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3)  # Fixed
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7) 
  )
plot.zerosum9_raincloud

```



```{r}
#| label: ZEROSUM_9 ANOVA

# Two-way ANOVA
model.ZEROSUM_9.PE <- lm(ZEROSUM_9 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)
Anova(model.ZEROSUM_9.PE)

## There is a statistically significant difference in ZEROSUM_9 scores between political parties. 
## There is a statistically significant difference in ZEROSUM_9 scores between ethnicities.

```

```{r}
#| label: APA format ZEROSUM_9 ANOVA

# Two-way ANOVA
model.ZEROSUM_9.PE <- lm(ZEROSUM_9 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)

anova_df <- as.data.frame(Anova(model.ZEROSUM_9.PE, type = 2))%>%
  tibble::rownames_to_column(var = "Predictor")
names(anova_df)[names(anova_df) == "Pr(>F)"] <- "p"

# APA format
nice_table(anova_df, 
           title = c("Table 9", "ZeroSum_9 ANOVA Table"), 
           highlight = 0.05, 
           stars = TRUE)

## There is a statistically significant difference in ZEROSUM_9 scores between political parties.
## There is a statistically significant difference in ZEROSUM_9 scores between ethnicities.
```

```{r}
#| label: model.ZEROSUM_9.PE Tukey HSD test
#| include: false
#| eval: false
# Perform Tukey HSD test
tukey_result <- TukeyHSD(model.ZEROSUM_9.PE)

# Print the Tukey HSD results
print(tukey_result)

```

```{r}
#| label: ZEROSUM_9 ANOVA plot

plot(model.ZEROSUM_9.PE)  # diagnostic plots
shapiro.test(residuals(model.ZEROSUM_9.PE))  # normality test

```

```{r}
#| label: model.ZEROSUM_9.PE on residuals

# Check how bad the non-normality is
hist(residuals(model.ZEROSUM_9.PE))
# check: for normal distribution
# results: not normal, but it's close (and likely would be with larger sample size)
# interpretation: this is okay

qqnorm(residuals(model.ZEROSUM_9.PE))
# check: points should fall near the diagonal line 
# results: potential outliers at the extremes, some of the point in the middle are far away the line

qqline(residuals(model.ZEROSUM_9.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation


```

```{r}
#| label: model.ZEROSUM_9.PE boxplot to check for outliers

boxplot(residuals(model.ZEROSUM_9.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

## overall: the visual inspection across the previous plots suggests we are okay to move forward with the ANOVA test, but must report the shapiro.test results here. 

```

### Zero Sum 10

```{r}
#| label: groupstats for ZERO_SUM10 

library(dplyr)
group_stats_10 <- select_data %>%
  group_by(POLITICALPARTY, RACIALIDENTITY.4) %>%
  summarise(
    mean_ZEROSUM_10 = mean(ZEROSUM_10, na.rm = TRUE),
    se = sd(ZEROSUM_10, na.rm = TRUE) / sqrt(n()),
    .groups = 'drop'
  )

# Check the structure
str(group_stats_10)
head(group_stats_10)
```

```{r}
#| label: zerosum10 pointrange

# Define the zesty_four palette
zesty_four <- c("#F5793A", "#A95AA1", "#85C0F9", "#0F2080") 

# Create the plot with pointrange
plot.zerosum10 <- ggplot(group_stats_10, aes(POLITICALPARTY, mean_ZEROSUM_10)) +
  geom_pointrange(aes(color = RACIALIDENTITY.4, 
                      ymin = mean_ZEROSUM_10 - se,  # Use standard error instead of sd
                      ymax = mean_ZEROSUM_10 + se), # Use standard error instead of sd
                  position = position_dodge(0.3)) +
  theme_bw() +
  labs(title = "Accessible healthcare for people with disabilities means \n longer wait times for non-disabled patients.",
       x = "Political Affiliation",
       y = "Agreement with Zero Sum Beliefs",
       color = "RACIALIDENTITY.4") + 
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8)
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7)  # Set the y-axis limits to match the scale
  )+
  scale_color_manual(values = zesty_four)

# Print and save the plots
print(plot.zerosum10)
ggsave("plot10:disabilities-non.png", 
       plot = plot.zerosum10, 
       width = 10, 
       height = 8, 
       dpi = 300)

```

```{r}
#| label: ZEROSUM_10 average score with 95% CI

# ZEROSUM_10 average score with 95% CI
group_stats_10_avg <- select_data %>%
  group_by(POLITICALPARTY) %>%
  summarise(
    mean_ZEROSUM_10 = mean(ZEROSUM_10, na.rm = TRUE),
    se = sd(ZEROSUM_10, na.rm = TRUE) / sqrt(n()),
    n = n(),
    .groups = 'drop'
  ) %>%
  mutate(
    # Calculate 95% confidence interval using t-distribution
    t_value = qt(0.975, df = n - 1),  # 95% CI, two-tailed
    ci_lower = mean_ZEROSUM_10 - t_value * se,
    ci_upper = mean_ZEROSUM_10 + t_value * se
  )
```



```{r}
#| label: plot.zerosum10_enhanced and Raincloud plot

plot.zerosum10_raincloud <- select_data %>%
  filter(!is.na(ZEROSUM_10), !is.na(POLITICALPARTY)) %>%
  ggplot(aes(x = POLITICALPARTY, y = ZEROSUM_10, fill = POLITICALPARTY, color = POLITICALPARTY)) +
  stat_halfeye(
    adjust = 0.5,
    width = 0.6,
    .width = 0,
    justification = -0.2,
    point_colour = NA,
    alpha = 0.7
  ) +
  geom_pointrange(
      data = group_stats_10_avg,
      aes(
        x = POLITICALPARTY,
        y = mean_ZEROSUM_10,
        ymin = ci_lower,
        ymax = ci_upper
      ),
      color = "black",
      size = 0.9,
      shape = 21,
      fill = "white",
      stroke = 1.2,
      inherit.aes = FALSE
    ) +
  geom_point(
    aes(shape = POLITICALPARTY),
    size = 1.5,
    alpha = 0.5,
    position = position_jitter(
      seed = 1, width = 0.1
    )
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3,
    color = "white",
    stroke = 1.5,
    shape = 21
  ) +
  scale_fill_manual(values = party_colors) +
  scale_color_manual(values = party_colors) +
  theme_bw() +
  labs(
    title = "Accessible healthcare for people with disabilities means \n longer wait times for non-disabled patients.",
    x = "Political Affiliation",
    y = "Level of Agreement",
    caption = "White dots = means; colored dots = individual responses"
  ) +
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    plot.caption = element_text(size = 9, color = "gray40", hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    legend.position = "none",
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3)  # Fixed
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7) 
  )
plot.zerosum10_raincloud

```




```{r}
#| label: ZEROSUM_10 ANOVA

# Two-way ANOVA
model.ZEROSUM_10.PE <- lm(ZEROSUM_10 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)
Anova(model.ZEROSUM_10.PE)

## There is a statistically significant difference in ZEROSUM_10 scores between political parties.
```

```{r}
#| label: APA format ZEROSUM_10 ANOVA

# Two-way ANOVA
model.ZEROSUM_10.PE <- lm(ZEROSUM_10 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)

anova_df <- as.data.frame(Anova(model.ZEROSUM_10.PE, type = 2))%>%
  tibble::rownames_to_column(var = "Predictor")
names(anova_df)[names(anova_df) == "Pr(>F)"] <- "p"

# APA format
nice_table(anova_df, 
           title = c("Table 10", "ZeroSum_10 ANOVA Table"), 
           highlight = 0.05, 
           stars = TRUE)

## There is a statistically significant difference in ZEROSUM_10 scores between political parties.

```

```{r}
#| label: model.ZEROSUM_10.PE Tukey HSD test
#| include: false
#| eval: false
# Perform Tukey HSD test
tukey_result <- TukeyHSD(model.ZEROSUM_10.PE)

# Print the Tukey HSD results
print(tukey_result)
```

```{r}
#| label: ZEROSUM_10 ANOVA plot

plot(model.ZEROSUM_10.PE)  # diagnostic plots
shapiro.test(residuals(model.ZEROSUM_10.PE))  # normality test

```

```{r}
#| label: model.ZEROSUM_10.PE on residuals

# Check how bad the non-normality is
hist(residuals(model.ZEROSUM_10.PE))
# check: for normal distribution
# results: not normal, but it's close (and likely would be with larger sample size)
# interpretation: this is okay

qqnorm(residuals(model.ZEROSUM_10.PE))
# check: points should fall near the diagonal line 
# results: potential outliers at the extremes

qqline(residuals(model.ZEROSUM_10.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

```

```{r}
#| label: model.ZEROSUM_10.PE boxplot to check for outliers

boxplot(residuals(model.ZEROSUM_10.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

## overall: the visual inspection across the previous plots suggests we are okay to move forward with the ANOVA test, but must report the shapiro.test results here. 

```

### Zero Sum 11

```{r}
#| label: groupstats for ZERO_SUM11 

library(dplyr)
group_stats_11 <- select_data %>%
  group_by(POLITICALPARTY, RACIALIDENTITY.4) %>%
  summarise(
    mean_ZEROSUM_11 = mean(ZEROSUM_11, na.rm = TRUE),
    se = sd(ZEROSUM_11, na.rm = TRUE) / sqrt(n()),
    .groups = 'drop'
  )

# Check the structure
str(group_stats_11)
head(group_stats_11)
```

```{r}
#| label: zerosum11 pointrange

# Define the zesty_four palette
zesty_four <- c("#F5793A", "#A95AA1", "#85C0F9", "#0F2080") 

# Create the plot with pointrange
plot.zerosum11 <- ggplot(group_stats_11, aes(POLITICALPARTY, mean_ZEROSUM_11)) +
  geom_pointrange(aes(color = RACIALIDENTITY.4, 
                      ymin = mean_ZEROSUM_11 - se,  # Use standard error instead of sd
                      ymax = mean_ZEROSUM_11 + se), # Use standard error instead of sd
                  position = position_dodge(0.3)) +
  theme_bw() +
  labs(title = "Universal healthcare means worse healthcare for those who can afford private insurance.",
       x = "Political Affiliation",
       y = "Agreement with Zero Sum Beliefs",
       color = "RACIALIDENTITY.4") + 
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8)
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7)  # Set the y-axis limits to match the scale
  )+
  scale_color_manual(values = zesty_four)

# Print and save the plots
print(plot.zerosum11)
ggsave("plot11:healthcare.png", 
       plot = plot.zerosum11, 
       width = 10, 
       height = 8, 
       dpi = 300)


```


```{r}
#| label: ZEROSUM_11 average score with 95% CI

# ZEROSUM_11 average score with 95% CI
group_stats_11_avg <- select_data %>%
  group_by(POLITICALPARTY) %>%
  summarise(
    mean_ZEROSUM_11 = mean(ZEROSUM_11, na.rm = TRUE),
    se = sd(ZEROSUM_11, na.rm = TRUE) / sqrt(n()),
    n = n(),
    .groups = 'drop'
  ) %>%
  mutate(
    # Calculate 95% confidence interval using t-distribution
    t_value = qt(0.975, df = n - 1),  # 95% CI, two-tailed
    ci_lower = mean_ZEROSUM_11 - t_value * se,
    ci_upper = mean_ZEROSUM_11 + t_value * se
  )
```



```{r}
#| label: plot.zerosum11_enhanced and Raincloud plot

plot.zerosum11_raincloud <- select_data %>%
  filter(!is.na(ZEROSUM_11), !is.na(POLITICALPARTY)) %>%
  ggplot(aes(x = POLITICALPARTY, y = ZEROSUM_11, fill = POLITICALPARTY, color = POLITICALPARTY)) +
  stat_halfeye(
    adjust = 0.5,
    width = 0.6,
    .width = 0,
    justification = -0.2,
    point_colour = NA,
    alpha = 0.7
  ) +
  geom_pointrange(
      data = group_stats_11_avg,
      aes(
        x = POLITICALPARTY,
        y = mean_ZEROSUM_11,
        ymin = ci_lower,
        ymax = ci_upper
      ),
      color = "black",
      size = 0.9,
      shape = 21,
      fill = "white",
      stroke = 1.2,
      inherit.aes = FALSE
    ) +
  geom_point(
    aes(shape = POLITICALPARTY),
    size = 1.5,
    alpha = 0.5,
    position = position_jitter(
      seed = 1, width = 0.1
    )
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3,
    color = "white",
    stroke = 1.5,
    shape = 21
  ) +
  scale_fill_manual(values = party_colors) +
  scale_color_manual(values = party_colors) +
  theme_bw() +
  labs(
    title = "Universal healthcare means worse healthcare for those \n who can afford private insurance.",
    x = "Political Affiliation",
    y = "Level of Agreement",
    caption = "White dots = means; colored dots = individual responses"
  ) +
  theme(
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    plot.caption = element_text(size = 9, color = "gray40", hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    legend.position = "none",
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3)  # Fixed
  ) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("1: Strongly Disbelieve", 
               "2: Disbelieve", 
               "3: Somewhat Disbelieve", 
               "4: Neither", 
               "5: Somewhat Believe", 
               "6: Believe", 
               "7: Strongly Believe"),
    limits = c(1, 7) 
  )
plot.zerosum11_raincloud

```


```{r}
#| label: ZEROSUM_11 ANOVA

# Two-way ANOVA
model.ZEROSUM_11.PE <- lm(ZEROSUM_11 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)
Anova(model.ZEROSUM_11.PE)

## There is a statistically significant difference in ZEROSUM_11 scores between political parties.

```

```{r}
#| label: APA format ZEROSUM_11 ANOVA

# Two-way ANOVA
model.ZEROSUM_11.PE <- lm(ZEROSUM_11 ~ POLITICALPARTY * RACIALIDENTITY.4, data = select_data)

anova_df <- as.data.frame(Anova(model.ZEROSUM_11.PE, type = 2))%>%
  tibble::rownames_to_column(var = "Predictor")
names(anova_df)[names(anova_df) == "Pr(>F)"] <- "p"

# APA format
nice_table(anova_df, 
           title = c("Table 11", "ZeroSum_11 ANOVA Table"), 
           highlight = 0.05, 
           stars = TRUE)

## There is a statistically significant difference in ZEROSUM_11 scores between political parties.

```

```{r}
#| label: model.ZEROSUM_11.PE Tukey HSD test
#| include: false
#| eval: false
# Perform Tukey HSD test
tukey_result <- TukeyHSD(model.ZEROSUM_11.PE)

# Print the Tukey HSD results
print(tukey_result)

```

```{r}
#| label: ZEROSUM_11 ANOVA plot

plot(model.ZEROSUM_11.PE)  # diagnostic plots
shapiro.test(residuals(model.ZEROSUM_11.PE))  # normality test
```

```{r}
#| label: model.ZEROSUM_11.PE on residuals

# Check how bad the non-normality is
hist(residuals(model.ZEROSUM_11.PE))
# check: for normal distribution
# results: not normal, but it's close (and likely would be with larger sample size)
# interpretation: this is okay

qqnorm(residuals(model.ZEROSUM_11.PE))
# check: points should fall near the diagonal line 
# results: potential outliers at the extremes

qqline(residuals(model.ZEROSUM_11.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation


```

```{r}
#| label: model.ZEROSUM_11.PE boxplot to check for outliers

boxplot(residuals(model.ZEROSUM_11.PE))
# check: for outliers
# results: no outliers
# interpretation: even though we violate assumptions previously, there are not outliers to drop to improve interpretation

## overall: the visual inspection across the previous plots suggests we are okay to move forward with the ANOVA test, but must report the shapiro.test results here. 
```

# Discussion
