Load the data

This data is part of the openintro textbook and we can load and inspect it. There are observations on 13 different variables, some categorical and some numerical. The meaning of each variable can be found by bringing up the help file:

?yrbss

data(yrbss)
glimpse(yrbss)

## Rows: 13,583
## Columns: 13
## $ age                      <int> 14, 14, 15, 15, 15, 15, 15, 14, 15, 15, 15, 1~
## $ gender                   <chr> "female", "female", "female", "female", "fema~
## $ grade                    <chr> "9", "9", "9", "9", "9", "9", "9", "9", "9", ~
## $ hispanic                 <chr> "not", "not", "hispanic", "not", "not", "not"~
## $ race                     <chr> "Black or African American", "Black or Africa~
## $ height                   <dbl> NA, NA, 1.73, 1.60, 1.50, 1.57, 1.65, 1.88, 1~
## $ weight                   <dbl> NA, NA, 84.4, 55.8, 46.7, 67.1, 131.5, 71.2, ~
## $ helmet_12m               <chr> "never", "never", "never", "never", "did not ~
## $ text_while_driving_30d   <chr> "0", NA, "30", "0", "did not drive", "did not~
## $ physically_active_7d     <int> 4, 2, 7, 0, 2, 1, 4, 4, 5, 0, 0, 0, 4, 7, 7, ~
## $ hours_tv_per_school_day  <chr> "5+", "5+", "5+", "2", "3", "5+", "5+", "5+",~
## $ strength_training_7d     <int> 0, 0, 0, 0, 1, 0, 2, 0, 3, 0, 3, 0, 0, 7, 7, ~
## $ school_night_hours_sleep <chr> "8", "6", "<5", "6", "9", "8", "9", "6", "<5"~

Before we carry on with our analysis, it’s is always a good idea to check with skimr::skim() to get a feel for missing values, summary statistics of numerical variables, and a very rough histogram.

skim(yrbss)

Variable type: character

Variable type: numeric

Exploratory Data Analysis

We first start with analysing the weight of participants in kilograms. Using visualization and summary statistics, we can describe the distribution of weights:

The distribution of weight is right(positive) skewed. 1004 observations are missed in the weight from.

skim(yrbss)

Variable type: character

Variable type: numeric

ggplot(data = yrbss)+
  geom_density(aes(x = weight), kernel = "gaussian") 

Next, consider the possible relationship between a high schooler’s weight and their physical activity. Plotting the data is a useful first step because it helps us quickly visualize trends, identify strong associations, and develop research questions.

Let’s create a new variable in the dataframe yrbss, called physical_3plus , which will be yes if they are physically active for at least 3 days a week, and no otherwise.

yrbss <- yrbss %>% 
  # count(physically_active_7d) %>% 
  mutate(yrbss,physical_3plus = ifelse(physically_active_7d >= 3,"yes","no"))


yrbss_count <- yrbss %>% 
  count(physical_3plus) %>% 
  mutate(perc = n/sum(n))
yrbss_count

## # A tibble: 3 x 3
##   physical_3plus     n   perc
##   <chr>          <int>  <dbl>
## 1 no              4404 0.324 
## 2 yes             8906 0.656 
## 3 <NA>             273 0.0201

There is a slightly positive relationship between thse two variables. I expected that there would be positive relationship between physical_3plus and weight because if people did more execrises per week, they would have more muscles which would be much heavier that fat.

ybrss_physical3minus <- 
  filter(yrbss, physical_3plus == "no")
# ybrss_physical3minus
# yrbss
ggplot(data = yrbss,aes(x = physical_3plus, y =weight))+
  geom_boxplot()

Confidence Interval

Boxplots show how the medians of the two distributions compare, but we can also compare the means of the distributions using either a confidence interval or a hypothesis test. Note that when we calculate the mean, SD, etc. weight in these groups using the mean function, we must ignore any missing values by setting the na.rm = TRUE.

# t.test(weight ~ physical_3plus, data = yrbss)
yrbss %>% 
  group_by(physical_3plus) %>% 
  summarise(mean_weight = mean(weight,na.rm = TRUE))

## # A tibble: 3 x 2
##   physical_3plus mean_weight
##   <chr>                <dbl>
## 1 no                    66.7
## 2 yes                   68.4
## 3 <NA>                  69.9

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 test with formula

The null hypotheses : The mean weights of people who exercise at least times a week is the same as the mean weights of people who don’t. The alternative hypothesesz: The mean weights of people who exercise at least times a week is significantly different from the mean weights of people who don’t.

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 between group no and group yes 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

Hypothesis test with infer

Next, we will introduce a new function, hypothesize, that falls into the infer workflow. We will use this method for conducting hypothesis tests.

But first, we need to initialize the test, which we will save as obs_diff.

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

## Response: weight (numeric)
## Explanatory: physical_3plus (factor)
## # A tibble: 1 x 1
##    stat
##   <dbl>
## 1  1.77

Notice how we can use the functions specify and calculate again like we did for calculating confidence intervals. Here, though, the statistic we are searching for is the difference in means, with the order being yes - no != 0.

After we have initialized the test, we need to simulate the test on the null distribution, which we will save as null.

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"))
null_dist

## Response: weight (numeric)
## Explanatory: physical_3plus (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 x 2
##    replicate    stat
##        <int>   <dbl>
##  1         1 -0.0818
##  2         2  0.218 
##  3         3 -0.0576
##  4         4  0.312 
##  5         5  0.229 
##  6         6  0.397 
##  7         7 -0.0132
##  8         8  0.156 
##  9         9  0.0601
## 10        10 -0.0137
## # ... with 990 more rows

Here, hypothesize is used to set the null hypothesis as a test for independence, i.e., that there is no difference between the two population means. In one sample cases, the null argument can be set to point to test a hypothesis relative to a point estimate.

Also, note that the type argument within generate is set to permute, which is the argument when generating a null distribution for a hypothesis test.

We can visualize this null distribution with the following code:

ggplot(data = null_dist, aes(x = stat)) +
  geom_histogram()

Now that the test is initialized and the null distribution formed, we can visualise to see how many of these null permutations have a difference of at least obs_stat of 1.77?

We can also calculate the p-value for our hypothesis test using the function infer::get_p_value().

null_dist %>% visualize() +
  shade_p_value(obs_stat = obs_diff, direction = "two-sided")

null_dist %>%
  infer::get_p_value(obs_stat = obs_diff, direction = "two_sided")

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

This the standard workflow for performing hypothesis tests.

Loading the data

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

## 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 ?

The data frame omega contains the salaries for the sample of 50 executives in the company. Can we conclude that there is a significant difference between the salaries of the male and female executives?

Note that we can perform different types of analyses, and check whether they all lead to the same conclusion

. Confidence intervals . Hypothesis testing . Correlation analysis . Regression

Calculate summary statistics on salary by gender. Also, create and print a dataframe where, for each gender, we show the mean, SD, sample size, the t-critical, the SE, the margin of error, and the low/high endpoints of a 95% condifence interval

# 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

We 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.

# 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 between group female and group male 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_2 <- omega %>%
  specify(salary ~ gender) %>%
  calculate(stat = "diff in means", order = c("female", "male"))

obs_diff_2

## Response: salary (numeric)
## Explanatory: gender (factor)
## # A tibble: 1 x 1
##     stat
##    <dbl>
## 1 -8696.

null_dist_2 <- omega %>%
  # specify variables
  specify(salary ~ gender) %>%
  
  # 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("female", "male"))

null_dist_2

## Response: salary (numeric)
## Explanatory: gender (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 x 2
##    replicate    stat
##        <int>   <dbl>
##  1         1 -1704. 
##  2         2  2626. 
##  3         3   469. 
##  4         4   519. 
##  5         5 -1423. 
##  6         6   824. 
##  7         7   -92.9
##  8         8  -357. 
##  9         9  -126. 
## 10        10 -1374. 
## # ... with 990 more rows

null_dist_2 %>% visualize()

ggplot(data = null_dist_2, aes(x = stat)) +
  geom_histogram()

null_dist_2 %>% visualize() +
  shade_p_value(obs_stat = obs_diff_2, direction = "two-sided")

null_dist_2 %>%
  infer::get_p_value(obs_stat = obs_diff_2, direction = "two_sided")

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

According to our calculation by t.test, our p-value for the hypothesis test is 0.0001651, which is smaller than 0.05 and hence statistically significant. The p-value we retrieved from our simulations with the infer package is also infinitesmal and statistically significant. There is strong evidence against the null hypothesis and the -8696.29 years (64542.84-73239.13) of difference we estimated in our sample means is really different from zero, therefore there is a significant difference.

We can safely conclude the salary levels of male executives is, on average, 8696.29 units more than female counterparts.

Relationship Experience - Gender?

At the board meeting, someone raised the issue that there was indeed a substantial difference between male and female salaries, but that this was attributable to other reasons such as differences in experience. A questionnaire send out to the 50 executives in the sample reveals that the average experience of the men is approximately 21 years, whereas the women only have about 7 years experience on average (see table below).

# 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

Based on this evidence, can we conclude that there is a significant difference between the experience of the male and female executives?

# hypothesis testing using t.test() 
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 between group female and group male 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

# hypothesis testing using infer package
obs_diff_3 <- omega %>%
  specify(experience ~ gender) %>%
  calculate(stat = "diff in means", order = c("female", "male"))

obs_diff_3

## Response: experience (numeric)
## Explanatory: gender (factor)
## # A tibble: 1 x 1
##    stat
##   <dbl>
## 1 -13.7

null_dist_3 <- omega %>%
  # specify variables
  specify(experience ~ gender) %>%
  
  # 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("female", "male"))

null_dist_3

## Response: experience (numeric)
## Explanatory: gender (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 x 2
##    replicate   stat
##        <int>  <dbl>
##  1         1 -3.80 
##  2         2 -4.12 
##  3         3 -2.12 
##  4         4  3.33 
##  5         5  0.923
##  6         6 -1.64 
##  7         7  3.09 
##  8         8  0.763
##  9         9 -0.920
## 10        10  0.442
## # ... with 990 more rows

null_dist_3 %>% visualize()

ggplot(data = null_dist_3, aes(x = stat)) +
  geom_histogram()

null_dist_3 %>% visualize() +
  shade_p_value(obs_stat = obs_diff_3, direction = "two-sided")

null_dist_3 %>%
  infer::get_p_value(obs_stat = obs_diff_3, direction = "two_sided")

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

According to our calculation by t.test, our p-value for the hypothesis test is 1.225e-05, which is smaller than 0.05 and hence statistically significant. The p-value we retrieved from our simulations with the infer package is also infinitesmal and statistically significant. There is strong evidence against the null hypothesis and the -13.7 years (7.38-21.13) of difference we estimated in our sample means is really different from zero, therefore there is a significant difference.

We can safely conclude the experience levels of male executives is, on average, 13.7 years more than female counterparts.

Relationship Salary - Experience ?

Someone at the meeting argues that clearly, a more thorough analysis of the relationship between salary and experience is required before any conclusion can be drawn about whether there is any gender-based salary discrimination in the company.

Analyse the relationship between salary and experience. Draw a scatterplot to visually inspect the data

ggplot(omega, aes(x = experience, y = salary)) + geom_point() +geom_smooth(se = FALSE)

In general, there is a moderate positive relationship between executives’ salary and experience. However, the strength of this relationship becomes weaker as the experience level increases, suggesting that experience may be an important factor influencing the salary levels for early-careers and not so much important for the already experienced professionals.

Check correlations between the data

We can use GGally:ggpairs() to create a scatterplot and correlation matrix. Essentially, we change the order our variables will appear in and have the dependent variable (Y), salary, as last in our list. We then pipe the dataframe to ggpairs() with aes arguments to colour by gender and make ths plots somewhat transparent (alpha = 0.3).

omega %>% 
  select(gender, experience, salary) %>% #order variables they will appear in ggpairs()
  ggpairs(aes(colour=gender, alpha = 0.3))+
  theme_bw()

The salary vs experience scatterplot shows that female executives in this sample are more concentrated on the left side of the plot, which means that generally they have lower experience levels than their male counterparts. Given the 0.812 correlation between salary and experience levels for female executives, which is higher than the 0.661 correlation for males, it is reasonable that we see a stronger relationship between salary and experience levels for early career professionals.

However, we cannot safely conclude whether this difference in salary level between female and male executives is due to gender difference or difference in experience levels. More research is required to investigate the causal relationship on salary levels and the two different factors.