library(ggplot2)
library(dplyr)
library(statsr)load("gss.Rdata")Since 1972, the General Social Survey (GSS) has been monitoring societal change and studying the growing complexity of American society. The GSS aims to gather data on contemporary American society in order to monitor and explain trends and constants in attitudes, behaviors, and attributes; to examine the structure and functioning of society in general as well as the role played by relevant subgroups; to compare the United States to other societies in order to place American society in comparative perspective and develop cross-national models of human society; and to make high-quality data easily accessible to scholars, students, policy makers, and others, with minimal cost and waiting.
GSS questions include such items as national spending priorities, marijuana use, crime and punishment, race relations, quality of life, and confidence in institutions.
Since this survey is done by random sampling and the results can be “generalized” but not “causality” due to no random assignment.
As the survey concentrates in studying the complexity of American Society, we can aim to find out if the people’s confidence level in the society is based on their working status. If so, then government can work strongly towards improving them. It could solve questions like if the house keeping people are losing their confidence in the society.
Variables usedFor the confidence level, we will aggregate the 6 confidence variables confed, conbus, conedu, conpress, confinan, conmedic and conjudge.
Since it has 3 categorical values, we’ll assign a probability value as
| Category | value |
|---|---|
| A Great Deal | 1 |
| Only Some | 0.5 |
| Hardly Any | 0 |
We’ll then create new columns confidence_fed, confidence_bus, confidence_educ, confidence_press, confidence_finan, confidence_medic and confidence_judge by the below R code
gss <- gss %>% mutate(confidence_fed = ifelse(confed == "A Great Deal", 1,
ifelse(confed == "Only Some", 0.5, 0)));
gss <- gss %>% mutate(confidence_bus = ifelse(conbus == "A Great Deal", 1,
ifelse(conbus == "Only Some", 0.5, 0)));
gss <- gss %>% mutate(confidence_educ = ifelse(coneduc == "A Great Deal", 1,
ifelse(coneduc == "Only Some", 0.5, 0)));
gss <- gss %>% mutate(confidence_press = ifelse(conpress == "A Great Deal", 1,
ifelse(conpress == "Only Some", 0.5, 0)));
gss <- gss %>% mutate(confidence_finan = ifelse(confinan == "A Great Deal", 1,
ifelse(confinan == "Only Some", 0.5, 0)));
gss <- gss %>% mutate(confidence_medic = ifelse(conmedic == "A Great Deal", 1,
ifelse(conmedic == "Only Some", 0.5, 0)));
gss <- gss %>% mutate(confidence_judge = ifelse(conjudge == "A Great Deal", 1,
ifelse(conjudge == "Only Some", 0.5, 0)));Now that we have numerical columns, we will create an aggregate column confidence which adds up all the above new columns. So higher the value in confidence, the more the confidence in the society.
gss <- gss %>% mutate(confidence = confidence_fed+confidence_bus+confidence_educ+confidence_press+confidence_finan+confidence_medic+confidence_judge)With all the required columns, let’s compare the confidence in the box plot
ggplot(data=gss, aes(x=factor(wrkstat), y=confidence)) + geom_boxplot()## Warning: Removed 24784 rows containing non-finite values (stat_boxplot).
gss %>% filter(!is.na(confidence)) %>% group_by(wrkstat) %>% summarise(iqr = IQR(confidence), mean=mean(confidence), sd=sd(confidence), count=n()) %>% arrange(desc(iqr))## # A tibble: 9 x 5
## wrkstat iqr mean sd count
## <fct> <dbl> <dbl> <dbl> <int>
## 1 Unempl, Laid Off 2 3.62 1.38 1076
## 2 Retired 2 3.78 1.43 3967
## 3 School 2 4.00 1.34 1005
## 4 Other 2 3.47 1.48 602
## 5 Working Fulltime 1.5 3.77 1.30 16569
## 6 Working Parttime 1.5 3.82 1.32 3359
## 7 Temp Not Working 1.5 3.68 1.34 680
## 8 Keeping House 1.5 3.83 1.34 5013
## 9 <NA> 1.25 3.5 1.14 6
From the graph and the summary value, we can infer that all the working groups are nearly confidence as every one else with similar IQR and median. But is it the case? Check out the below inference.
With the mean/IQR values being some close to each other, let’s do the analysis of the data by ANOVA as we have the variable wrkstat with 8 different values.
Ho : nothing is going on. All the mean values are equal. So all working group people are in same confidence level
Ha : Confidence level is affected by the working status
1.1 Since the survey is from random sampling, we can assume independence
1.2 And the sample must be n < 10% the population
Check the distributions below
ggplot(data=gss %>% filter(wrkstat=="Working Fulltime", !is.na(confidence)), aes(x=confidence)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data=gss %>% filter(wrkstat=="Working Parttime", !is.na(wrkstat)), aes(x=confidence)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 2483 rows containing non-finite values (stat_bin).
ggplot(data=gss %>% filter(wrkstat=="Temp Not Working", !is.na(wrkstat)), aes(x=confidence)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 533 rows containing non-finite values (stat_bin).
ggplot(data=gss %>% filter(wrkstat=="Unempl, Laid Off", !is.na(wrkstat)), aes(x=confidence)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 797 rows containing non-finite values (stat_bin).
ggplot(data=gss %>% filter(wrkstat=="Retired", !is.na(wrkstat)), aes(x=confidence)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 3675 rows containing non-finite values (stat_bin).
ggplot(data=gss %>% filter(wrkstat=="School", !is.na(wrkstat)), aes(x=confidence)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 746 rows containing non-finite values (stat_bin).
ggplot(data=gss %>% filter(wrkstat=="Keeping House", !is.na(wrkstat)), aes(x=confidence)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 4374 rows containing non-finite values (stat_bin).
All the distributions are nearly normal within each group
var((gss %>% filter(wrkstat == "Working Fulltime" , !is.na(confidence)))$confidence)## [1] 1.691377
var((gss %>% filter(wrkstat == "Working Parttime" , !is.na(confidence)))$confidence)## [1] 1.735324
var((gss %>% filter(wrkstat == "Temp Not Working" , !is.na(confidence)))$confidence)## [1] 1.782723
var((gss %>% filter(wrkstat == "Unempl, Laid Off" , !is.na(confidence)))$confidence)## [1] 1.909654
var((gss %>% filter(wrkstat == "Retired" , !is.na(confidence)))$confidence)## [1] 2.056441
var((gss %>% filter(wrkstat == "School" , !is.na(confidence)))$confidence)## [1] 1.786105
var((gss %>% filter(wrkstat == "Keeping House" , !is.na(confidence)))$confidence)## [1] 1.807612
All the values are very close to each other.
gss %>% filter(!is.na(confidence)) %>% group_by(wrkstat) %>% summarise(mean=mean(confidence), count=n())## # A tibble: 9 x 3
## wrkstat mean count
## <fct> <dbl> <int>
## 1 Working Fulltime 3.77 16569
## 2 Working Parttime 3.82 3359
## 3 Temp Not Working 3.68 680
## 4 Unempl, Laid Off 3.62 1076
## 5 Retired 3.78 3967
## 6 School 4.00 1005
## 7 Keeping House 3.83 5013
## 8 Other 3.47 602
## 9 <NA> 3.5 6
We will do the inference for Anova
inference(y=confidence, x=wrkstat, data=gss %>% filter(!is.na(confidence)), type="ht", statistic = "mean", null=0, sig_level = 0.05, method="theoretical", alternative="greater")## Warning: Ignoring null value since it's undefined for ANOVA
## Response variable: numerical
## Explanatory variable: categorical (8 levels)
## n_Working Fulltime = 16569, y_bar_Working Fulltime = 3.7671, s_Working Fulltime = 1.3005
## n_Working Parttime = 3359, y_bar_Working Parttime = 3.8247, s_Working Parttime = 1.3173
## n_Temp Not Working = 680, y_bar_Temp Not Working = 3.6779, s_Temp Not Working = 1.3352
## n_Unempl, Laid Off = 1076, y_bar_Unempl, Laid Off = 3.6176, s_Unempl, Laid Off = 1.3819
## n_Retired = 3967, y_bar_Retired = 3.7779, s_Retired = 1.434
## n_School = 1005, y_bar_School = 3.9995, s_School = 1.3365
## n_Keeping House = 5013, y_bar_Keeping House = 3.8336, s_Keeping House = 1.3445
## n_Other = 602, y_bar_Other = 3.4693, s_Other = 1.4765
##
## ANOVA:
## df Sum_Sq Mean_Sq F p_value
## wrkstat 7 165.9201 23.7029 13.3153 < 0.0001
## Residuals 32263 57432.3299 1.7801
## Total 32270 57598.25
##
## Pairwise tests - t tests with pooled SD:
## # A tibble: 28 x 3
## group1 group2 p.value
## * <chr> <chr> <dbl>
## 1 Working Parttime Working Fulltime 0.0226
## 2 Temp Not Working Working Fulltime 0.0878
## 3 Unempl, Laid Off Working Fulltime 0.000369
## 4 Retired Working Fulltime 0.645
## 5 School Working Fulltime 0.0000000826
## 6 Keeping House Working Fulltime 0.00197
## 7 Other Working Fulltime 0.0000000752
## 8 Temp Not Working Working Parttime 0.00893
## 9 Unempl, Laid Off Working Parttime 0.00000942
## 10 Retired Working Parttime 0.135
## # ... with 18 more rows
With this higher F value, p-value is
pf(13.3153, 7, 32263, lower.tail = FALSE);## [1] 2.879938e-17
So reject the null hypothesis and conclude that confidence level in society is affected by the work status
As there is a relation between work status and confidence level, let’s check the Confidence Interval of each work status as it could be helpful in exactly understanding the mindset of the work group.
Let’s create some data which will contain only the data of each work status category as below
work_full_data <- gss %>% filter(!is.na(confidence), wrkstat == "Working Fulltime") %>% select(confidence)
work_full_result <- inference(y=confidence, data=work_full_data, type="ci", statistic = "mean", null=0, sig_level = 0.05, method="theoretical")## Single numerical variable
## n = 16569, y-bar = 3.7671, s = 1.3005
## 95% CI: (3.7473 , 3.7869)
work_part_data <- gss %>% filter(!is.na(confidence), wrkstat == "Working Parttime") %>% select(confidence)
work_part_data_result <- inference(y=confidence, data=work_part_data, type="ci", statistic = "mean", null=0, sig_level = 0.05, method="theoretical")## Single numerical variable
## n = 3359, y-bar = 3.8247, s = 1.3173
## 95% CI: (3.7801 , 3.8692)
temp_not_working_data <- gss %>% filter(!is.na(confidence), wrkstat == "Temp Not Working") %>% select(confidence)
temp_not_working_result <- inference(y=confidence, data=temp_not_working_data, type="ci", statistic = "mean", null=0, sig_level = 0.05, method="theoretical")## Single numerical variable
## n = 680, y-bar = 3.6779, s = 1.3352
## 95% CI: (3.5774 , 3.7785)
unemp_data <- gss %>% filter(!is.na(confidence), wrkstat == "Unempl, Laid Off") %>% select(confidence)
unemp_result <- inference(y=confidence, data=unemp_data, type="ci", statistic = "mean", null=0, sig_level = 0.05, method="theoretical")## Single numerical variable
## n = 1076, y-bar = 3.6176, s = 1.3819
## 95% CI: (3.5349 , 3.7002)
retired_data <- gss %>% filter(!is.na(confidence), wrkstat == "Retired") %>% select(confidence)
retired_result <- inference(y=confidence, data=retired_data, type="ci", statistic = "mean", null=0, sig_level = 0.05, method="theoretical")## Single numerical variable
## n = 3967, y-bar = 3.7779, s = 1.434
## 95% CI: (3.7333 , 3.8226)
school_data <- gss %>% filter(!is.na(confidence), wrkstat == "School") %>% select(confidence)
school_result <- inference(y=confidence, data=school_data, type="ci", statistic = "mean", null=0, sig_level = 0.05, method="theoretical")## Single numerical variable
## n = 1005, y-bar = 3.9995, s = 1.3365
## 95% CI: (3.9168 , 4.0822)
keeping_house_data <- gss %>% filter(!is.na(confidence), wrkstat == "Keeping House") %>% select(confidence)
keeping_house_result <- inference(y=confidence, data=keeping_house_data, type="ci", statistic = "mean", null=0, sig_level = 0.05, method="theoretical")## Single numerical variable
## n = 5013, y-bar = 3.8336, s = 1.3445
## 95% CI: (3.7964 , 3.8709)
And now we compare on where each work status category stand
ci_id <- c("Working Fulltime","Working Parttime", "Temp Not Working","Unempl, Laid Off", "Retired", "School", "Keeping House")
ci_bounds <- c(work_full_result$CI, work_part_data_result$CI, temp_not_working_result$CI, unemp_result$CI, retired_result$CI, school_result$CI, keeping_house_result$CI)
ci_data <- data.frame(ci_id,ci_bounds)
ggplot(data = ci_data, aes(x = ci_bounds, y = ci_id, group = ci_id)) + geom_point(size = 2) + geom_line() + geom_vline(xintercept = (gss %>% filter(!is.na(confidence)) %>% summarise(mean=mean(confidence)))$mean, color = "darkgray")