Brasil_Layla_Homework_3
2024-03-11
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.00Including 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", "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.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] "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"ess <- read_fst("All-ESS-Data.fst")france_data <- read.fst("france_data.fst")Homework Start
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.
france_data <- ess %>%
# Filter for France
filter(cntry == "FR") %>%
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% c(77, 88, 99) ~ NA_character_ # Values 77, 88, 99 in 'lrscale' are set as NA
),
# 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
)france_data <- france_data %>% filter(!is.na(ID))table(france_data$ID)##
## Left Right
## 4709 4302table(france_data$educ.ba)##
## BA or more No BA
## 2261 6750unique(france_data$ID)## [1] "Left" "Right"unique(france_data$educ.ba)## [1] "No BA" "BA or more"france_data <- france_data %>% filter(!is.na(educ.ba))## Not necessary, but always good to rename when doing operations so you can backtrack to clean dataset
tab_dat <- france_data
# Create a crosstable and save it into the object mytab
mytab <- table(tab_dat$ID, 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.1311243 0.3914592
## [2,] 0.1197912 0.3576253# Table of expected frequencies
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2]
## [1,] 1181.561 3527.439
## [2,] 1079.439 3222.561# 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: 100.8551Calculating the p-value
# 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 = 100.86, df = 1, p-value < 2.2e-16Interpretation: I calculated a critical value of 3.841, the threshold for my test statisitc. I also calculated the chi-squared test statistic as 100.8551 which measures how much the observed data deviates from the data expected if education and political positioning were independent. Since the x-squared value of 100.8551 is larger than the critical value of 3.841 it means education level (BA or not) and political positioning (left or right side) are not independent at the 0.05 significant level. This means that the probability of politically leaning towards the left or right side is not the same across different educational levels (BA or no BA). The extremely small p value of < 2.2e-16 is lower than 0.05 there is strong evidence to reject the null hypothesis. In my data, there is an association between political positioning and highest educational level achieved.
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).
table(france_data$sbbsntx)##
## 1 2 3 4 5 7 8
## 389 681 350 370 202 2 32france_data <- france_data %>%
mutate(
sbbsntxrecoded = recode(as.character(sbbsntx),
'1' = "Strongly Agree",
'2' = "Agree",
'3' = "Neither disagree or agree",
'4' = "Disagree",
'5' = "Strongly disagree",
'7' = NA_character_,
'8' = NA_character_,
'9' = NA_character_)
) %>%
filter(!is.na(sbbsntxrecoded)) # Removing rows with NA in clsprty or geotable(france_data$sbbsntxrecoded)##
## Agree Disagree Neither disagree or agree
## 681 370 350
## Strongly Agree Strongly disagree
## 389 202france_data <- france_data %>%
mutate(
geo = recode(as.character(domicil),
'1' = "Urban",
'2' = "Urban",
'3' = "Rural",
'4' = "Rural",
'5' = "Rural",
'7' = NA_character_,
'8' = NA_character_,
'9' = NA_character_)
) %>%
filter(!is.na(geo)) # Removing rows with NA in clsprty or geotable(france_data$geo)##
## Rural Urban
## 1376 616Step 1
test_stat <- france_data %>%
specify(explanatory = geo, # change variable name for explanatory variable
response = sbbsntxrecoded) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "Chisq") # replace in between quotation marks appropriate test statistic
print(test_stat$stat)## X-squared
## 17.31799step 2
null_distribution_one <- france_data %>%
specify(explanatory = geo,
response = sbbsntxrecoded) %>%
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 = "Chisq")Step 3
p_val <- null_distribution_one %>% # 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
p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0.008Answer:I calculated a chi-squared test statistic of 11.16404 which is high, suggesting there might be observed differences between where you live (urban or rural) and whether one thinks social benefits/services cost businesses too much in taxes/charges. A relationship is suggested between these two variables. Furthermore, the p value is under 0.05 which makes it statistically significant, and gives evidence to reject the null hypothesis.
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_one %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided")## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
null_distribution_one## Response: sbbsntxrecoded (factor)
## Explanatory: geo (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 1.89
## 2 2 7.50
## 3 3 2.35
## 4 4 0.862
## 5 5 6.49
## 6 6 4.76
## 7 7 5.94
## 8 8 6.45
## 9 9 1.53
## 10 10 7.56
## # ℹ 990 more rowsconf_int <- null_distribution_one %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_distribution_one %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int)## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
null_distribution_one## Response: sbbsntxrecoded (factor)
## Explanatory: geo (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 1.89
## 2 2 7.50
## 3 3 2.35
## 4 4 0.862
## 5 5 6.49
## 6 6 4.76
## 7 7 5.94
## 8 8 6.45
## 9 9 1.53
## 10 10 7.56
## # ℹ 990 more rowsAnswer: the test statistic seems to fall outside the null distribution given the red line’s position in the rejection area. Given the graph it is made clear that we can reject the null hypothesis. This means from this sample we have evidence to prove a significant relationship between between where you live (urban or rural) and whether one thinks social benefits/services cost businesses too much in taxes/charges.
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 %>%
mutate(
ccrdprs = ifelse(ccrdprs %in% c(66, 77, 88, 99), NA, ccrdprs),
)
france_data <- france_data %>% filter(!is.na(ccrdprs))table(france_data$ccrdprs)##
## 0 1 2 3 4 5 6 7 8 9 10
## 26 9 21 21 34 123 94 186 231 90 151test_stat_two <- france_data %>%
specify(explanatory = ID, # change variable name for explanatory variable
response = ccrdprs) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "t") # replace in between quotation marks appropriate test statistic## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "Left" - "Right", or divided in the order "Left" / "Right" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("Left", "Right")` to the calculate() function.print(test_stat_two$stat)## t
## 4.855399Step 2: Simulate the null distribution
null_distribution <- france_data %>%
specify(explanatory = ID,
response = ccrdprs) %>%
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 = "t")## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "Left" - "Right", or divided in the order "Left" / "Right" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("Left", "Right")` to the calculate() function.Step 3: Calculate the p-value of your sample
p_val <- null_distribution %>% # replace name here if you assigned something other than null_distribution above
get_pvalue(obs_stat = test_stat_two, 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 0The t-statistic measures the difference between the observed sample mean and the population mean which helps provide support to reject the null hypothesis. I calculated a t-statistic of 4.855399 which is quite low which means there is must be little difference between ones political position (left or right) and the extent one feels personal responsibility to reduce climate change. The P value is 0 which also means the test statistic is statistically significant.
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.
Step 4: Visualize
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat_two, direction = "two-sided")
null_distribution## Response: ccrdprs (numeric)
## Explanatory: ID (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 -0.125
## 2 2 -0.318
## 3 3 0.728
## 4 4 1.06
## 5 5 0.480
## 6 6 0.0666
## 7 7 -1.41
## 8 8 -0.290
## 9 9 -0.154
## 10 10 -1.83
## # ℹ 990 more rowsconf_int <- null_distribution %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat_two, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int)
null_distribution## Response: ccrdprs (numeric)
## Explanatory: ID (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 -0.125
## 2 2 -0.318
## 3 3 0.728
## 4 4 1.06
## 5 5 0.480
## 6 6 0.0666
## 7 7 -1.41
## 8 8 -0.290
## 9 9 -0.154
## 10 10 -1.83
## # ℹ 990 more rowsAnswer: it seems that the test stastic (red line) falls in the null distribution which means we fail to reject the null hypothesis. The red areas seem to be more extreme than the test statistic. This suggests that the mean extent one feels personal responsibility to reduce climate change is not much different for the those who are politically leaning left or right.