1.0.1 - Beginner - Introduction to Statistics using R

This notebook is an introduction to basic statistics using Jupyter and R, and some fundamental data analysis. It is a high-level review of the most important applied tools from earlier units.

Outline

Prerequisites

• Introduction to Jupyter
• Introduction to R

References

• Esteban Ortiz-Ospina and Max Roseqoifhoihr (2018) - “Economic inequality by gender”. Published online at OurWorldInData.org. Retrieved from: https://ourworldindata.org/economic-inequality-by-gender [Online Resource]

Outcomes

In this notebook, you will learn how to:

• Import data from the Survey of Financial Security (Statistics Canada, 2019)
• Wrangle, reshape and visualize SFS_data as part of an Exploratory Data Analysis (EDA)
• Run statistical tests, such as the $t$-test, to compare mean income of male-led vs. female-led households
• Generate summary statistics tables and other data-representations of the data using group_by()
• Optional: Run a formal two sample t-test to check for heterogeneity in how gender affects income and compare the returns to education

Part 1: Import Data into R

The data we use comes from the 2019 Survey of Financial Security released by Statistics Canada ¹.

# run this cell to load necessary packages for this tutorial
# install.packages('vtable')
# install.packages('viridis')
library(tidyverse)
library(haven)
library(dplyr)
library(vtable)
library(viridis)


source("beginner_intro_to_statistics2_tests.r")
source("beginner_intro_to_statistics2_functions.r")
# warning messages are okay

The tidyverse is a collection of R packages developed by Hadley Wickham and his colleagues as a cohesive set of tools for data manipulation, visualization, and analysis. In a tidy data set, each variable forms a column and each observation forms a row. tidyverse packages such as the tidyr and dplyr are recommended for cleaning and transforming your data into tidy formats.

Let’s import the .dta file from Statistics Canada using the read_dta function.

# if this is your first time using Jupyter Lab, the shortcut to run a cell is `Shift + Enter`
SFS_data <- read_dta("../datasets_beginner/SFS_2019_Eng.dta")

Here are some of the common file extensions and import functions in R:

• .dta and read_dta() for STATA files
• .csv and read_csv() for data stored as comma-separated values
• .Rda and load() for RStudio files and other files formatted for R

head(SFS_data, 5)

Note: head(df, n) displays the first n rows of the data frame. Other popular methods include glance() and print().

# you can read the documentation for a given function by adding a question-mark before its name
?head

Part 2: Exploratory Data Analysis in R

There are a routine of steps you should generally follow as part of your EDA or Exploratory Data Analysis. Normally, you would analyze and visualize the variation, correlation, and distribution of your variables of interest. We do this to gain an intuitive understanding of the data before we undertake any formal hypothesis tests or model-fitting.

Let’s think of our key variables of interest. We’re interested in estimating the effect of gender on differences in earnings.

• Independent variable: gender of the highest income earner

• Variable of interest: income after tax for each individual

• Variable of interest: income before tax for each individual

• Control: wealth for the household

• Control: level of education

Cleaning and Reshaping SFS_data

For now, it’d be convenient to work with a new data frame containing only the key variables (columns) listed above. Moreover, the columns need to be renamed so they are easier for the reader to remember.

# rename columns
SFS_data <- SFS_data %>%
    rename(income_before_tax = pefmtinc) %>% 
    rename(income_after_tax = pefatinc) %>%
    rename(wealth = pwnetwpg) %>%
    rename(gender = pgdrmie) %>%
    rename(education = peducmie) 

# drop rows where tax info is missing, ie. pefmtinc = 'NA'.
SFS_data <- filter(SFS_data, !is.na(SFS_data$income_before_tax))

keep <- c("pefamid", "gender", "education", "wealth", "income_before_tax", "income_after_tax")

# new df with chosen columns
df_gender_on_wealth <- SFS_data[keep]

# preview
head(df_gender_on_wealth, 5)

Note: This is another tidy representation of the original data but with less variables. The original data set is still stored as SFS_data.

Ensuring correct data-types

Notice that education is stored as chr but we want to keep it as a factor. The variable education came encoded as it is from a set of values {1, 2, 3, 4, 9}, each of which represent a level of education obtained.

df_gender_on_wealth <- df_gender_on_wealth %>%
         mutate(education = as.factor(education), 
         gender = as.factor(gender),
         income_before_tax = as.numeric(income_before_tax),
         income_after_tax = as.numeric(income_after_tax))

head(df_gender_on_wealth, 2)

All good! Let’s use descriptive statistics to understand how each of the numbers in the set {1, 2, 3, 9} represent an individual’s educational background.

Computing Descriptive Statistics using vtable in R

Let’s calculate the summary statistics of our dataset.

Note: the sumtable method from the vtable package can be used to display the table in different formats including LaTeX, HTML, and data.frame.

# out = "kable" tells it to return a knitr::kable()
# replace "kable" with "latex" and see what happens!
sumtbl <- sumtable(df_gender_on_wealth, out = "kable")
sumtbl

This is like having a birds-eye view of our data. As a researcher, we should take note of outliers and other irregularities and ask how those issues might affect the validity of our models and tests.

Note: see Appendix for a common method to remove outliers using Z-score thresholds.

Grouping observations

Wouldn’t it be neat to see how mean or median incomes for male and female-led households look like based on the level of education obtained by the main income-earner?

by_gender_education <- df_gender_on_wealth %>%
  group_by(gender, education) %>%
  summarise(mean_income = mean(income_before_tax, na.rm = TRUE),
            median_income = median(income_before_tax, na.rm = TRUE),
            mean_wealth = mean(wealth, na.rm = TRUE),
            median_wealth = median(wealth, na.rm = TRUE))

by_gender_education

Note: this is again a tidy representation of SFS_data. Grouping observations by gender and education makes it a bit easier to make comparisons across groups.

We can take this chain-of-thought further and generate a heatmap using the ggplot package.

library(ggplot2)
library(viridis)

# Create the heatmap with an accessible color palette
heatmap_plot <- ggplot(by_gender_education, aes(x = education, y = gender, fill = mean_income)) +
  geom_tile() +
  scale_fill_viridis_c(option = "plasma", na.value = "grey", name = "Mean Income") +
  labs(x = "Education", y = "Gender")

# Display the heatmap
heatmap_plot

Note: we use scale_fill_viridis_c() from the viridis package to ensure that the color palette follows the standards of DS.

Now, what does this tell you about how male-led households (gender = 1) compare with female-led households in terms of the mean household income? Does this tell if education widens the income gap between male-led and female-led households with the same level of education?

We can infer from the visualization that the female-led households with the same level of education have different mean incomes as compared to male-led households. This smells of heterogeneity and we can explore regression and other empirical methods to formally test this claim.

However, we shouldn’t yet draw any conclusive statements about the relationships between gender (of the main income earner), income, education and other variables such as wealth.

As researchers, we should ask if the differences in the mean or median incomes for the two groups are significant at all. We can then go a bit further and test if education indeed widens the gap or not.

Think Deeper: how would you specify the null and alternative hypotheses?

Test your knowledge

Match the function with the appropriate description. Enter your answer as a long string with the letter choices in order.

1. Order rows using column values
2. Keep distinct/unique rows
3. Keep rows that match a condition
4. Get a glimpse of your data
5. Create, modify, and delete columns
6. Keep or drop columns using their names and types
7. Count the observations in each group
8. Group by one or more variables
9. A general vectorised if-else

1. mutate()
2. glimpse()
3. filter()
4. case_when()
5. select()
6. group_by()
7. distinct()
8. arrange()
9. count()

Note: it’s fine if you don’t know all those functions yet! Match the functions you know and run code to figure out the rest.

# Enter your answer as a long string ex: if you think the matches are 1-B, 2-C, 3-A, enter answer as "BCA"
answer_1 <- ""

test_1()

Part 3: Running $t$-tests in R

Let’s run a t-test for a comparison of means.

# performs a t-test for means comparison
t_test_result <- t.test(income_before_tax ~ gender, data = df_gender_on_wealth)
print(t_test_result)

The 95% confidence interval does not include 0 and we can confirm that the male-led households on average earn more as income before tax than the female-led households, and the gap is statistically significant.

Let’s now run a test to compare the medians of both groups.

# perform a Mann-Whitney U test for median comparison
mannwhitneyu_test_result <- wilcox.test(income_before_tax ~ gender, data = df_gender_on_wealth)
print(mannwhitneyu_test_result)

This p-value is again highly significant, and based on our data, the median incomes for the two groups are not equal.

Our variable of interest is income, and so far, we have provided statistical evidence for the case that the gender of the main income-earner is correlated with the household’s income.

We are however more interested in the causal mechanisms through which education and wealth determine how gender affects household income.

Think Deeper: According to Ortiz-Ospina and Roser (2018), women are overrepresented in low-paying jobs and are underrepresented in high-paying ones. What role does the attainment of education play in sorting genders into high vs. low-paying jobs? Can we test this formally with the data?

Studying how wealth and education might impact the income-gap

There are multiple reasons to study the links between wealth and the income gap. For instance, we might want to answer whether having more wealth affects an individual’s income.

We can use some of the methods we have learned in R to analyze and visualize relationships between income, gender, education and wealth.

Let’s see if having a university degree widens the gender income gap.

SFS_data <- SFS_data %>% 
            mutate(university = case_when(     # create a new variable with mutate
                            education == "4" ~ "Yes",    # use case_when and ~ operator to applt `if else` conditions 
                            TRUE ~ "No")) %>% 
            mutate(university = as_factor(university)) #remember, it's a factor!

head(SFS_data$university, 10)

Let’s visualize how the mean wealth compares for male-led vs. female-led households, conditional on whether the main-income earner went to university.

results <- SFS_data %>%
           group_by(university,gender) %>%
           summarize(m_wealth = mean(wealth), sd_wealth = sd(wealth))

results 

f <- ggplot(data = SFS_data, aes(x = gender, y = wealth)) + xlab("Gender") + ylab("Wealth")    # label and define our x and y axis
f <- f + geom_bar(stat = "summary", fun = "mean", fill = "lightblue")    # produce a summary statistic, the mean
f <- f + facet_grid(. ~ university)    # add a grid by education

f

It smells like the wealth gap between the two types of households widens for groups that have obtained a university degree.

Similarly, let’s look at the difference in wealth gap in percentage terms. We use results generated in the previous cell (the $4\times 4$ table) as the inputs this time. We need to load the package scales to use the function percent.

library(scales)

percentage_table <- SFS_data %>%
                    group_by(university) %>%
                    group_modify(~ data.frame(wealth_gap = mean(filter(., gender == 2)$wealth)/mean(filter(., gender == 1)$wealth) - 1)) %>%
                    mutate(wealth_gap = scales::percent(wealth_gap))

percentage_table

Notice the signs are both negative. Hence, on average, female-led households have less wealth regardless of whether they have a university degree or not.

More importantly, based on our data, female-led households with university degrees on average have 28% less wealth than male-led households with university degrees. Comparing the two groups given they don’t have university degrees, the gap is quite smaller: 18%.

So, we have shown that the gap widens by about 10% when conditioned for a university degree.

Let’s test this further by creating sub-samples of “university degree” and “no university degree” respectively and then running a formal two sample t-test.

university_data <- filter(SFS_data, university == "Yes") # university only data 
nuniversity_data <- filter(SFS_data, university == "No") # non university data

t2 = t.test(
       x = filter(university_data, gender == 1)$wealth,
       y = filter(university_data, gender == 2)$wealth,
       alternative = "two.sided",
       mu = 0,
       conf.level = 0.95)

t2  # test for the wealth gap in university data

round(t2$estimate[1] - t2$estimate[2],2) # rounds our estimate


t3 = t.test(
       x = filter(nuniversity_data, gender == 1)$wealth,
       y = filter(nuniversity_data, gender == 2)$wealth,
       alternative = "two.sided",
       mu = 0,
       conf.level = 0.95)

t3 # test for the wealth gap in non-university data

round(t3$estimate[1] - t3$estimate[2],2) # rounds our estimate

In both tests, the p-values are very small, indicating strong statistical evidence to reject the null hypothesis. The confidence intervals also provide a range of plausible values for the difference in means, further supporting the alternative hypothesis.

Based on these results, there appears to be a significant difference in wealth between the two gender groups regardless of university-status, with males consistently having higher mean wealth compared to females.

Optional: Returns to HS diploma

Next, examine whether returns to education differ between genders. For our purposes, we will define returns to education as the difference in average income before tax between two subsequent education levels.

The following t-test finds the returns to education of a high school diploma for males (retHS) and for females(retHSF).

# Returns to education: High school diploma

less_than_high_school_data <- filter(SFS_data, education == 1) # Less than high school
high_school_data <- filter(SFS_data, education == 2) # High school
post_secondary_data <- filter(SFS_data, education == 3) # Non-university post-secondary
university_data <- filter(SFS_data, education == 4) # University


retHS = t.test(
       x = filter(high_school_data, gender == 1)$income_before_tax,
       y = filter(less_than_high_school_data, gender == 1)$income_before_tax,
       alternative = "two.sided",
       mu = 0,
       conf.level = 0.95)
retHS_ans=round(retHS$estimate[1] - retHS$estimate[2],2)

retHSF = t.test(
       x = filter(high_school_data, gender == 2)$income_before_tax,
       y = filter(less_than_high_school_data, gender == 2)$income_before_tax,
       alternative = "two.sided",
       mu = 0,
       conf.level = 0.95)

retHS
retHSF
retHS_ans=round(retHS$estimate[1] - retHS$estimate[2],2)
retHSF_ans=round(retHSF$estimate[1] - retHSF$estimate[2],2)

We have found statistically significant evidence for the case that returns to graduating with a high school diploma are indeed positive for individuals living in both male-led and female-led households.

Test your knowledge

As an exercise, create a copy of the cell above and try to calculate the returns of a university degree for males and females.

# your code here

Now let’s work with a simulated dataset of mutual fund performance. Interpret the data below as the yearly returns for a sample of 300 mutual funds from 2010 to 2015.

fund_performance

Create a subset of the data with the returns for 2015. Rename the column to investment_returns. Store your answer in fp_15.

# feel free to use this cell if you need 

fp_15 <- fund_performance %>%
          ...(fund, ...)

answer_2 <- fp_15

test_2()

Calculate the mean and median return of funds in 2015. Store your answers in mean_ret and median_ret, respectively.

# feel free to use this cell if you need 

mean_ret <- ...(...)

answer_3 <- mean_ret

test_3()

median_ret <- ...(...)

answer_4 <- median_ret

test_4()

Let’s suppose the market return (average return of investments available) was 5%. Run a 95% confidence level t-test on the returns of fp_15 to find whether the funds outperformed the market or not. Complete the code below.

t_stat = ...( 
       ...,
       mu = ...,
       alternative = "two.sided",
       conf.level = ...)

answer_5 <- t_stat$conf.int

test_5()

Do we have statistical evidence to believe that the funds outperformed the market?

1. Yes
2. No

# enter your answer as either "A" or "B"
answer_6 <- ""

test_6()

But wait! Do you notice anything interesting about this dataset? Investigate the dataset with special attention to the NAs. Do you notice a pattern?

fund_performance

There are no funds with negative performance in the dataset! It’s likely that the NAs have replaced the observations with negative returns. How might that affect our analysis of fund performance? Think about the biases that could have been introduced in our mean and statistical test calculations.

Appendix

Removing outliers is a common practice in data analysis. The code below removes outliers based on a custom Z-score threshold.

Note: here we use the 95th percentile but you should first visualize your data with box plots and then find a convenient threshold to remove outliers in the variables of interest.

# function to remove outliers based on z-score
remove_outliers_zscore <- function(data, variable, threshold) {
  z_scores <- scale(data[[variable]])
  data_without_outliers <- data[abs(z_scores) <= threshold, ]
  return(data_without_outliers)
}

# set the threshold for z-score outlier removal
zscore_threshold <- 1.645    # Adjust as needed

# remove outliers based on z-score for the desired variable
df_filtered <- remove_outliers_zscore(df_gender_on_wealth, "wealth", zscore_threshold)

df_filtered

Footnotes

1. Statistics Canada, Survey of Financial Security, 2019, 2021. Reproduced and distributed on an “as is” basis with the permission of Statistics Canada. Adapted from Statistics Canada, Survey of Financial Security, 2019, 2021. This does not constitute an endorsement by Statistics Canada of this product.