final

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

# List of packages
packages <- c("tidyverse", "fst", "modelsummary") # add any you need here

# Install packages if they aren't installed already
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)

# Load the packages
lapply(packages, library, character.only = TRUE)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

## [[1]]
##  [1] "lubridate" "forcats"   "stringr"   "dplyr"     "purrr"     "readr"    
##  [7] "tidyr"     "tibble"    "ggplot2"   "tidyverse" "stats"     "graphics" 
## [13] "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[2]]
##  [1] "fst"       "lubridate" "forcats"   "stringr"   "dplyr"     "purrr"    
##  [7] "readr"     "tidyr"     "tibble"    "ggplot2"   "tidyverse" "stats"    
## [13] "graphics"  "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[3]]
##  [1] "modelsummary" "fst"          "lubridate"    "forcats"      "stringr"     
##  [6] "dplyr"        "purrr"        "readr"        "tidyr"        "tibble"      
## [11] "ggplot2"      "tidyverse"    "stats"        "graphics"     "grDevices"   
## [16] "utils"        "datasets"     "methods"      "base"

library(tidyverse)
library(modelsummary)

# Load the data
germany_data <- read_fst("germany_data.fst")

# preparing the data for the table1
germany_data <- germany_data %>%
  mutate(
    Ideology = case_when(
      lrscale >= 1 & lrscale <= 2 ~ "Extreme Left",
      lrscale > 2 & lrscale <= 4 ~ "Left",
      lrscale > 4 & lrscale <= 6 ~ "Moderate",
      lrscale > 6 & lrscale <= 8 ~ "Right",
      lrscale > 8 & lrscale <= 10 ~ "Extreme Right",
      TRUE ~ NA_character_  # Keeps NA values as is, ensuring character type
    ),
    progressiveness = case_when(
      hmsfmlsh == 1 ~ "extremely ashamed",
      hmsfmlsh == 2 ~ "ashamed",
      hmsfmlsh == 3 ~ "indifferent",
      hmsfmlsh == 4 ~ "supportive",
      hmsfmlsh == 5 ~ "extremely supportive",
      TRUE ~ NA_character_  # Keeps other values as NA
    ),
    religious_belief = case_when(
      rlgblg == 1 ~ "religious",
      rlgblg == 2 ~ "atheist",
      TRUE ~ NA_character_  # Keeps other values as NA
    )
  ) %>%
  # Drops the rows where any of these is NA
  drop_na(Ideology, progressiveness, religious_belief) %>%
  # Converting the cleaned variables to factors
  mutate(
    Ideology = factor(Ideology, levels = c("Extreme Left", "Left", "Moderate", "Right", "Extreme Right")),
    progressiveness = factor(progressiveness, levels = c("extremely ashamed", "ashamed", "indifferent", "supportive", "extremely supportive")),
    religious_belief = factor(religious_belief, levels = c("religious", "atheist"))
  )

# Ensuring that Ideology, progressiveness, and religious_belief are factors with the correct levels
germany_data <- germany_data %>%
  mutate(
    Ideology = factor(Ideology, levels = c("Extreme Left", "Left", "Moderate", "Right", "Extreme Right")),
    progressiveness = factor(progressiveness, levels = c("extremely ashamed", "ashamed", "indifferent", "supportive", "extremely supportive")),
    religious_belief = factor(religious_belief, levels = c("religious", "atheist"))
  )

# Filtering out NAs
germany_data <- germany_data %>%
  filter(!is.na(Ideology), !is.na(progressiveness), !is.na(religious_belief))

# Create a new table1 with the title as descriptive statistics for main variables
table1 <- datasummary_skim(
  germany_data %>% 
    select(Ideology, progressiveness, religious_belief),
  type = "categorical",
  title = "Table 1. Descriptive statistics for main variables",
  output = "flextable"
)

# printing the table
table1

Table 1. Descriptive statistics for main variables
		N	%
Ideology	Extreme Left	528	11.0
	Left	1432	29.9
	Moderate	2276	47.5
	Right	483	10.1
	Extreme Right	73	1.5
progressiveness	extremely ashamed	95	2.0
	ashamed	251	5.2
	indifferent	578	12.1
	supportive	1322	27.6
	extremely supportive	2546	53.1
religious_belief	religious	2681	55.9
	atheist	2111	44.1

table 1 is done

library(infer)
library(tidyverse)

# Step 1: Specifying and calculating the test statistic
test_result <- germany_data %>%
  specify(response = progressiveness, explanatory = Ideology) %>%
  hypothesize(null = "independence") %>%
  calculate(stat = "Chisq")  # Use 'Chisq' for categorical vs. categorical data

# Step 2: Generating the null distribution and calculate p-value
null_distribution <- germany_data %>%
  specify(response = progressiveness, explanatory = Ideology) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "Chisq")

# Step 3: Calculating the p-value
p_value <- null_distribution %>%
  get_pvalue(obs_stat = test_result, direction = "greater")

## Warning: Please be cautious in reporting a p-value of 0. This result is an approximation
## based on the number of `reps` chosen in the `generate()` step.
## ℹ See `get_p_value()` (`?infer::get_p_value()`) for more information.

# Step 4: Visualize the null distribution and test statistic
null_distribution %>%
  visualize() +
  geom_vline(aes(xintercept = test_result$stat), color = "red", linetype = "dashed") +
  shade_pvalue(obs_stat = test_result, direction = "greater")

## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf

# Printing the p-value and the visualization so we know what it is 
print(p_value)

## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1       0

library(dplyr)
library(ggplot2)
#now for the graph of our choice, we chose to use a graph that'll show us the relationship between shame towards homosexuality and political scale/ideology.
germany_data <- germany_data %>%
  mutate(
    Progressiveness = case_when(
      hmsfmlsh == 1 ~ "Extremely Ashamed",
      hmsfmlsh == 2 ~ "Ashamed",
      hmsfmlsh == 3 ~ "Indifferent",
      hmsfmlsh == 4 ~ "Supportive",
      hmsfmlsh == 5 ~ "Extremely Supportive",
      TRUE ~ NA_character_  # Set anything that is not 1, 2, 3, 4, or 5 to NA
    ),
    Ideology = case_when(
      lrscale >= 1 & lrscale <= 2 ~ "Extreme Left",
      lrscale > 2 & lrscale <= 4 ~ "Left",
      lrscale > 4 & lrscale <= 6 ~ "Moderate",
      lrscale > 6 & lrscale <= 8 ~ "Right",
      lrscale > 8 & lrscale <= 10 ~ "Extreme Right",
      TRUE ~ NA_character_  # For anything outside the 1-10 range or missing
    )
  ) %>%
  filter(!is.na(Progressiveness), !is.na(Ideology))  # Filter out missing data

# Calculate the percentages
ideology_progressiveness_percentages <- germany_data %>%
  group_by(Ideology, Progressiveness) %>%
  summarise(count = n(), .groups = 'drop') %>%
  mutate(percentage = count / sum(count) * 100) 

# Plotting
ideology_progressiveness_plot <- ggplot(ideology_progressiveness_percentages, aes(x = Ideology, y = percentage, fill = Progressiveness)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Ideology", y = "Percentage", title = "Ideological Views on Homosexuality in Germany") +
  scale_fill_brewer(palette = "Spectral") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Display the data and plot
print(ideology_progressiveness_percentages)

## # A tibble: 25 × 4
##    Ideology      Progressiveness      count percentage
##    <chr>         <chr>                <int>      <dbl>
##  1 Extreme Left  Ashamed                 23      0.480
##  2 Extreme Left  Extremely Ashamed        7      0.146
##  3 Extreme Left  Extremely Supportive   336      7.01 
##  4 Extreme Left  Indifferent             37      0.772
##  5 Extreme Left  Supportive             125      2.61 
##  6 Extreme Right Ashamed                  6      0.125
##  7 Extreme Right Extremely Ashamed        8      0.167
##  8 Extreme Right Extremely Supportive    28      0.584
##  9 Extreme Right Indifferent             13      0.271
## 10 Extreme Right Supportive              18      0.376
## # ℹ 15 more rows

print(ideology_progressiveness_plot)

###this is the graph of our choice.

germany_data <- germany_data %>%
  mutate(
    Ideology = factor(Ideology, levels = c("Extreme Left", "Left", "Moderate", "Right", "Extreme Right")),
    progressiveness = factor(progressiveness, levels = c("extremely ashamed", "ashamed", "indifferent", "supportive", "extremely supportive")),
    religious_belief = factor(religious_belief, levels = c("religious", "atheist"))
  )

# First model: Single predictor (Ideology)
model1 <- lm(progressiveness ~ Ideology, data = germany_data)

## Warning in model.response(mf, "numeric"): using type = "numeric" with a factor
## response will be ignored

## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors

# Second model: Two predictors (Ideology and ReligiousBelief)
model2 <- lm(progressiveness ~ Ideology + religious_belief, data = germany_data)

## Warning in model.response(mf, "numeric"): using type = "numeric" with a factor
## response will be ignored

## Warning in model.response(mf, "numeric"): '-' not meaningful for factors

# Third model: Interaction model (Ideology and ReligiousBelief)
model3 <- lm(progressiveness ~ Ideology * religious_belief, data = germany_data)

## Warning in model.response(mf, "numeric"): using type = "numeric" with a factor
## response will be ignored

## Warning in model.response(mf, "numeric"): '-' not meaningful for factors

# Using the modelsummary package to create a table
table <- modelsummary(list(model1, model2, model3),
                      title = "Regression Model Summary",
                      fmt = 2,
                      gof_omit = "Observations")

## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful for
## factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful for
## factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful for
## factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors

## Warning in Ops.factor(res, 2): '^' not meaningful for factors

# Output the table
table

Regression Model Summary
	(1)	(2)	(3)
(Intercept)	4.44	4.32	4.38
IdeologyLeft	−0.05	−0.03	−0.05
IdeologyModerate	−0.24	−0.21	−0.28
IdeologyRight	−0.54	−0.49	−0.64
IdeologyExtreme Right	−0.73	−0.68	−0.79
religious_beliefatheist		0.22	0.11
IdeologyLeft × religious_beliefatheist			0.02
IdeologyModerate × religious_beliefatheist			0.14
IdeologyRight × religious_beliefatheist			0.38
IdeologyExtreme Right × religious_beliefatheist			0.26
Num.Obs.	4792	4792	4792

###this is our interaction plot chart

library(ggplot2)
library(dplyr)

# Fit the model
model3 <- lm(progressiveness ~ Ideology * religious_belief, data = germany_data)

## Warning in model.response(mf, "numeric"): using type = "numeric" with a factor
## response will be ignored

## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors

# Create a new data frame for predictions
prediction_data <- expand.grid(
  Ideology = levels(germany_data$Ideology),
  religious_belief = levels(germany_data$religious_belief),
  progressiveness = mean(as.numeric(levels(germany_data$progressiveness)))  # Use mean progressiveness or another representative value
)

## Warning in mean(as.numeric(levels(germany_data$progressiveness))): NAs
## introduced by coercion

# Predict values
prediction_data$predicted_progressiveness <- predict(model3, newdata = prediction_data)

# Plot the interaction
interaction_plot <- ggplot(prediction_data, aes(x = Ideology, y = predicted_progressiveness, color = religious_belief, group = religious_belief)) +
  geom_line() +
  geom_point() +
  labs(
    title = "Interaction of Ideology and Religious Belief on Progressiveness",
    x = "Ideology",
    y = "Predicted Progressiveness",
    color = "Religious Belief"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))  # Improve label readability

# Print the plot
print(interaction_plot)

##this is the interaction plot visualization with ggplot

final

lacey vila

2024-04-20

R Markdown

Including Plots

table 1 is done