Research Report SOC222
Quincy Hou
2024-04-12
Initial set-up
options(repos = c(CRAN = "https://cran.rstudio.com"))
# Install necessary packages if you haven't already
install.packages(c("haven", "fst"), repos='https://cran.rstudio.com')## Installing packages into 'C:/Users/quinc/AppData/Local/R/win-library/4.3'
## (as 'lib' is unspecified)## package 'haven' successfully unpacked and MD5 sums checked## Warning: cannot remove prior installation of package 'haven'## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying
## C:\Users\quinc\AppData\Local\R\win-library\4.3\00LOCK\haven\libs\x64\haven.dll
## to C:\Users\quinc\AppData\Local\R\win-library\4.3\haven\libs\x64\haven.dll:
## Permission denied## Warning: restored 'haven'## package 'fst' successfully unpacked and MD5 sums checked## Warning: cannot remove prior installation of package 'fst'## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying
## C:\Users\quinc\AppData\Local\R\win-library\4.3\00LOCK\fst\libs\x64\fst.dll to
## C:\Users\quinc\AppData\Local\R\win-library\4.3\fst\libs\x64\fst.dll: Permission
## denied## Warning: restored 'fst'##
## The downloaded binary packages are in
## C:\Users\quinc\AppData\Local\Temp\RtmpMDROqd\downloaded_packages# Load the packages
library(haven)## Warning: package 'haven' was built under R version 4.3.3library(fst)## Warning: package 'fst' was built under R version 4.3.3library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "MASS", "aod", "interactions", "kableExtra", "flextable", "scales") # 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 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.0
## ✔ readr 2.1.5## ── 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## Warning: package 'infer' was built under R version 4.3.3## Warning: package 'effects' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3## lattice theme set by effectsTheme()
## See ?effectsTheme for details.## Warning: package 'survey' was built under R version 4.3.3## Loading required package: grid
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Loading required package: survival
##
## Attaching package: 'survey'
##
## The following object is masked from 'package:graphics':
##
## dotchart
##
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select## Warning: package 'aod' was built under R version 4.3.3##
## Attaching package: 'aod'
##
## The following object is masked from 'package:survival':
##
## rats## Warning: package 'interactions' was built under R version 4.3.3##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rows## Warning: package 'flextable' was built under R version 4.3.3##
## Attaching package: 'flextable'
##
## The following objects are masked from 'package:kableExtra':
##
## as_image, footnote
##
## The following object is masked from 'package:purrr':
##
## compose
##
##
## Attaching package: 'scales'
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factor## [[1]]
## [1] "lubridate" "forcats" "stringr" "purrr" "readr" "tidyr"
## [7] "tibble" "tidyverse" "ggplot2" "dplyr" "fst" "haven"
## [13] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [19] "base"
##
## [[2]]
## [1] "infer" "lubridate" "forcats" "stringr" "purrr" "readr"
## [7] "tidyr" "tibble" "tidyverse" "ggplot2" "dplyr" "fst"
## [13] "haven" "stats" "graphics" "grDevices" "utils" "datasets"
## [19] "methods" "base"
##
## [[3]]
## [1] "infer" "lubridate" "forcats" "stringr" "purrr" "readr"
## [7] "tidyr" "tibble" "tidyverse" "ggplot2" "dplyr" "fst"
## [13] "haven" "stats" "graphics" "grDevices" "utils" "datasets"
## [19] "methods" "base"
##
## [[4]]
## [1] "modelsummary" "infer" "lubridate" "forcats" "stringr"
## [6] "purrr" "readr" "tidyr" "tibble" "tidyverse"
## [11] "ggplot2" "dplyr" "fst" "haven" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"
##
## [[5]]
## [1] "effects" "carData" "modelsummary" "infer" "lubridate"
## [6] "forcats" "stringr" "purrr" "readr" "tidyr"
## [11] "tibble" "tidyverse" "ggplot2" "dplyr" "fst"
## [16] "haven" "stats" "graphics" "grDevices" "utils"
## [21] "datasets" "methods" "base"
##
## [[6]]
## [1] "survey" "survival" "Matrix" "grid" "effects"
## [6] "carData" "modelsummary" "infer" "lubridate" "forcats"
## [11] "stringr" "purrr" "readr" "tidyr" "tibble"
## [16] "tidyverse" "ggplot2" "dplyr" "fst" "haven"
## [21] "stats" "graphics" "grDevices" "utils" "datasets"
## [26] "methods" "base"
##
## [[7]]
## [1] "MASS" "survey" "survival" "Matrix" "grid"
## [6] "effects" "carData" "modelsummary" "infer" "lubridate"
## [11] "forcats" "stringr" "purrr" "readr" "tidyr"
## [16] "tibble" "tidyverse" "ggplot2" "dplyr" "fst"
## [21] "haven" "stats" "graphics" "grDevices" "utils"
## [26] "datasets" "methods" "base"
##
## [[8]]
## [1] "aod" "MASS" "survey" "survival" "Matrix"
## [6] "grid" "effects" "carData" "modelsummary" "infer"
## [11] "lubridate" "forcats" "stringr" "purrr" "readr"
## [16] "tidyr" "tibble" "tidyverse" "ggplot2" "dplyr"
## [21] "fst" "haven" "stats" "graphics" "grDevices"
## [26] "utils" "datasets" "methods" "base"
##
## [[9]]
## [1] "interactions" "aod" "MASS" "survey" "survival"
## [6] "Matrix" "grid" "effects" "carData" "modelsummary"
## [11] "infer" "lubridate" "forcats" "stringr" "purrr"
## [16] "readr" "tidyr" "tibble" "tidyverse" "ggplot2"
## [21] "dplyr" "fst" "haven" "stats" "graphics"
## [26] "grDevices" "utils" "datasets" "methods" "base"
##
## [[10]]
## [1] "kableExtra" "interactions" "aod" "MASS" "survey"
## [6] "survival" "Matrix" "grid" "effects" "carData"
## [11] "modelsummary" "infer" "lubridate" "forcats" "stringr"
## [16] "purrr" "readr" "tidyr" "tibble" "tidyverse"
## [21] "ggplot2" "dplyr" "fst" "haven" "stats"
## [26] "graphics" "grDevices" "utils" "datasets" "methods"
## [31] "base"
##
## [[11]]
## [1] "flextable" "kableExtra" "interactions" "aod" "MASS"
## [6] "survey" "survival" "Matrix" "grid" "effects"
## [11] "carData" "modelsummary" "infer" "lubridate" "forcats"
## [16] "stringr" "purrr" "readr" "tidyr" "tibble"
## [21] "tidyverse" "ggplot2" "dplyr" "fst" "haven"
## [26] "stats" "graphics" "grDevices" "utils" "datasets"
## [31] "methods" "base"
##
## [[12]]
## [1] "scales" "flextable" "kableExtra" "interactions" "aod"
## [6] "MASS" "survey" "survival" "Matrix" "grid"
## [11] "effects" "carData" "modelsummary" "infer" "lubridate"
## [16] "forcats" "stringr" "purrr" "readr" "tidyr"
## [21] "tibble" "tidyverse" "ggplot2" "dplyr" "fst"
## [26] "haven" "stats" "graphics" "grDevices" "utils"
## [31] "datasets" "methods" "base"data <- read_spss("romania_data.sav")
write_fst(data, "romania_data.fst")
romania_data <- read.fst("romania_data.fst")romania_data <- read.fst("romania_data.fst")
romania_data_cleaned <- romania_data %>%
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
TRUE ~ NA_character_
),
lrscale = as.numeric(lrscale),
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
TRUE ~ "No BA"
)
) %>%
filter(
!is.na(gndr), !is.na(educ.ba), !is.na(trstplt), !is.na(lrscale) # Now filter out the NA values
)
num_obs <- nrow(romania_data_cleaned )
print(paste("Number of observations", num_obs))## [1] "Number of observations 47015"romania_data <- read.fst("romania_data.fst")
## Calculating Descriptive Statistics for Trust in Parliament (trstplt)
trstplt_stats <- romania_data %>%
filter(!is.na(trstplt)) %>% # Ensure you're only including non-NA values
summarise(
mean_trstplt = mean(trstplt), # Calculate mean
sd_trstplt = sd(trstplt) # Calculate standard deviation
)
## Output trstplt statistics
print(trstplt_stats)## mean_trstplt sd_trstplt
## 1 3.275472 2.423905## Calculating Descriptive Statistics for Left-Right Scale (lrscale)
# First, ensure the lrscale is properly cleaned and numeric
romania_data <- romania_data %>%
mutate(lrscale = as.numeric(lrscale)) %>% # Convert to numeric if not already
filter(!is.na(lrscale) & lrscale <= 10) # Exclude any NA or invalid values
lrscale_stats <- romania_data %>%
summarise(
mean_lrscale = mean(lrscale), # Calculate mean
sd_lrscale = sd(lrscale) # Calculate standard deviation
)
## Output lrscale statistics
print(lrscale_stats)## mean_lrscale sd_lrscale
## 1 5.1991 2.278228##LRSCALE##
romania_data <- romania_data %>%
# Filter for France
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
TRUE ~ NA_character_ # Set anything that is not 1 or 2 to NA
),
lrscale = case_when(
lrscale %in% 0:3 ~ "Left", # Left-wing (0 to 3)
lrscale %in% 7:10 ~ "Right", # Right-wing (7 to 10)
TRUE ~ NA_character_ # Moderate (4, 5, 6) and special codes (77, 88, 99) set to NA
)
) TASK 4 EXPLANATION
Calculate the left right party percentages by filtering the left right scale and gender.
lrscale_percentages <- romania_data %>% # Begin with the dataset 'france_data'
filter(!is.na(lrscale), !is.na(gndr)) %>% # Filter out rows where 'lrscale' or 'gender' is NA (missing data)
group_by(gndr, lrscale) %>% # Group the data by 'gender' and 'lrscale' categories
summarise(count = n(), .groups = 'drop') %>% # Summarise each group to get counts, and then drop groupings
mutate(percentage = count / sum(count) * 100) # Calculate percentage for each group by dividing count by total count and multiplying by 100
lrscale_percentages # The resulting dataframe## # A tibble: 4 × 4
## gndr lrscale count percentage
## <chr> <chr> <int> <dbl>
## 1 Female Left 5283 23.0
## 2 Female Right 6444 28.1
## 3 Male Left 4800 20.9
## 4 Male Right 6433 28.0TASK 4 ANSWER
Visualization
# Create a ggplot object for horizontal bar chart with the specified style
lrscale_plot_v2 <- ggplot(lrscale_percentages,
aes(x = percentage, # Use percentage directly
y = reorder(gndr, -percentage), # Order bars within each gender
fill = gndr)) + # Fill color based on Gender
# Create horizontal bar chart
geom_col() + # Draws the bars using the provided data
coord_flip() + # Flip coordinates to make bars horizontal
# Remove fill color legend
guides(fill = "none") + # Removes legend for the fill aesthetic
# Split the plot based on Political Orientation
facet_wrap(~ lrscale, nrow = 1) + # Separate plots for Left/Right
# Labels and titles for the plot
labs(x = "Percentage of Respondents", # X-axis label
y = NULL, # Remove Y-axis label
title = "Political Orientation by Gender", # Main title
subtitle = "Comparing the percentage distribution of left vs. right for Romania (2008)") + # Subtitle
# Adjust visual properties of the plot
theme(plot.title = element_text(size = 16, face = "bold"), # Format title
plot.subtitle = element_text(size = 12), # Format subtitle
axis.title.y = element_blank(), # Remove Y-axis title
legend.position = "bottom") # Position the legend at the bottom
# Display the ggplot object
lrscale_plot_v2
##TRUST IN PARLIMENT##
romania_data <- read.fst("romania_data.fst")romania_data <- romania_data %>%
# Filter for France
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
TRUE ~ NA_character_ # Set anything that is not 1 or 2 to NA
),
trstplt = ifelse(trstplt %in% c(77, 88, 99), NA, trstplt), # set values 77, 88, and 99 to NA.
) table(romania_data$trstplt) # if there are values that are not supposed to be there (e.g., 77, 88, 99 in this case), then we need to deal with it##
## 0 1 2 3 4 5 6 7 8 9 10
## 10387 5651 6563 7123 6395 9101 4565 3387 1562 388 332trust_percentages <- romania_data %>%
filter(!is.na(trstplt), !is.na(gndr)) %>% # Filter out missing data
count(gndr, trstplt) %>% # Count the number of responses by gender and trust in politics
group_by(gndr) %>% # Group by gender
mutate(percentage = n / sum(n) * 100) # Calculate percentage within each gendertrust_plot <- ggplot(trust_percentages, aes(x = trstplt, y = percentage, fill = gndr)) +
geom_col() +
coord_flip() + # Flip coordinates to make bars horizontal
facet_wrap(~ gndr, nrow = 1, scales = "free_y") + # Separate plots for each gender
labs(x = "Trust in Politics Level", y = "Percentage of Respondents", title = "Trust in Parliament by Gender in Romania (2008)") +
theme_minimal() +
scale_fill_manual(values = c("Male" = "blue", "Female" = "red")) # Set colors for genders
# Display the plot
trust_plot
##BA OR NO BA##
romania_data <- read.fst("romania_data.fst")romania_data <- romania_data %>%
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
TRUE ~ NA_character_
),
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
TRUE ~ "No BA"
)
) %>%
filter(!is.na(gndr) & !is.na(educ.ba))education_by_gender <- romania_data %>%
group_by(gndr, educ.ba) %>%
summarise(count = n(), .groups = 'drop') %>%
mutate(percentage = count / sum(count) * 100)education_plot <- ggplot(education_by_gender, aes(x = educ.ba, y = count, fill = educ.ba)) +
geom_bar(stat = "identity", position = "dodge") +
facet_wrap(~gndr) +
labs(title = "Educational Attainment by Gender in Romania (2008)",
x = "Educational Level",
y = "Count of Respondents") +
scale_fill_manual(values = c("BA or more" ="Green", "No BA" = "Yellow")) +
theme_minimal()
education_plot
##GNDR X TRSTPLT##
romania_data <- read.fst("romania_data.fst")Step 1: Calculate the test statistic of your sample
romania_data <- romania_data %>%
mutate(
# Clean 'trstplt' values
trstplt = ifelse(trstplt %in% c(77, 88, 99), NA, trstplt),
# Recode 'gndr'
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_, # Assuming '9' is a code for missing data
TRUE ~ as.character(gndr)
),
# Recode 'edulvla' to 'educ.ba'
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more", # Assuming '5' indicates bachelor's degree or higher for rounds < 5
TRUE ~ "No BA" # All other cases are considered 'No BA'
),
# Clean 'edulvla' values
edulvla = ifelse(edulvla %in% c(77, 88), NA_integer_, edulvla)
) %>%
mutate(
# Force 'educ.ba' to be treated as character to avoid factor conversion issues in ggplot
educ.ba = as.character(educ.ba)
) %>%
# Remove rows with missing data for the analysis
filter(!is.na(gndr), !is.na(educ.ba), !is.na(trstplt))Step 1: Calculate the observed test statistic
test_stat <- romania_data %>%
specify(explanatory = gndr, # change variable name for explanatory variable
response = trstplt) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "diff in means", order = c("Female", "Male")) ## Message: The independence null hypothesis does not inform calculation of the
## observed statistic (a difference in means) and will be ignored.print(test_stat$stat)## [1] -0.04490246Step 2: Simulate the null distribution
null_distribution <- romania_data %>%
specify(response = trstplt, explanatory = gndr) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "diff in means", order = c("Female", "Male"))Step 3: Calculate the p-value of your sample
p_val <- null_distribution %>%
get_pvalue(obs_stat = test_stat, direction = "two-sided")
p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0.024Step 4: Visualize
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided")
null_distribution## Response: trstplt (numeric)
## Explanatory: gndr (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 -0.00202
## 2 2 0.0324
## 3 3 0.0266
## 4 4 -0.0417
## 5 5 -0.0136
## 6 6 0.00445
## 7 7 0.00329
## 8 8 -0.00231
## 9 9 -0.0182
## 10 10 0.0433
## # ℹ 990 more rows##GNDR X LRSCALE##
Step 1: Calculate the test statistic of your sample
romania_data <- read.fst("romania_data.fst")romania_data <- romania_data %>%
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
TRUE ~ NA_character_ # Set anything that is not 1 or 2 to NA
),
LR = case_when(
lrscale %in% 0:3 ~ "Left",
lrscale %in% 7:10 ~ "Right",
lrscale %in% 4:6 ~ "Moderate",
lrscale %in% c(77, 88, 99) ~ NA_character_
)
) %>%
# Remove rows with missing data for the analysis
filter(!is.na(gndr), !is.na(lrscale))##MODELS for GNDR & TRSPLT##
romania_data <- read.fst("romania_data.fst")Renaming dataset to df.
df <- romania_dataNote: you have both df and france_data which contain the french data (same observations, under different names). So you can always revert back to france_data if needed
Survey weights
From ess website:
For data collection the ess uses strictly probability-based samples. Every element in the ess target population should therefore have a greater than zero probability of being included into the sample. When analysing ess data estimates, the likelihood of each respondent to be part of the sample should also be taken into account - which means that the most accurate estimates will be obtained only after weighting the data.
df$weight <- df$dweight * df$pweight
survey_design <- svydesign(ids = ~1, data = df, weights = ~weight)Multiple linear regression
We will work with our populist attitudes scale based on Norris and Inglehart’s conceptualization and measurement again. So let’s run our code:
Setting up our populist scale
# Setting values greater than 10 to NA
df$trstplt <- ifelse(df$trstplt > 10, NA, df$trstplt)
df$trstprl <- ifelse(df$trstprl > 10, NA, df$trstprl)
df$trstprt <- ifelse(df$trstprt > 10, NA, df$trstprt)df <- df %>%
# Filter for France
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
TRUE ~ NA_character_ # Set anything that is not 1 or 2 to NA
),
lrscale = case_when(
lrscale %in% 0:3 ~ "Left", # Left-wing (0 to 3)
lrscale %in% 7:10 ~ "Right", # Right-wing (7 to 10)
TRUE ~ NA_character_ # Moderate (4, 5, 6) and special codes (77, 88, 99) set to NA
),
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
TRUE ~ "No BA"
)
) model1 <- lm(trstplt ~ gndr, data = df, weights = weight)
model2 <- lm(trstplt ~ gndr + educ.ba + lrscale, data = df, weights = weight)
model3 <- lm(trstplt ~ gndr*educ.ba, data = df, weights = weight)
modelsummary(
list(model1, model2, model3),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
statistic = NULL
)| (1) | (2) | (3) | |
|---|---|---|---|
| (Intercept) | 3.2 (0.0)*** | 3.4 (0.0)*** | 3.2 (0.0)*** |
| gndrMale | 0.0 (0.0) | −0.1 (0.0)* | 0.1 (0.0)** |
| educ.baNo BA | −0.3 (0.0)*** | −0.1 (0.0)* | |
| lrscaleRight | 0.4 (0.0)*** | ||
| gndrMale × educ.baNo BA | −0.2 (0.0)*** | ||
| Num.Obs. | 55424 | 22699 | 55424 |
| R2 | 0.000 | 0.010 | 0.001 |
| R2 Adj. | 0.000 | 0.010 | 0.001 |
| AIC | 294237.7 | 120029.7 | 294164.9 |
| BIC | 294264.5 | 120069.8 | 294209.5 |
| Log.Lik. | −147115.864 | −60009.840 | −147077.454 |
| RMSE | 2.43 | 2.45 | 2.42 |
interaction_plot <- effect("gndr*educ.ba", model3, na.rm=TRUE)
plot(interaction_plot,
main="Interaction Effect of Gender and Education on Trust in Parliament",
xlab="Education Level",
ylab="Trust in Parliament")
interaction_plot##
## gndr*educ.ba effect
## educ.ba
## gndr BA or more No BA
## Female 3.236436 3.168294
## Male 3.348610 3.077693df$gndr <- factor(df$gndr, levels = c("Male", "Female"))
df$educ.ba <- factor(df$educ.ba, levels = c("No BA", "BA or more"))
ggplot(df, aes(x=educ.ba, y=trstplt, color=gndr)) +
stat_summary(fun=mean, geom="point", size=3) +
stat_summary(fun=mean, geom="line", aes(group=gndr), size=1) +
labs(title="Interaction Effect of Gender and Education on Trust in Parliament",
x="Education Level", y="Mean Trust in Parliament") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Removed 1298 rows containing non-finite values (`stat_summary()`).
## Removed 1298 rows containing non-finite values (`stat_summary()`).