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

# 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.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ 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

## Warning: package 'fst' was built under R version 4.3.2

## Warning: package 'modelsummary' was built under R version 4.3.2

## Version 2.0.0 of `modelsummary`, to be released soon, will introduce a
##   breaking change: The default table-drawing package will be `tinytable`
##   instead of `kableExtra`. All currently supported table-drawing packages
##   will continue to be supported for the foreseeable future, including
##   `kableExtra`, `gt`, `huxtable`, `flextable, and `DT`.
##   
##   You can always call the `config_modelsummary()` function to change the
##   default table-drawing package in persistent fashion. To try `tinytable`
##   now:
##   
##   config_modelsummary(factory_default = 'tinytable')
##   
##   To set the default back to `kableExtra`:
##   
##   config_modelsummary(factory_default = 'kableExtra')

## [[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"

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

## Warning: package 'infer' was built under R version 4.3.3

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

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

getwd()

## [1] "C:/Users/hassa947/Downloads"

ess <- read_fst("All-ESS-Data.fst")

## Warning: package 'fstcore' was built under R version 4.3.2

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

france_data <- ess %>%
  # Filter for France
  filter(cntry == "FR") %>%
  
  # Recoding clsprty and cleaning trstplt (will use later)
  mutate(
    clsprty = case_when(
      clsprty == 1 ~ "Yes",  # Recode 1 to "Yes"
      clsprty == 2 ~ "No",   # Recode 2 to "No"
      clsprty %in% c(7, 8, 9) ~ NA_character_,  # Handle other specific cases where you want to set it as NA
      TRUE ~ as.character(clsprty)  # Keep other values as-is but ensure they are characters
    ),
    trstplt = ifelse(trstplt %in% c(77, 88, 99), NA, trstplt),
    
    # Recoding gender
    gndr = case_when(
      gndr == 1 ~ "Male",
      gndr == 2 ~ "Female",
      gndr == 9 ~ NA_character_,
      TRUE ~ as.character(gndr)
    ),

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

unique(france_data$gndr)

## [1] "Female" "Male"

unique(france_data$clsprty)

## [1] "Yes" "No"  NA

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

# double-checking + quick distribution
unique(france_data$clsprty)

## [1] "Yes" "No"

table(france_data$clsprty)

## 
##   No  Yes 
## 9098 9480

france_clean <- france_data %>% 
  mutate(
    clsprty = ifelse(clsprty == 2, 0, ifelse(clsprty %in% c(7, 8, 9), NA, clsprty)), # If 'clsprty' is 2, set it to 0. If it's 7, 8, or 9, set it to NA.
    stfdem = ifelse(stfdem %in% c(77, 88, 99), NA, stfdem), # For 'stfdem', set values 77, 88, and 99 to NA.
    trstplt = ifelse(trstplt %in% c(77, 88, 99), NA, trstplt), # For 'trstplt', set values 77, 88, and 99 to NA.
  )

france_clean <- france_clean %>% filter(!is.na(clsprty), !is.na(stfdem), !is.na(trstplt))

france_clean <- france_clean %>%
  mutate(
    # Recoding 'lrscale' into a new variable 'ID'
    ID = case_when(
      lrscale %in% 0:3 ~ "Left",                    # Values 0 to 3 in 'lrscale' are categorized as "Left"
      lrscale %in% 7:10 ~ "Right",                  # Values 7 to 10 in 'lrscale' are categorized as "Right"
      lrscale %in% 4:6 ~ "Moderate",                # Values 4 to 6 in 'lrscale' are categorized as "Moderate"
      lrscale %in% c(77, 88, 99) ~ NA_character_    # Values 77, 88, 99 in 'lrscale' are set as NA
    ),
    # Recoding 'yrbrn' (year born) into a new variable 'Gen'
    Gen = case_when(
      yrbrn %in% 1900:1945 ~ "1",                   # Birth years from 1900 to 1945 categorized as "1"
      yrbrn %in% 1946:1964 ~ "2",                   # Birth years from 1946 to 1964 categorized as "2"
      yrbrn %in% 1965:1979 ~ "3",                   # Birth years from 1965 to 1979 categorized as "3"
      yrbrn %in% 1980:1996 ~ "4",                   # Birth years from 1980 to 1996 categorized as "4"
      TRUE ~ as.character(yrbrn)                    # Other years remain unchanged
    ),
    # Converting 'Gen' to a factor with meaningful labels
    Gen = factor(Gen, levels = c("1", "2", "3", "4"),
                 labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials")) # Assigning labels to the levels in 'Gen'
  )

tab_dat <- france_data

# Create a crosstable and save it into the object mytab
mytab <- table(tab_dat$clsprty, 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.1099746 0.3797444
## [2,] 0.1145921 0.3956889

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

## Table of Expected Frequencies:

print(ftab)

##          [,1]     [,2]
## [1,] 2043.108 7054.892
## [2,] 2128.892 7351.108

# 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

#interpreting the critical value: at the 0.05 significance level, we reject the null hypothesis if the chi-squareed test statistic, which we calculate from the chisq, test function is greater than 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: 160.4057

#interpreation of chi-square: a value of 160.4057 suggests a more significant difference between the actual and preddicted frequencies, therefore making the evidence strogner to reject the null hypothesis.

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

## p-value: 0

#interpreation of p-value: Since our p-value is smaller than 0.05, it suggests strong evidence against the null hypothesis, therefore we reject it.

# 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 = 160.41, df = 1, p-value < 2.2e-16

#interpretation: the p-value denotes a small number almost equal to 0. This small p-value indicates a substantial correlation betwen the two categorical variables examined. This terefore prompts us to reject the null hypothesis. This is further supported by the large chi-squared test statistic of 160.41. an increased chi-square statistic suggests a more significant difference between the actual and preddicted frequencies, therefore making the evidence even strogner to reject the null hypothesis.

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

test_stat <- france_data %>%
  specify(explanatory = domicil, # change variable name for explanatory variable
          response = sbbsntx) %>% # change variable name for outcome of interest
  hypothesize(null = "independence") %>%
  calculate(stat = "slope") # replace in between quotation marks appropriate test statistic

## Warning: Removed 14528 rows containing missing values.

## Message: The independence null hypothesis does not inform calculation of the
## observed statistic (a slope) and will be ignored.

print(test_stat$stat)

##    domicil 
## -0.0768133

#Interpretation: A slope of -0.0788 suggests that there isn’t enough informstion about the relationship between the variables. The negavtive number represents a negative relationship between the two variables. If one variable increases, the other decreases.

null_distribution <- france_data %>%
  specify(explanatory = domicil,
          response = sbbsntx) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>% # only line we are adding, use reps = 1000 as standard but change if too computationally demanding
  calculate(stat = "slope")

## Warning: Removed 14528 rows containing missing values.

p_val <- null_distribution %>% # replace name here if you assigned something other than null_distribution above
  get_pvalue(obs_stat = test_stat, direction = "two-sided") # would only replace test_stat if assigned another name in Step 1

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

#interpreation of p-value: Since our p-value is smaller than 0.05, it suggests strong evidence against the null hypothesis, therefore we reject it.

#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")

null_distribution

## Response: sbbsntx (numeric)
## Explanatory: domicil (numeric)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate      stat
##        <int>     <dbl>
##  1         1 -0.00183 
##  2         2 -0.00406 
##  3         3 -0.00276 
##  4         4 -0.0328  
##  5         5  0.0288  
##  6         6 -0.00517 
##  7         7 -0.000345
##  8         8 -0.0295  
##  9         9 -0.0176  
## 10        10  0.00652 
## # ℹ 990 more rows

#Interpreation: The values here indicate clear evidence, under assumptions of frequentist hypothesis testing, to reject the null. If the observed statistic would fall within the null distribution, we would see shaded red regions of the null distribution that are as extreme than our observed statistic

#visualizing using confidence intervals

conf_int <- null_distribution %>%
  get_confidence_interval(level = 0.95, type = "percentile")


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

null_distribution

## Response: sbbsntx (numeric)
## Explanatory: domicil (numeric)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate      stat
##        <int>     <dbl>
##  1         1 -0.00183 
##  2         2 -0.00406 
##  3         3 -0.00276 
##  4         4 -0.0328  
##  5         5  0.0288  
##  6         6 -0.00517 
##  7         7 -0.000345
##  8         8 -0.0295  
##  9         9 -0.0176  
## 10        10  0.00652 
## # ℹ 990 more rows

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

table(france_data$ccrdprs)

## 
##   0   1   2   3   4   5   6   7   8   9  10  66  88  99 
##  70  25  60  74 107 468 391 714 928 355 692  25  23   1

table(france_data$lrscale)

## 
##    0    1    2    3    4    5    6    7    8    9   10   77   88 
##  930  493 1253 1974 1806 5167 1560 1723 1439  401  676  313  843

france_data <- france_data %>%
  mutate(
    ccrdprs = ifelse(ccrdprs %in% c(66:99), NA, ccrdprs)
  )

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

france_data <- france_data %>%
  mutate(
    lrscale = factor(lrscale)
  )

france_data <- france_data %>%
    mutate(lrscale = ifelse(lrscale %in% 77:99, NA, lrscale)) %>%
    filter(!lrscale %in% 4:6) %>%
  mutate(
    lrscale = case_when(
      lrscale %in% 0:3 ~ "Left",
      lrscale %in% 7:10 ~ "Right",
      lrscale %in% c(77, 88, 99) ~ NA_character_
    ))

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

table(france_data$ccrdprs)

## 
##   0   1   2   3   4   5   6   7   8   9  10 
##  23   8  22  32  50 182 161 308 418 145 257

table(france_data$lrscale)

## 
##  Left Right 
##   471  1135

test_stat <- france_data %>%
  specify(explanatory = lrscale, 
          response = ccrdprs) %>% 
  hypothesize(null = "independence") %>%
  calculate(stat = "t", order = c("Left", "Right"
  ))

print(test_stat$stat)

##        t 
## 4.466582

#interpreation: The t-statistic measures the difference between the observed sample mean and the population mean in units of standard error. The larger the magnitude of t, the greater the evidence against the null hypothesis. A t-statistic of 4.35 is not substantially large indicating a smaller difference between the two variables ccrdrps and lrscale. However our p-value can give us more evidence on whether to reject the null hypothesis or not.

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

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

#Intepretation of p-value: Since our p-value is smaller than 0.05, it suggests strong evidence against the null hypothesis, therefore we reject it.

#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")

null_distribution

## Response: ccrdprs (numeric)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate    stat
##        <int>   <dbl>
##  1         1 -0.910 
##  2         2  2.67  
##  3         3  0.859 
##  4         4 -0.479 
##  5         5 -0.150 
##  6         6 -0.0986
##  7         7  0.378 
##  8         8 -1.74  
##  9         9  0.966 
## 10        10  0.322 
## # ℹ 990 more rows

#visualizing using confidence intervals

conf_int <- null_distribution %>%
  get_confidence_interval(level = 0.95, type = "percentile")


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

null_distribution

## Response: ccrdprs (numeric)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate    stat
##        <int>   <dbl>
##  1         1 -0.910 
##  2         2  2.67  
##  3         3  0.859 
##  4         4 -0.479 
##  5         5 -0.150 
##  6         6 -0.0986
##  7         7  0.378 
##  8         8 -1.74  
##  9         9  0.966 
## 10        10  0.322 
## # ℹ 990 more rows

Hassan_Aleena_Homework_3

2024-03-09

R Markdown

Including Plots