Muskal_Emily_Homework_3
Emily Muskal
2024-03-06
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.
packages <- c("tidyverse", "infer", "fst")
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_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## Warning: package 'infer' was built under R version 4.3.3## [[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"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.
packages <- c("tidyverse", "infer", "fst")france_data <- read.fst("france_data.fst")lapply(packages, library, character.only = TRUE)## [[1]]
## [1] "fstcore" "fst" "infer" "lubridate" "forcats" "stringr"
## [7] "dplyr" "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [13] "tidyverse" "stats" "graphics" "grDevices" "utils" "datasets"
## [19] "methods" "base"
##
## [[2]]
## [1] "fstcore" "fst" "infer" "lubridate" "forcats" "stringr"
## [7] "dplyr" "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [13] "tidyverse" "stats" "graphics" "grDevices" "utils" "datasets"
## [19] "methods" "base"
##
## [[3]]
## [1] "fstcore" "fst" "infer" "lubridate" "forcats" "stringr"
## [7] "dplyr" "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [13] "tidyverse" "stats" "graphics" "grDevices" "utils" "datasets"
## [19] "methods" "base"france_data <- read_fst("france_data.fst")france_data <- france_data %>%
filter(cntry == "FR") %>%
mutate(
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
)
) %>%
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)
)france_data <- france_data %>%
filter(cntry == "FR") %>%
mutate(
ID = 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_
),
Gen = case_when(
yrbrn %in% 1900:1945 ~ "1",
yrbrn %in% 1946:1964 ~ "2",
yrbrn %in% 1965:1979 ~ "3",
yrbrn %in% 1980:1996 ~ "4",
TRUE ~ as.character(yrbrn)
),
Gen = factor(Gen, levels = c("1", "2", "3", "4"),
labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
)table(france_data$lrscale)##
## 0 1 2 3 4 5 6 7 8 9 10 77 88
## 945 499 1269 1996 1838 5274 1594 1749 1458 409 686 381 940table(france_data$educ.ba)##
## BA or more No BA
## 4252 14786tab_dat <- france_data
mytab <- table(tab_dat$ID, tab_dat$educ.ba)
rsums <- rowSums(mytab)
csums <- colSums(mytab)
N <- sum(mytab)
ptab <- tcrossprod(rsums/N, csums/N)
cat("Table of Expected Proportions:\n")## Table of Expected Proportions:print(ptab)## [,1] [,2]
## [1,] 0.06135806 0.2044319
## [2,] 0.11343879 0.3779537
## [3,] 0.05605487 0.1867628ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2]
## [1,] 1087.0808 3621.919
## [2,] 2009.7951 6696.205
## [3,] 993.1241 3308.876alpha <- 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#The critical value for independence, set at 3.841, marks the point at which we would reject the null hypothesis, indicating a significant association between the variables lrscale and educ.ba. However, since our computed Chi-squared test statistic falls below this threshold, we fail to reject the null hypothesis. This outcome implies that we lack adequate evidence to assert an association between lrscale and educ.ba in our sample. Consequently, we cannot infer that one variable significantly influences the other based on our analysis.
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 <- read.fst("france_data.fst")
table(france_data$domicil)##
## 1 2 3 4 5 7 8
## 3584 2361 6257 5668 1162 3 3table(france_data$sbbsntx)##
## 1 2 3 4 5 7 8
## 757 1538 804 640 301 3 100france_clean <- france_data %>%
mutate(
domicil = ifelse(domicil == 2, 0, ifelse(domicil %in% c(7, 8, 9), NA, domicil)),
sbbsntx = ifelse(sbbsntx %in% c(77, 88, 99), NA, sbbsntx),) france_clean <- france_clean %>% filter(!is.na(domicil), !is.na(sbbsntx),) table(france_clean$domicil)##
## 0 1 3 4 5
## 487 699 1462 1249 245table(france_clean$sbbsntx)##
## 1 2 3 4 5 7 8
## 756 1538 804 640 301 3 100france_data <- france_data %>%
mutate(
sbbsntx = recode(as.character(sbbsntx),
'1' = "agree strongly",
'2' = "agree",
'3' = "neither agree or disagree",
'4' = "dissagree",
'5' = "strongly disagree",
'7' = NA_character_,
'8' = NA_character_,
'9' = NA_character_)
)%>%
filter(!is.na(sbbsntx))france_data <- france_data %>%
mutate(
domicil = 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(domicil))france_data <- france_data %>%
filter(!is.na(domicil) & !is.na(sbbsntx))
france_data <- france_data %>%
mutate(
domicil = factor(domicil),
sbbsntx = factor(sbbsntx)
) test_stat <- france_data %>%
specify(explanatory = domicil,
response = sbbsntx) %>%
hypothesize(null = "independence") %>%
calculate(stat = "chisq")
print(test_stat$stat)## X-squared
## 19.97187null_distribution <- france_data %>%
specify(explanatory = domicil,
response = sbbsntx) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "chisq")p_val1 <- 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_val1## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0#A p-value of 0.142 suggests that if the null hypothesis were true (i.e., if there were truly no relationship between the variables being tested), there would be a 14.2% chance of observing the data or more extreme results.
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.
null_distribution## Response: sbbsntx (factor)
## Explanatory: domicil (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 6.33
## 2 2 1.38
## 3 3 0.854
## 4 4 1.87
## 5 5 1.84
## 6 6 3.03
## 7 7 3.19
## 8 8 4.56
## 9 9 4.06
## 10 10 0.332
## # ℹ 990 more rowsthe data has a left skew and is also shifted to the left making it negative.
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 <- read.fst("france_data.fst")
table(france_data$lrscale)##
## 0 1 2 3 4 5 6 7 8 9 10 77 88
## 945 499 1269 1996 1838 5274 1594 1749 1458 409 686 381 940table(france_data$ccrdprs)##
## 0 1 2 3 4 5 6 7 8 9 10 66 77 88 99
## 72 26 64 78 114 478 403 736 952 358 708 27 1 29 1france_clean <- france_data %>%
mutate(
lrscale = ifelse(lrscale == 2, 0, ifelse(lrscale %in% c(7, 8, 9), NA, lrscale)),
ccrdprs = ifelse(ccrdprs %in% c(77, 88, 99), NA, ccrdprs),) france_clean <- france_clean %>% filter(!is.na(lrscale), !is.na(ccrdprs),) table(france_clean$lrscale)##
## 0 1 3 4 5 6 10 77 88
## 395 87 412 340 1098 327 156 110 257table(france_clean$ccrdprs)##
## 0 1 2 3 4 5 6 7 8 9 10 66
## 55 21 50 58 85 378 318 554 743 295 604 21france_data <- france_data %>%
mutate(
ccrdprs = recode(as.character(ccrdprs),
'1' = "not at all",
'2' = "agree",
'3' = "neither agree or disagree",
'4' = "dissagree",
'5' = "strongly disagree",
'7' = NA_character_,
'8' = NA_character_,
'9' = NA_character_)
)%>%
filter(!is.na(ccrdprs))france_data <- france_data %>%
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% 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'
)%>%
filter(!is.na(lrscale))france_data <- france_data %>%
filter(!is.na(lrscale) & !is.na(ccrdprs))
france_data <- france_data %>%
mutate(
lrscale = factor(lrscale),
ccrdprs = factor(ccrdprs)
) test_stat <- france_data %>%
specify(explanatory = lrscale,
response = ccrdprs) %>%
hypothesize(null = "independence") %>%
calculate(stat = "chisq")
print(test_stat$stat)## X-squared
## 215.049null_distribution <- france_data %>%
specify(explanatory = lrscale,
response = ccrdprs) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "chisq")p_val1 <- 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_val1## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0#A p-value of 0.142 suggests that if the null hypothesis were true (i.e., if there were truly no relationship between the variables being tested), there would be a 14.2% chance of observing the data or more extreme results.
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.
null_distribution## Response: ccrdprs (factor)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 141.
## 2 2 162.
## 3 3 155.
## 4 4 127.
## 5 5 133.
## 6 6 141.
## 7 7 114.
## 8 8 156.
## 9 9 149.
## 10 10 134.
## # ℹ 990 more rows