Harang_Dagmar_Homework

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.

rm(list=ls()); gc()

##           used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells  545962 29.2    1225182 65.5         NA   669417 35.8
## Vcells 1014116  7.8    8388608 64.0      16384  1851733 14.2

# List of packages
packages <- c("tidyverse", "infer", "fst") # 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] "infer"     "lubridate" "forcats"   "stringr"   "dplyr"     "purrr"    
##  [7] "readr"     "tidyr"     "tibble"    "ggplot2"   "tidyverse" "stats"    
## [13] "graphics"  "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[3]]
##  [1] "fst"       "infer"     "lubridate" "forcats"   "stringr"   "dplyr"    
##  [7] "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"   "tidyverse"
## [13] "stats"     "graphics"  "grDevices" "utils"     "datasets"  "methods"  
## [19] "base"

france_data <- read.fst("france_data.fst")

When you are done, click on “Knit” above, then “Knit to Html”. Wait for everything to compile. If you get an error like “Execution halted”, it means there are issues with your code you must fix. When all issues are fixed, it will prompt a new window. Then click on “Publish” in the top right, and then Rpubs (the first option) and follow the instructions to create your Rpubs account and get your Rpubs link for your document (i.e., html link as I provide for the tutorial).

Note: Make sure to provide both your markdown file and R pubs link. If you do not submit both, you will be penalized 2 pts. out of the 5 pts. total.

Task 1 Based on the recoding of lrscale for “left” and “right” and omitting “moderates” (see Tutorial 4), and educ.ba from this tutorial, do the long form coding for Chi-Square. Using the steps outlined in the tutorial, first generate the tables of expected proportions and frequencies. Determine and interpret the critical value for independence. Finally, determine and interpret both the Pearson’s chi-squared statistic and the p-value. Discuss main takeaways.

unique(france_data$lrscale)

##  [1]  3  2  1  4  5  6 10  8 88  7  9  0 77

unique(france_data$educ.ba)

## NULL

france_data <- france_data %>% filter(!is.na(lrscale))


# double-checking + quick distribution
unique(france_data$educ.ba)

## NULL

table(france_data$educ.ba)

## < table of extent 0 >

# Now you can proceed with your data manipulation code
france_data <- france_data %>%
  # Filter for France
  filter(cntry == "FR") %>%
  
  # Recoding clsprty and cleaning trstplt (will use later)
  mutate(
    lrscale = case_when(
      lrscale == 1 ~ "Yes",  
      lrscale == 2 ~ "No",   
      lrscale %in% c(7, 8, 9) ~ NA_character_,  
      TRUE ~ as.character(lrscale)  
    ),

  
    # Recoding education with case_when for consistency
    educ.ba = case_when(
      essround < 5 & edulvla == 5 ~ "BA or more",
      essround >= 5 & edulvlb > 600 ~ "BA or more",
      TRUE ~ "No BA"
    )
  ) %>%
  # Ensure all NAs are handled uniformly
  mutate(
    edulvla = ifelse(edulvla %in% c(77, 88), NA_integer_, edulvla),
    edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
    educ.ba = as.character(educ.ba)  # Force educ.ba to be treated as character
  )

tab_dat <- france_data

# Create a crosstable and save it into the object mytab
mytab <- table(tab_dat$lrscale, tab_dat$educ.ba)

# Calculate row-wise sums of counts in table 
rsums <- rowSums(mytab)

# Calculate column-wise sums of counts in table
csums <- colSums(mytab)

# Get the total count of the table (N)
N <- sum(mytab)

# Generate table of expected proportions
ptab <- tcrossprod(rsums/N, csums/N)
cat("Table of Expected Proportions:\n")

## Table of Expected Proportions:

print(ptab)

##              [,1]       [,2]
##  [1,] 0.013715802 0.04756030
##  [2,] 0.009956656 0.03452525
##  [3,] 0.028970095 0.10045540
##  [4,] 0.026676871 0.09250352
##  [5,] 0.076547236 0.26543175
##  [6,] 0.023135437 0.08022340
##  [7,] 0.005529863 0.01917510
##  [8,] 0.013643231 0.04730866
##  [9,] 0.018418362 0.06386669
## [10,] 0.007242524 0.02511385

ftab <- N * ptab
cat("Table of Expected Frequencies:\n")

## Table of Expected Frequencies:

print(ftab)

##             [,1]      [,2]
##  [1,]  211.52509  733.4749
##  [2,]  153.55155  532.4485
##  [3,]  446.77681 1549.2232
##  [4,]  411.41071 1426.5893
##  [5,] 1180.51148 4093.4885
##  [6,]  356.79471 1237.2053
##  [7,]   85.28155  295.7185
##  [8,]  210.40591  729.5941
##  [9,]  284.04798  984.9520
## [10,]  111.69420  387.3058

# Critical Value for Independence:
alpha <- 0.05 # significance level
c_val <- qchisq(alpha, df = 1, lower.tail = FALSE)
cat("Critical Value:", round(c_val, 3), "\n")

## Critical Value: 3.841

# Pearson's chi-squared statistic:
test_stat <- sum((mytab - ftab)^2 / ftab)
cat("Pearson's X^2:", round(test_stat, 4), "\n")

## Pearson's X^2: 470.4402

# p-value:
p_val <- pchisq(test_stat, df = 1, lower.tail = F)
cat("p-value:", round(p_val, 4), "\n")

## p-value: 0

# Using chisq.test for validation:
cat("Chi-squared test result:\n")

## Chi-squared test result:

print(chisq.test(mytab, correct = F))

## 
##  Pearson's Chi-squared test
## 
## data:  mytab
## X-squared = 470.44, df = 9, p-value < 2.2e-16

answer: The difference between the observed data and what we’d expect under the assumption of independence is statistically significant. Thus, there is a statistically significant association between being on the left or right side of the political scale and one’s educational level. This means that the probability of feeling close to a party is not the same across different education levels.The extremely small p-value (essentially 0) means that we reject the null hypothesis that the two variables are independent of each other. In other words, there’s strong evidence to suggest that there’s a significant association between the two variables in mytab.The large chi-squared statistic of 481.47 further reinforces this point. If the two variables were independent, we’d expect a chi-squared statistic closer to 0.

Task 2 Do Steps 1 to 3 using the infer package (as we did in the tutorial), for the variables sbbsntx (recoded to the corresponding 5 categories) as the response and domicil (recoded to urban and rural, setting 2 to urban instead of peri-urban) as the explanatory variable for the country of France.

See variable info here: https://ess.sikt.no/en/datafile/ffc43f48-e15a-4a1c-8813-47eda377c355/92?tab=1&elements=[%22874a86e2-a0f6-40b4-aef7-51c8eed98a7d/1%22])

Provide interpretations of the output (consider the variable info).

france_data <- france_data %>%
  filter(cntry == "FR") %>%
  mutate(
    sbbsntx = case_when(
      sbbsntx == 1 ~ "Agree strongly",
      sbbsntx == 2 ~ "Agree",
      sbbsntx == 3 ~ "Neither agree nor disagree",
      sbbsntx == 4 ~ "Disagree",
      sbbsntx == 5 ~ "Disagree strongly",
      sbbsntx %in% c(7, 8, 9) ~ NA_character_,
      TRUE ~ as.character(sbbsntx)
    ),
    domicil = case_when(
      domicil == 1 ~ "Big city",
      domicil == 2 ~ "Suburbs",
      domicil == 3 ~ "Town",
      domicil == 4 ~ "Village",
      domicil == 5 ~ "Farm",
      domicil %in% c(7, 8, 9) ~ NA_character_,
      TRUE ~ as.character(domicil)
    )
  )

test_stat <- france_data %>%
  specify(explanatory = domicil,
          response = sbbsntx) %>% 
  hypothesize(null = "independence") %>%
  calculate(stat = "Chisq")

## Warning: Removed 14999 rows containing missing values.

print(test_stat$stat)

## X-squared 
##  36.01302

null_distribution <- france_data %>%
  specify(explanatory = domicil,
          response = sbbsntx) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "Chisq")

## Warning: Removed 14999 rows containing missing values.

p_val <- null_distribution %>% # 
  get_pvalue(obs_stat = test_stat, direction = "two-sided") 

p_val

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

answer: Considering that the p-value stands at 0.006, it suggests that if the null was true, we would reject the null (because the sample is less than than 0.05)

Task 3 Now do Step 4 (i.e., the null distribution visualization with confidence intervals) for the same variables and country as Task 2, and interpret the output.

null_distribution %>%
  visualize() +
  shade_p_value(obs_stat = test_stat, direction = "two-sided")

## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.

Answer:

Task 4 Conduct Steps 1 to 3 using the infer package (as we did in the tutorial), for the variables ccrdprs (leaving it on the 0-10 numeric scale) as the response (or outcome) variable and lrscale (recoded as left and right, omitting “moderates” as we did in Task 1) as the explanatory variable.

Variable info here: https://ess.sikt.no/en/datafile/ffc43f48-e15a-4a1c-8813-47eda377c355/92?tab=1&elements=[%2283b08cf2-508e-49a0-9fc8-c3cc281290ec/4%22]

Provide interpretations of the output (consider the variable info).

france_data <- france_data %>%
  filter(cntry == "FR") %>%
  mutate(
    lrscale = case_when(
      lrscale == 1 ~ "Yes",  
      lrscale == 2 ~ "No",   
      lrscale %in% c(7, 8, 9) ~ NA_character_,  
      TRUE ~ as.character(lrscale)  
    ),

    ccrdprs = case_when(
      ccrdprs == 0 ~ "Not at all",
      ccrdprs == 1 ~ NA_character_,
      ccrdprs == 2 ~ NA_character_,
      ccrdprs == 3 ~ NA_character_,
      ccrdprs == 4 ~ NA_character_,
      ccrdprs == 5 ~ NA_character_,
      ccrdprs == 6 ~ NA_character_,
      ccrdprs == 7 ~ NA_character_,
      ccrdprs == 8 ~ NA_character_,
      ccrdprs == 9 ~ NA_character_,
      ccrdprs == 10 ~ "A great deal",
      ccrdprs %in% c(7, 8, 9) ~ NA_character_,
      TRUE ~ as.character(ccrdprs)
    )
  )

test_stat <- france_data %>%
  specify(explanatory = lrscale,
          response = ccrdprs) %>% 
  hypothesize(null = "independence") %>%
  calculate(stat = "Chisq")

## Warning: Removed 18333 rows containing missing values.

print(test_stat$stat)

## X-squared 
##  100.0201

null_distribution <- france_data %>%
  specify(explanatory = lrscale,
          response = ccrdprs) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "Chisq")

## Warning: Removed 18333 rows containing missing values.

p_val <- null_distribution %>% # 
  get_pvalue(obs_stat = test_stat, direction = "two-sided")

## 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.

p_val

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

Answer: Considering that the p-value stands at 0.004, it suggests that if the null was true, we would reject the null (because the sample is less than than 0.05)

Task 5 Now do Step 4 (i.e., the null distribution visualization with confidence intervals) for the same variables and country as Task 4, and interpret the output.

null_distribution %>%
  visualize() +
  shade_p_value(obs_stat = test_stat, direction = "two-sided")

## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.

Harang_Dagmar_Homework_3

2024-03-11

R Markdown

Including Plots