Loading packages
# 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"
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.
Copying dataframe and cleaning ‘lrscale’ and ‘educ.ba’
variables
france_data_lr <- france_data |>
mutate(
lrscale = case_when(
lrscale %in% 0:3 ~ "Left",
lrscale %in% 7:10 ~ "Right",
TRUE ~ NA_character_ # Moderate (4, 5, 6) and special codes (77, 88, 99) set to NA
),
# Re-coding education with case_when
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
) ) |>
# Cleaning NA's
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
)
Double-checking that variables are clean
unique(france_data_lr$lrscale)
## [1] "Left" NA "Right"
unique(france_data_lr$educ.ba)
## [1] "No BA" "BA or more"
Removing ‘lrscale’ NA values
france_data_lr <- france_data_lr |> filter(!is.na(lrscale))
unique(france_data_lr$lrscale)
## [1] "Left" "Right"
table(france_data_lr$lrscale)
##
## Left Right
## 4709 4302
Long Form Coding for Chi-square
Calculating and Interpreting Critical value for Independence
tab_dat <- france_data_lr # making copy of data
# Creating a crosstable and save it into the object mytab
mytab <- table(tab_dat$lrscale, tab_dat$educ.ba) # tabulates two columns of left or right and whether have BA
# Calculate row-wise sums of counts in table
rsums <- rowSums(mytab) # row sums of table - left or right
# Calculate column-wise sums of counts in table
csums <- colSums(mytab) # column sums - having bachelors vs. not
# Get the total count of the table (N)
N <- sum(mytab)
# Generate table of expected proportions - take row sum and column sum and create proportions if expected no relationship between variables
ptab <- tcrossprod(rsums/N, csums/N) # tcrossprod calculates expected proportions
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) #qchisq = chi-square test - gives critical value of alpha with degrees of freedom as 1 - only two possible options you can choose for variables which is why degrees of freedom is 1
cat("Critical Value:", round(c_val, 3), "\n") # cat = print this text "critical value" and round it to 3 decimal points
## Critical Value: 3.841
# \n = print on new line
Interpreting Critical Value for Independence
The critical value for ‘lrscale’ and ‘educ.ba’ variables is 3.841.
The critical value denotes the rejection region or threshold of
rejecting the null hypothesis depending on which statistical test is
being used.
A value of 3.841 means that if our observed statistic value is above
3.841, it would be within the rejection region, meaning that political
orientation and education level are not independent at the 0.05
significance level. As a result, the null hypothesis will be
rejected.
Calculating and interpreting Pearson’s chi-squared statistic
# Pearson's chi-squared statistic:
test_stat_lr <- sum((mytab - ftab)^2 / ftab) # calculating observed chi-square statistic (sum of) observed - expected squared over expected
cat("Pearson's X^2:", round(test_stat_lr, 4), "\n")
## Pearson's X^2: 100.8551
Interpreting Pearson’s chi-squared statistic
Comparison to the critical value: The observed test statistic (x^2)
of 100.8551 is much larger than our critical value of 3.841. This means
that our observed test statistic is within the rejection region, so the
null hypothesis of no difference between political orientation and
education level is rejected. The difference between our observed data
and what we’d expect under the assumption of independence is
statistically significant.
Thus, the probability of being left or right for political
orientation is not the same across different education levels.
Calculating the p-value
# p-value:
p_val_lr <- pchisq(test_stat_lr, df = 1, lower.tail = F) # p instead of q in chisq - saying get me p value instead
cat("p-value:", round(p_val_lr, 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)) # running chi-square test
##
## Pearson's Chi-squared test
##
## data: mytab
## X-squared = 100.86, df = 1, p-value < 2.2e-16
Interpreting the p-value
p-value < 2.2e-16, which means that is is extremely close to 0.
Thus, we can conclude there is strong evidence to reject the null
hypothesis that political orientation and education level are
independent from each other because the p-value shows an almost 0%
chance that the observed statistic will occur under the null
distribution.
The large chi-squared statistic of 100.86 further reinforces this
point. If the two variables were independent, we’d expect a chi-squared
statistic closer to 0.
Main Takeaways
The chi-square test determines if two categorical variables are
independent of each other. According to our observed statistic of 100.86
in comparison to our critical value of 3.841, our test statistic
suggests that the relationship between political orientation and
education level is not statistically independent, as it exceeds the
critical value. The p-value reinforces this conclusion because it is
almost 0, which is less than our 0.05 alpha level of significance. Thus,
whether an individual is right or left on the political orientation
scale is dependent or explained by whether they have BA or not.
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).
Copying dataframe and re-coding ‘sbbsntx’ and ‘domicil’
france_data_sb <- france_data |>
mutate(
sbbsntx = case_when(
sbbsntx == 1 ~ "Agree strongly",
sbbsntx == 2 ~ "Agree",
sbbsntx == 3 ~ "Neither agree nor disagree",
sbbsntx == 4 ~ "Disagree",
sbbsntx == 5 ~ "Disagree strongly",
TRUE ~ NA_character_
),
sbbsntx = fct_relevel(factor(sbbsntx),
"Agree strongly",
"Agree",
"Neither agree nor disagree",
"Disagree",
"Disagree strongly")
) |>
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))
Double-checking if re-coding is successful
table(france_data_sb$sbbsntx)
##
## Agree strongly Agree
## 756 1538
## Neither agree nor disagree Disagree
## 804 640
## Disagree strongly
## 301
table(france_data_sb$domicil)
##
## Rural Urban
## 13087 5945
unique(france_data_sb$domicil)
## [1] "Rural" "Urban"
unique(france_data_sb$sbbsntx)
## [1] <NA> Agree strongly
## [3] Neither agree nor disagree Agree
## [5] Disagree strongly Disagree
## 5 Levels: Agree strongly Agree Neither agree nor disagree ... Disagree strongly
Cleaning ‘sbbsntx’ NA values again
france_data_sb <- france_data_sb |> filter(!is.na(sbbsntx))
unique(france_data_sb$sbbsntx)
## [1] Agree strongly Neither agree nor disagree
## [3] Agree Disagree strongly
## [5] Disagree
## 5 Levels: Agree strongly Agree Neither agree nor disagree ... Disagree strongly
Step 1: Calculate the test statistic of your sample
test_stat_sb <- france_data_sb |>
specify(explanatory = domicil, # change variable name for explanatory variable
response = sbbsntx) |> # change variable name for outcome of interest
hypothesize(null = "independence") |> # store data until we do a test
calculate(stat = "Chisq") # replace in between quotation marks appropriate test statistic
print(test_stat_sb$stat)
## X-squared
## 19.97187
Interpreting Value
19.97187 is the observed test statistic (X^2) that we will compare to
the null distribution to determine the significance between geographic
location and believing that social benefits/services cost businesses too
much in taxes. This is because we don’t have a critical value to compare
the observed test statistic to with the infer package.
Step 2: Simulate the null distribution
null_distribution_sb <- france_data_sb |>
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 = "Chisq") # 1000 different times it calculates t statistic
Interpreting Null Simulation
This is creating a simulation of the distribution of expected values
of ‘domicil’ and ‘sbbsntx’ when they are cross-tabulated, under the
hypothetical assumption that they are independent. This means assuming
there is no difference in people’s geographic location and whether they
feel that social benefits/services cost businesses too much in
taxes.
Step 3: Calculate the p-value of your sample
p_val_sb <- null_distribution_sb |> # replace name here if you assigned something other than null_distribution above
get_pvalue(obs_stat = test_stat_sb, direction = "greater") # 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.
# given that 13 = calculated observed value how many times does the expected test statistics exceed the critical value
p_val_sb
## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0
Interpreting P-Value
The p-value is very close to 0, which indicates an almost 0% chance
that the observed test statistic of 19.97 will be captured in the null
distribution, suggesting that there is a significant association between
a person’s geographic location and whether they believe social
benefits/services cost businesses too much in taxes.
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.
Step 4: Visualize without confidence intervals
null_distribution_sb |>
visualize() +
shade_p_value(obs_stat = test_stat_sb, direction = "greater") # add on observed statistic and direction - two sided test
## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
null_distribution_sb
## Response: sbbsntx (factor)
## Explanatory: domicil (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 1.52
## 2 2 0.0557
## 3 3 5.67
## 4 4 6.91
## 5 5 1.44
## 6 6 2.46
## 7 7 6.24
## 8 8 2.81
## 9 9 7.55
## 10 10 1.62
## # ℹ 990 more rows
Interpretation
According to the null distribution, the observed test statistic of
19.97 falls outside of the null distribution. This suggests there is
significant evidence to reject the null hypothesis under assumptions of
frequentist hypothesis testing because if the null hypothesis is true,
our observed test result is significantly different from that
expectation.
Thus, it suggests that whether a person believes that social
benefits/services cost businesses too much in taxes is dependent or can
be explained by the person’s geographic location.
Step 4: Visualize with confidence intervals
conf_int_sb <- null_distribution_sb |>
get_confidence_interval(level = 0.95, type = "percentile") # showing confidence interval - if red line was in green boundary - cannot reject null
null_distribution_sb |>
visualize() +
shade_p_value(obs_stat = test_stat_sb, direction = "greater") +
shade_confidence_interval(endpoints = conf_int_sb)
## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
null_distribution_sb
## Response: sbbsntx (factor)
## Explanatory: domicil (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 1.52
## 2 2 0.0557
## 3 3 5.67
## 4 4 6.91
## 5 5 1.44
## 6 6 2.46
## 7 7 6.24
## 8 8 2.81
## 9 9 7.55
## 10 10 1.62
## # ℹ 990 more rows
Interpreting Output
The confidence interval showcases the range of not rejecting the null
in the simulated null distribution. Thus, anything outside of the two
boundaries of the confidence interval indicates the rejection region of
the chi-square test. Since our observed test statistic (X^2) of 19.97 is
outside of the confidence interval range, it means it is within the
rejection region, so the null hypothesis is rejected.
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).
Re-coding ‘lrscale’ variable
france_data_ccr <- france_data |>
mutate(
lrscale = case_when(
lrscale %in% 0:3 ~ "Left", # Left-wing (0 to 3)
lrscale %in% 7:10 ~ "Right", # Right-wing (7 to 10)
TRUE ~ NA_character_ # Moderate (4, 5, 6) and special codes (77, 88, 99) set to NA
))
Checking ‘lrscale’ and cleaning it
unique(france_data_ccr$lrscale)
## [1] "Left" NA "Right"
france_data_ccr <- france_data_ccr |> filter(!is.na(lrscale))
Checking ‘ccrdprs’ variable and cleaning it
table(france_data_ccr$ccrdprs)
##
## 0 1 2 3 4 5 6 7 8 9 10 66 88
## 33 12 28 38 57 205 173 344 472 174 335 13 13
france_data_ccr$ccrdprs[france_data_ccr$ccrdprs %in% 66:88] <- NA
table(france_data_ccr$ccrdprs)
##
## 0 1 2 3 4 5 6 7 8 9 10
## 33 12 28 38 57 205 173 344 472 174 335
unique(france_data_ccr$ccrdprs)
## [1] NA 10 8 7 6 5 9 0 4 2 3 1
france_data_ccr <- france_data_ccr |> filter(!is.na(ccrdprs))
Double-checking ‘lrscale’ and ‘ccrdprs’ variables
unique(france_data_ccr$lrscale)
## [1] "Left" "Right"
unique(france_data_ccr$ccrdprs)
## [1] 10 8 7 6 5 9 0 4 2 3 1
Step 1: Calculate the test statistic of your sample
test_stat_ccr <- france_data_ccr |>
specify(explanatory = lrscale, # change variable name for explanatory variable
response = ccrdprs) |> # change variable name for outcome of interest
hypothesize(null = "independence") |> # store data until we do a test
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_ccr$stat)
## t
## 5.482311
Interpreting Value
5.482311 is the observed test statistic that we will compare to the
null distribution to determine the significance between the extent
people feel a personal responsibility to reduce climate change and their
political orientation. This is because we don’t have a critical value to
compare the observed test statistic to with the infer package.
Step 2: Simulate the null distribution
null_distribution_ccr <- france_data_ccr |>
specify(explanatory = lrscale,
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") # 1000 different times it calculates t 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.
Interpreting Null Simulation
This is creating a simulation of the distribution of expected values
of ‘ccrdprs’ and ‘lrscale’ when they are cross-tabulated, under the
hypothetical assumption that they are independent. This means assuming
there is no difference in the extent that people feel a personal
responsibility to reduce climate change and their political
orientation.
Step 3: Calculate the p-value of your sample
p_val_ccr <- null_distribution_ccr |> # replace name here if you assigned something other than null_distribution above
get_pvalue(obs_stat = test_stat_ccr, 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.
# given that 13 = calculated observed value how many times does the expected test statistics exceed the critical value
p_val_ccr
## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0
Interpreting P-Value
The p-value is very close to 0, which indicates an almost 0% chance
that the observed test statistic of 5.48 will be captured in the null
distribution, suggesting that there is a significant association between
the extent that people feel a personal responsibility to reduce climate
change and their political orientation.
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 without confidence intervals
null_distribution_ccr |>
visualize() +
shade_p_value(obs_stat = test_stat_ccr, direction = "two-sided") # add on observed statistic and direction - two sided test
null_distribution_ccr
## Response: ccrdprs (numeric)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 -0.512
## 2 2 -0.00361
## 3 3 0.0176
## 4 4 0.123
## 5 5 -0.215
## 6 6 1.78
## 7 7 1.18
## 8 8 0.123
## 9 9 0.803
## 10 10 0.887
## # ℹ 990 more rows
Interpretation
According to the null distribution, the observed test statistic of
5.48 falls outside of the null distribution. This suggests there is
significant evidence to reject the null hypothesis under assumptions of
frequentist hypothesis testing because if the null hypothesis is true,
our observed test result is significantly different from that
expectation.
Thus, it suggests that the extent a person feels personal
responsibility to reduce climate change is dependent or can be explained
by the person’s political orientation.
Step 4: Visualize with confidence intervals
conf_int_ccr <- null_distribution_ccr |>
get_confidence_interval(level = 0.95, type = "percentile") # showing confidence interval - if red line was in green boundary - cannot reject null
null_distribution_ccr |>
visualize() +
shade_p_value(obs_stat = test_stat_ccr, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int_ccr)
null_distribution_ccr
## Response: ccrdprs (numeric)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 -0.512
## 2 2 -0.00361
## 3 3 0.0176
## 4 4 0.123
## 5 5 -0.215
## 6 6 1.78
## 7 7 1.18
## 8 8 0.123
## 9 9 0.803
## 10 10 0.887
## # ℹ 990 more rows
Interpreting Output
The confidence interval showcases the range of not rejecting the null
in the simulated null distribution. Thus, anything outside of the two
boundaries of the confidence interval indicates the rejection region of
our t-test. Since our observed test statistic of 5.48 is outside of the
confidence interval range, it means it is within the rejection region,
so the null hypothesis is rejected. The extent that a person feels
personal responsibility to reduce climate change can be explained by
their political orientation.