R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

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.5
## ✔ 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"
getwd()
## [1] "C:/Users/User/Downloads"
ess <- read_fst("All-ESS-Data.fst")
unique(ess$cntry)
##  [1] "AT" "BE" "CH" "CZ" "DE" "DK" "ES" "FI" "FR" "GB" "GR" "HU" "IE" "IL" "IT"
## [16] "LU" "NL" "NO" "PL" "PT" "SE" "SI" "EE" "IS" "SK" "TR" "UA" "BG" "CY" "RU"
## [31] "HR" "LV" "RO" "LT" "AL" "XK" "ME" "RS" "MK"

#TASK1 #Prompt and question: calculate the average for the variable ‘happy’ for the country of Norway. On average, based on the ESS data, who reports higher levels of happiness: Norway or Belgium? #based on the ESS data norway reports higher level of happiness because the mean for norway is 7.9 whereas for belgium the mean is 7.7

belgium_happy <- ess %>% # note: if you work from belgium_data replace "ess" with belgium_data
  filter(cntry == "BE") %>% 
  select(happy)
belgium_happy$y <- belgium_happy$happy

table(belgium_happy$y)
## 
##    0    1    2    3    4    5    6    7    8    9   10   77   88   99 
##   50   27  104  194  234  830  999 3503 6521 3402 1565    3   16    3
# need to remove 77, 88, 99 or else will alter results. See data portal for what they represent (e.g. DK, Refusal, etc.)

# Recode values 77 through 99 to NA
belgium_happy$y[belgium_happy$y %in% 77:99] <- NA

# checking again
table(belgium_happy$y)
## 
##    0    1    2    3    4    5    6    7    8    9   10 
##   50   27  104  194  234  830  999 3503 6521 3402 1565

###Calculating Mean and Median

mean_y <- mean(belgium_happy$y, na.rm = TRUE)
cat("Mean of 'y' is:", mean_y, "\n")
## Mean of 'y' is: 7.737334
# Now, let's determine the median of 'y'.
median_y <- median(belgium_happy$y, na.rm = TRUE)
cat("Median of 'y' is:", median_y, "\n")
## Median of 'y' is: 8
# Using tidyverse syntax principles, we can get both the mean and median at once.
belgium_happy %>%
  summarize(
    mean_y = mean(y, na.rm = TRUE),
    median_y = median(y, na.rm = TRUE)
  ) %>%
  print()
##     mean_y median_y
## 1 7.737334        8

###finding the mode and standard deviation

mode_y <- belgium_happy %>%
  filter(!is.na(y)) %>%
  count(y) %>%
  arrange(desc(n)) %>%
  slice(1) %>%
  pull(y)

cat("\nMode of Y:", mode_y, "\n")
## 
## Mode of Y: 8
sd_y <- sd(belgium_happy$y, na.rm = TRUE)
cat("Standard Deviation of 'y':", sd_y, "\n")
## Standard Deviation of 'y': 1.52045
# Let's compare if had not cleaned properly
belgium_happy2 <- ess %>%
  filter(cntry == "BE") %>%
  select(happy)

sd_happy2 <- sd(belgium_happy2$happy, na.rm = TRUE)
cat("Standard Deviation of 'happy not cleaned':", sd_happy2, "\n")
## Standard Deviation of 'happy not cleaned': 3.234511

###the summary function

summary(belgium_happy)
##      happy              y         
##  Min.   : 0.000   Min.   : 0.000  
##  1st Qu.: 7.000   1st Qu.: 7.000  
##  Median : 8.000   Median : 8.000  
##  Mean   : 7.839   Mean   : 7.737  
##  3rd Qu.: 9.000   3rd Qu.: 9.000  
##  Max.   :99.000   Max.   :10.000  
##                   NA's   :22

###NORWAY

ess <- read_fst("All-ESS-Data.fst")
unique(ess$cntry)
##  [1] "AT" "BE" "CH" "CZ" "DE" "DK" "ES" "FI" "FR" "GB" "GR" "HU" "IE" "IL" "IT"
## [16] "LU" "NL" "NO" "PL" "PT" "SE" "SI" "EE" "IS" "SK" "TR" "UA" "BG" "CY" "RU"
## [31] "HR" "LV" "RO" "LT" "AL" "XK" "ME" "RS" "MK"
norway_happy <- ess %>% # note: if you work from norway_data replace "ess" with norway_data
  filter(cntry == "NO") %>% 
  select(happy)
norway_happy$y <- norway_happy$happy

table(norway_happy$y)
## 
##    0    1    2    3    4    5    6    7    8    9   10   77   88 
##   15   29   59  163  238  730  817 2617 5235 3796 2344   12   10
# need to remove 77, 88, 99 or else will alter results. See data portal for what they represent (e.g. DK, Refusal, etc.)

# Recode values 77 through 99 to NA
norway_happy$y[norway_happy$y %in% 77:99] <- NA

# checking again
table(norway_happy$y)
## 
##    0    1    2    3    4    5    6    7    8    9   10 
##   15   29   59  163  238  730  817 2617 5235 3796 2344

###Calculating Mean and Mediam NORWAY

mean_y <- mean(norway_happy$y, na.rm = TRUE)
cat("Mean of 'y' is:", mean_y, "\n")
## Mean of 'y' is: 7.975005
#Now, let's determine the median of 'y'.
median_y <- median(norway_happy$y, na.rm = TRUE)
cat("Median of 'y' is:", median_y, "\n")
## Median of 'y' is: 8
# Using tidyverse syntax principles, we can get both the mean and median at once.
norway_happy %>%
  summarize(
    mean_y = mean(y, na.rm = TRUE),
    median_y = median(y, na.rm = TRUE)
  ) %>%
  print()
##     mean_y median_y
## 1 7.975005        8

###Finding the Mode and Standard Deviation NORWAY

mode_y <- norway_happy %>%
  filter(!is.na(y)) %>%
  count(y) %>%
  arrange(desc(n)) %>%
  slice(1) %>%
  pull(y)

cat("\nMode of Y:", mode_y, "\n")
## 
## Mode of Y: 8
sd_y <- sd(norway_happy$y, na.rm = TRUE)
cat("Standard Deviation of 'y':", sd_y, "\n")
## Standard Deviation of 'y': 1.539186
# Let's compare if had not cleaned properly
norway_happy2 <- ess %>%
  filter(cntry == "NO") %>%
  select(happy)

sd_happy2 <- sd(norway_happy2$happy, na.rm = TRUE)
cat("Standard Deviation of 'happy not cleaned':", sd_happy2, "\n")
## Standard Deviation of 'happy not cleaned': 3.146621

###The Summary Function NORWAY

summary(norway_happy)
##      happy              y         
##  Min.   : 0.000   Min.   : 0.000  
##  1st Qu.: 7.000   1st Qu.: 7.000  
##  Median : 8.000   Median : 8.000  
##  Mean   : 8.076   Mean   : 7.975  
##  3rd Qu.: 9.000   3rd Qu.: 9.000  
##  Max.   :88.000   Max.   :10.000  
##                   NA's   :22

#TASK2 #Prompt and question: what is the most common category selected, for Irish respondents, for frequency of binge drinking? The variable of interest is: alcbnge. #the most common catergory selcted for Irish respondents is weekly for frequency of binge drinking ###counting categorical variables

ireland_bdf12 <- ess %>%
  filter(cntry == "IE") %>%
  select(alcbnge)

ireland_bdf12$y <- ireland_bdf12$alcbnge

table(ireland_bdf12$y)
## 
##   1   2   3   4   5   6   7   8 
##  65 650 346 417 239 641  26   6
# Recode values 6 through 8 to NA
ireland_bdf12$y[ireland_bdf12$y %in% 6:8] <- NA
# Compute mean and median
mean_y <- mean(ireland_bdf12$y, na.rm = TRUE)
median_y <- median(ireland_bdf12$y, na.rm = TRUE)

cat("Mean of 'y':", mean_y, "\n")
## Mean of 'y': 3.066977
cat("Median of 'y':", median_y, "\n")
## Median of 'y': 3
# Converting to categories to get mode as a category instead of a number
df <- ireland_bdf12 %>%
  mutate(
    y_category = case_when(
      y == 1 ~ "Daily or almost daily",
      y == 2 ~ "Weekly",
      y == 3 ~ "Monthly",
      y == 4 ~ "Less than monthly",
      y == 5 ~ "Never",
      TRUE ~ NA_character_
    ),
    y_category = fct_relevel(factor(y_category),  ### here you would put the categories in order you want them to appear or else it will appear alphabetically
                             "Daily or almost daily", 
                             "Weekly", 
                             "Monthly", 
                             "Less than monthly", 
                             "Never")
  )

# To confirm the conversion:
table(df$y_category)
## 
## Daily or almost daily                Weekly               Monthly 
##                    65                   650                   346 
##     Less than monthly                 Never 
##                   417                   239
# Let's determine the mode of our newly created category:

get_mode <- function(v) {
  tbl <- table(v)
  mode_vals <- as.character(names(tbl)[tbl == max(tbl)])
  return(mode_vals)
}

mode_values <- get_mode(df$y_category)
cat("Mode of y category:", paste(mode_values, collapse = ", "), "\n")
## Mode of y category: Weekly

#TASK3 #Prompt and question: when you use the summary() function for the variable plnftr (about planning for future or taking every each day as it comes from 0-10) for both the countries of Portugal and Serbia, what do you notice? What stands out as different when you compare the two countries (note: look up the variable information on the ESS website to help with interpretation)? Explain while referring to the output generated. # I notice the difference between the mean and median for each country. Portugal has a mean of 6.4 whereas Serbia has 4.9. The median for Portugal is 5 and for Serbia its 4. You can see a differences in the stats ###Comparing means with groupby

result <- ess %>%
  # Step 1: Filter for the countries of interest
  filter(cntry %in% c("RS", "PT")) %>%
  
  # Step 2: Recode 77, 88, 99 to NA for 'plnftr'
  mutate(plnftr = recode(plnftr, `77` = NA_real_, `88` = NA_real_, `99` = NA_real_)) %>%
  
  # Step 3: Compute mean by country
  group_by(cntry) %>%
  summarize(mean_plnftr = mean(plnftr, na.rm = TRUE))

print(result)
## # A tibble: 2 × 2
##   cntry mean_plnftr
##   <chr>       <dbl>
## 1 PT           5.42
## 2 RS           4.14
portugal_data <- filter(ess, cntry == "PT")
serbia_data <- filter(ess, cntry == "RS")
summary_portugal <- summary(portugal_data$plnftr)
summary_serbia <- summary(serbia_data$plnftr)
print("Summary for Portugal:")
## [1] "Summary for Portugal:"
print(summary_portugal)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   3.000   5.000   6.426   8.000  88.000   14604
print("Summary for Serbia:")
## [1] "Summary for Serbia:"
print(summary_serbia)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   4.000   4.983   8.000  88.000    1505

#average of Y by X #task4

italy_data <- ess %>% 
  filter(cntry == "IT")

# Convert gender and stfdem (representing pol ID self-placement from left to right)
italy_data <- italy_data %>%
  mutate(
    gndr = case_when(
      gndr == 1 ~ "Male",
      gndr == 2 ~ "Female",
      TRUE ~ as.character(gndr)
    ),
    stfdem = ifelse(stfdem %in% c(77, 88), NA, stfdem)  # Convert stfdem values
  )
# Compute mean for male
mean_male_stfdem <- italy_data %>%
  filter(gndr == "Male") %>%
  summarize(mean_stfdem_men = mean(stfdem, na.rm = TRUE))

print(mean_male_stfdem)
##   mean_stfdem_men
## 1        4.782646
# Compute mean for male
mean_Female_stfdem <- italy_data %>%
  filter(gndr == "Female") %>%
  summarize(mean_stfdem_Women = mean(stfdem, na.rm = TRUE))

print(mean_Female_stfdem)
##   mean_stfdem_Women
## 1          4.694652
# Compute average of stfdem by gender
means_by_gender <- italy_data %>%
  group_by(gndr) %>% # here you are "grouping by" your second variable
  summarize(stfdem = mean(stfdem, na.rm = TRUE)) # here you are summarizing your variable of interest

print(means_by_gender)
## # A tibble: 3 × 2
##   gndr   stfdem
##   <chr>   <dbl>
## 1 9        3.25
## 2 Female   4.69
## 3 Male     4.78

###boxplots #TASK5 #Prompt: Interpret the boxplot graph of stfedu and stfhlth that we generated already: according to ESS data, would we say that the median French person is more satisfied with the education system or health services? Explain. #education system

getwd()
## [1] "C:/Users/User/Downloads"
ess <- read_fst("All-ESS-Data.fst")
france_data <- ess
france_data %>%
  # Setting values to NA
  mutate(stfedu = ifelse(stfedu %in% c(77, 88, 99), NA, stfedu),
         stfhlth = ifelse(stfhlth %in% c(77, 88, 99), NA, stfhlth)) %>%
  # Reshaping the data
  select(stfedu, stfhlth) %>%
  gather(variable, value, c(stfedu, stfhlth)) %>%
  # Creating the boxplot
  ggplot(aes(x = variable, y = value)) +
  geom_boxplot() +
  labs(y = "Y-axis", x = "X-axis", title = "Boxplot of satisfaction with the state of education vs. health services") + theme_minimal()
## Warning: Removed 29886 rows containing non-finite values (`stat_boxplot()`).

france_data %>%
  # Setting values to NA
  mutate(stfedu = ifelse(stfedu %in% c(77, 88, 99), NA, stfedu),
         stfhlth = ifelse(stfhlth %in% c(77, 88, 99), NA, stfhlth)) %>%
  # Reshaping the data
  select(stfedu, stfhlth) %>%
  gather(variable, value, c(stfedu, stfhlth)) %>%
  # Creating the boxplot
  ggplot(aes(y = value)) +
  geom_boxplot() +
  labs(y = "Y-axis", title = "Boxplot of satisfaction with the state of education vs. health services") +
  theme_minimal()
## Warning: Removed 29886 rows containing non-finite values (`stat_boxplot()`).