Exploratory Data Analysis

yrbss_weight <- yrbss %>%
  select(weight) #selecting the weight column and allocating it 

yrbss_weight %>%
  summary #checking the descriptive statistics of the weight

##      weight    
##  Min.   : 30   
##  1st Qu.: 56   
##  Median : 64   
##  Mean   : 68   
##  3rd Qu.: 76   
##  Max.   :181   
##  NA's   :1004

yrbss_weight %>%
  ggplot(aes(x=weight)) +
  geom_histogram(binwidth=5) +
  labs(title="Histogram of weights", x="Weight", y="Count") +
  NULL #plotting a histogram of the weights

yrbss_weight %>%
  filter(is.na(weight)) %>%
  count() %>%
  paste('missing values') #tCounting the missing values in weight i.e. 1004 missing values

## [1] "1004 missing values"

Summarising physical exercise data

yrbss <- yrbss %>%
  mutate(physical_3plus=if_else(physically_active_7d >= 3, 'yes', 'no'))
#making a new column titles physical_3plus if the person is active for more than 3 days a week

yrbss %>%
  filter(!is.na(physical_3plus)) %>% #removing the N.A.'s
  group_by(physical_3plus) %>%
  summarise(n=n()) %>%
  mutate(prop=round(n/sum(n), 3)) %>% #creating a new column for the proportion of physical_3plus
  arrange(-prop) #arranging the proportion in descending order 

## # A tibble: 2 x 3
##   physical_3plus     n  prop
##   <chr>          <int> <dbl>
## 1 yes             8906 0.669
## 2 no              4404 0.331

yrbss %>%
  filter(!is.na(physical_3plus)) %>% #removing the N.A.'s in physical_3plus
  count(physical_3plus, sort=T) %>%
  mutate(prop=round(n/sum(n), 3)) #getting the number and proportion of people active for more than 3 days

## # A tibble: 2 x 3
##   physical_3plus     n  prop
##   <chr>          <int> <dbl>
## 1 yes             8906 0.669
## 2 no              4404 0.331

Calculating a confidence interval for the population proportion of high schools that are NOT active 3 or more days per week

yrbss %>%
  filter(!is.na(physical_3plus)) %>% #Removing the N.A.'s 
  count(physical_3plus, sort=T) %>%
  mutate(prop=n/sum(n)) %>% #Finding the proportion 
  filter(physical_3plus == 'no') %>%
  summarise(
    count=sum(n),
    se=sqrt(prop*(1-prop)/count), #getting the standard error (0.00709)
    t_critical=qt(0.975, count-1), #getting the t-critical value i.e. 1.96 (at 95% level)
    lower=prop-t_critical*se,
    upper=prop+t_critical*se) #getting the confidence interval (upper and lower CI)

## # A tibble: 1 x 5
##   count      se t_critical lower upper
##   <int>   <dbl>      <dbl> <dbl> <dbl>
## 1  4404 0.00709       1.96 0.317 0.345

Boxplot of physical_3plus vs. weight

yrbss %>%
  filter(!is.na(physical_3plus)) %>% #Removing the NA's  
  ggplot(aes(x=physical_3plus, y=weight)) +
  geom_boxplot() + #box plot of weights for people that are active for more or less than 3 days
  labs(title="Relationship between exercise frequency and weight", x="Exercise at least 3 days per week?", y="Weight") +
  NULL

Inference 1 Although, the weights of the people that exercise more or less than 3 days appear to be familiar, there is a relationship between the level of physical activity and respective weight. Weight appears to be more clustered for the people that exercise over 3 days and there are evidently more outliers (higher weight) for the people exercising less than 3 days. Therefore, we infer that exercising is helpful in maintaining weight. Nonetheless, it is important to keep in mind that there are some people that do weight training to put on muscle (in tuen increasing weight), weights may not be the best indicator of fitness.

Calculating the confidence Interval

yrbss %>% 
  filter(!is.na(physical_3plus)) %>% #Removing the NA's 
  group_by(physical_3plus) %>% #grouping by the people that exercise for more than 3 days
  summarise(
    mean=mean(weight, na.rm=T), #finding the descriptive statistics 
    sd=sd(weight, na.rm=T),
    count=n(),
    se=sd/sqrt(count), #finding the standard error, t critical value, and confidence interval
    t_critical=qt(0.975, count-1),
    lower=mean-t_critical*se,
    upper=mean+t_critical*se)

## # A tibble: 2 x 8
##   physical_3plus  mean    sd count    se t_critical lower upper
##   <chr>          <dbl> <dbl> <int> <dbl>      <dbl> <dbl> <dbl>
## 1 no              66.7  17.6  4404 0.266       1.96  66.2  67.2
## 2 yes             68.4  16.5  8906 0.175       1.96  68.1  68.8

There is an observed difference of about 1.77kg (68.44 - 66.67), and we notice that the two confidence intervals do not overlap. It seems that the difference is at least 95% statistically significant. Let us also conduct a hypothesis test.

Hypothesis testing

Null: There is no difference in weight means between people who are active more than or less than 3 days per week. Alternative: There is a difference in weight means between people who are active more than or less than 3 days per week.

t.test(weight ~ physical_3plus, data = yrbss) 

## 
##  Welch Two Sample t-test
## 
## data:  weight by physical_3plus
## t = -5, df = 7479, p-value = 9e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.42 -1.12
## sample estimates:
##  mean in group no mean in group yes 
##              66.7              68.4

#Running a t-test to check for the relation between weight and physical activity

obs_diff <- yrbss %>% 
  specify(weight ~ physical_3plus) %>%
  calculate(stat = "diff in means", order = c("yes", "no")) #initialising the test

null_dist <- yrbss %>%
  # specify variables
  specify(weight ~ physical_3plus) %>%
  
  # assume independence, i.e, there is no difference
  hypothesize(null = "independence") %>%
  
  # generate 1000 reps, of type "permute"
  generate(reps = 1000, type = "permute") %>%
  
  # calculate statistic of difference, namely "diff in means"
  calculate(stat = "diff in means", order = c("yes", "no"))

ggplot(data = null_dist, aes(x = stat)) +
  geom_histogram() + #Histogram for the null distribution 
  labs(title="Histogram for null distribution", x="Difference in means", y="Count") +
  NULL

null_dist %>% visualize() +
  #visualizing and shading for the p-value of our hypothesis
  shade_p_value(obs_stat = obs_diff, direction = "two-sided") 

null_dist %>%
  #calculating the p-value for our hypothesis
  get_p_value(obs_stat = obs_diff, direction = "two_sided")

## # A tibble: 1 x 1
##   p_value
##     <dbl>
## 1       0

Comparing directors by ratings of their movies

directors <- c('Tim Burton', 'Steven Spielberg')

movies %>%
  filter(director %in% directors) %>% #only keeping Tim Burton and Steven Spielberg
  select(director, rating) %>% 
  group_by(director) %>% 
  #selecting only the director and ratings then further grouping by director
  summarise(mean_rating = mean(rating),
            median_rating = median(rating),
            sd_rating = sd(rating),
            count = n(),
            t_critical = qt(0.975, count-1),
            se_rating = sd_rating/sqrt(count),
            margin_of_error = t_critical * se_rating,
            rating_low = mean_rating - margin_of_error,
            rating_high = mean_rating + margin_of_error) %>% 
  #getting the descriptive statistics, standard error, t-critical value, margin of error, and confidence interval
  mutate(
    xmin=max(rating_low),
    xmax=min(rating_high)) %>%
  #creating a plot of the mean rating of various directors and differentiating by the colors
  ggplot(aes(x=mean_rating, y=factor(director, level=directors), colour=director)) +
  #creating a scatter plot 
  geom_point() +
  geom_errorbarh(aes(xmin=rating_low, xmax=rating_high), width=0.1, size=2) +
  geom_rect(aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf, alpha=0.05)) +
  #adding text on the rectangle
  annotate('text', x=6.53, y=directors[1], label="6.53") +
  annotate('text', x=6.93, y=directors[1], label="6.93") +
  annotate('text', x=7.33, y=directors[1], label="7.33") +
  annotate('text', x=7.27, y=directors[2], label="7.27") +
  annotate('text', x=7.57, y=directors[2], label="7.57") +
  annotate('text', x=7.87, y=directors[2], label="7.87") +
  #adding title, subtitle, and the description of x and y axis
  labs(x="Mean IMDB Rating",
       y="", 
       title="Do Speilberg and Burton have the same mean IMDB ratings?",
      subtitle="95% confidence intervals overlap") +
  theme_bw() +
  theme(legend.position='none') +
  #removing the legend 
  NULL

Loading the data

omega <- read_csv(here::here("data", "omega.csv"))
glimpse(omega) #Exploring the omega dataset

## Rows: 50
## Columns: 3
## $ salary     <dbl> 81894, 69517, 68589, 74881, 65598, 76840, 78800, 70033, 635…
## $ gender     <chr> "male", "male", "male", "male", "male", "male", "male", "ma…
## $ experience <dbl> 16, 25, 15, 33, 16, 19, 32, 34, 1, 44, 7, 14, 33, 19, 24, 3…

Relationship Salary - Gender ?

# Summary Statistics of salary by gender
mosaic::favstats(salary ~ gender, data=omega)

##   gender   min    Q1 median    Q3   max  mean   sd  n missing
## 1 female 47033 60338  64618 70033 78800 64543 7567 26       0
## 2   male 54768 68331  74675 78568 84576 73239 7463 24       0

# Dataframe with two rows (male-female) and having as columns gender, mean, SD, sample size, 
# the t-critical value, the standard error, the margin of error, 
# and the low/high endpoints of a 95% condifence interval

omega %>%
  group_by(gender) %>%
  summarise(
    mean=mean(salary),
    sd=sd(salary),
    count=n(),
    se=sd/sqrt(count),
    t_critical=qt(0.975, count-1),
    lower=mean-t_critical*se,
    upper=mean+t_critical*se) %>%
  as.data.frame()

##   gender  mean   sd count   se t_critical lower upper
## 1 female 64543 7567    26 1484       2.06 61486 67599
## 2   male 73239 7463    24 1523       2.07 70088 76390

Inference 2 Evidently, the mean salary for females is lower than that of males, but the variation for both the gender’s is similar i.e. there is not a huge gap in the salary levels within each gender. Additionally, since the t-critical value is greater than 2, with 95% level of confidence, we can reject our null hypothesis that mean salary of males subtracted from the mean salary of females is zero (it is also worthwhile to notice that the confidence intervals do not overlap). You can also run a hypothesis testing, assuming as a null hypothesis that the mean difference in salaries is zero, or that, on average, men and women make the same amount of money. You should tun your hypothesis testing using t.test() and with the simulation method from the infer package.

Hypothesis testing for salary difference between genders

# hypothesis testing using t.test() 
t.test(salary ~ gender, data=omega)

## 
##  Welch Two Sample t-test
## 
## data:  salary by gender
## t = -4, df = 48, p-value = 2e-04
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -12973  -4420
## sample estimates:
## mean in group female   mean in group male 
##                64543                73239

# hypothesis testing using infer package
obs_diff <- omega %>%
  specify(response=salary) %>%
  calculate(stat='mean')

omega %>%
  specify(salary ~ gender) %>%
  hypothesize(null='independence') %>%
  generate(reps=1000, type='permute') %>%
  calculate(stat='diff in means', order=c('female', 'male')) %>%
  get_p_value(obs_stat=obs_diff, direction='both')

## # A tibble: 1 x 1
##   p_value
##     <dbl>
## 1       0

Inference 3: The p-value from the t-test is 0.0002 and from get_p_value is 0 (which is much lesser than the hurdle of 0.05). Therefore with 95% level of confidence, we can reject the null hypothesis that the mean difference in salaries is 0. Further, since this p-value seems improbable, it may be worthwhile to consider if there are floating-point numbers or there has been omission of important data points.

Relationship Experience - Gender?

# Summary Statistics of salary by gender
favstats(experience ~ gender, data=omega)

##   gender min    Q1 median   Q3 max  mean    sd  n missing
## 1 female   0  0.25    3.0 14.0  29  7.38  8.51 26       0
## 2   male   1 15.75   19.5 31.2  44 21.12 10.92 24       0

t.test(experience ~ gender, data=omega)

## 
##  Welch Two Sample t-test
## 
## data:  experience by gender
## t = -5, df = 43, p-value = 1e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -19.35  -8.13
## sample estimates:
## mean in group female   mean in group male 
##                 7.38                21.12

Inference 4: This p-value below 0.05 suggests that at 95% confidence level we can reject the null hypothesis that the mean difference between the work experience of males and females is 0. This finding definately endangers our earlier assumption that the difference in mean salaries of males and females is due to some inherent bias (gender disparity at workplace). However, after this analysis we can say that one of the reasons for lower salaries of females is the lesser work experience as comapred to that of men (a case of correlation not causation).

Relationship Salary - Experience ?

#making a scatter plot of the experience and salary of employees
omega %>%
  ggplot(aes(x=experience, y=salary)) +
  geom_point() +
  geom_smooth(se=F) +
  NULL

Checking the correlations between the data

#creating a correlation matrix to see the relationship between experience and salary seperately for men and women
omega %>%
  select(gender, experience, salary) %>%
  ggpairs(aes(colour=gender, alpha = 0.3)) +
  theme_bw()

Inference 5: On a closer look, some of the possible inferencs are: 1. The feamles appear to be at the lower level of jobs (interns, junior managers, etc) whereas the top managament is clustered with male. 2. Concievably, the correlation between salary and experience is more for the females than males i.e. despite lesser work experience it is possible for males to have higher paying jobs.