1.0.1 - Beginner - Introduction to Statistics using R
introduction
review
econ 226
econ 227
descriptive statistics
beginner
importing data
data wrangling
t-test
inferential statistics
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.
Author
COMET Team William Co, Anneke Dresselhuis, Jonathan Graves, Emrul Hasan, Jonah Heyl, Mridul Manas, Shiming Wu.
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)
# A tibble: 5 × 100
pefamid pweight pagemieg PAS1MRAG PAS1MRG1 PAS1MRG2 PASR1MFA PASR1MR pasrbuyg
<chr> <dbl> <chr> <dbl+lbl> <dbl+lb> <dbl+lb> <dbl+lb> <dbl+l> <dbl+lb>
1 00001 4415. 06 NA NA NA NA 2 [No] 7 [200…
2 00002 805. 14 NA NA NA NA NA NA
3 00003 1393. 12 NA NA NA NA NA 4 [198…
4 00004 217. 11 NA NA NA NA NA 5 [199…
5 00005 455. 11 NA NA NA NA NA 5 [199…
# ℹ 91 more variables: pasrcnmg <dbl+lbl>, pasrcon <dbl+lbl>,
# pasrcst <dbl+lbl>, pasrcurg <dbl+lbl>, PASRDPO1 <dbl+lbl>,
# PASRDPO2 <dbl+lbl>, PASRDPO3 <dbl+lbl>, PASRDPO4 <dbl+lbl>,
# PASRDPO5 <dbl+lbl>, pasrdwng <dbl+lbl>, pasrfnmg <dbl+lbl>,
# pasrint <dbl+lbl>, pasrintg <chr>, pasrmoag <dbl+lbl>, pasrmpfg <dbl+lbl>,
# pasrmryg <dbl+lbl>, pasrrntg <dbl+lbl>, pasrskp <dbl+lbl>,
# pattcrc <dbl+lbl>, pattcrlm <dbl+lbl>, pattcrr <dbl+lbl>, …
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
starting httpd help server ... done
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 columnsSFS_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 columnsdf_gender_on_wealth <- SFS_data[keep]# previewhead(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.
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
Summary Statistics
Variable
N
Mean
Std. Dev.
Min
Pctl. 25
Pctl. 75
Max
gender
10201
... 1
6132
60%
... 2
4069
40%
education
10201
... 1
1216
12%
... 2
2279
22%
... 3
2819
28%
... 4
3827
38%
... 9
60
1%
wealth
10201
1030591
1744118
-875500
110750
1270250
30055000
income_before_tax
10201
95625
125418
-47000
27000
130000
3300000
income_after_tax
10201
86570
80987
-479300
39000
114250
2159100
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?
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by gender and education.
ℹ Output is grouped by gender.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(gender, education))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
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 paletteheatmap_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 heatmapheatmap_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.
Order rows using column values
Keep distinct/unique rows
Keep rows that match a condition
Get a glimpse of your data
Create, modify, and delete columns
Keep or drop columns using their names and types
Count the observations in each group
Group by one or more variables
A general vectorised if-else
mutate()
glimpse()
filter()
case_when()
select()
group_by()
distinct()
arrange()
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 comparisont_test_result <-t.test(income_before_tax ~ gender, data = df_gender_on_wealth)print(t_test_result)
Welch Two Sample t-test
data: income_before_tax by gender
t = 15.057, df = 10167, p-value < 2.2e-16
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
30119.32 39135.46
sample estimates:
mean in group 1 mean in group 2
109437.54 74810.15
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 comparisonmannwhitneyu_test_result <-wilcox.test(income_before_tax ~ gender, data = df_gender_on_wealth)print(mannwhitneyu_test_result)
Wilcoxon rank sum test with continuity correction
data: income_before_tax by gender
W = 15307599, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0
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)
[1] No No No No No Yes No Yes Yes Yes
Levels: No Yes
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.
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by university and gender.
ℹ Output is grouped by university.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(university, gender))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
results
# A tibble: 4 × 4
# Groups: university [2]
university gender m_wealth sd_wealth
<fct> <dbl+lbl> <dbl> <dbl>
1 No 1 [Male] 817055. 1459054.
2 No 2 [Female] 668021. 1447673.
3 Yes 1 [Male] 1680511. 2346053.
4 Yes 2 [Female] 1213004. 1563383.
f <-ggplot(data = SFS_data, aes(x = gender, y = wealth)) +xlab("Gender") +ylab("Wealth") # label and define our x and y axisf <- f +geom_bar(stat ="summary", fun ="mean", fill ="lightblue") # produce a summary statistic, the meanf <- f +facet_grid(. ~ university) # add a grid by educationf
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)
Warning: package 'scales' was built under R version 4.4.2
# A tibble: 2 × 2
# Groups: university [2]
university wealth_gap
<fct> <chr>
1 No -18%
2 Yes -28%
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 datat2 =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
Welch Two Sample t-test
data: filter(university_data, gender == 1)$wealth and filter(university_data, gender == 2)$wealth
t = 7.3872, df = 3782.4, p-value = 1.836e-13
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
343428.9 591586.0
sample estimates:
mean of x mean of y
1680511 1213004
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
Welch Two Sample t-test
data: filter(nuniversity_data, gender == 1)$wealth and filter(nuniversity_data, gender == 2)$wealth
t = 3.9821, df = 5192, p-value = 6.925e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
75662.69 222403.86
sample estimates:
mean of x mean of y
817054.8 668021.5
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 diplomaless_than_high_school_data <-filter(SFS_data, education ==1) # Less than high schoolhigh_school_data <-filter(SFS_data, education ==2) # High schoolpost_secondary_data <-filter(SFS_data, education ==3) # Non-university post-secondaryuniversity_data <-filter(SFS_data, education ==4) # UniversityretHS =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
Welch Two Sample t-test
data: filter(high_school_data, gender == 1)$income_before_tax and filter(less_than_high_school_data, gender == 1)$income_before_tax
t = 11.44, df = 1968.1, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
28213.26 39888.30
sample estimates:
mean of x mean of y
81191.49 47140.72
retHSF
Welch Two Sample t-test
data: filter(high_school_data, gender == 2)$income_before_tax and filter(less_than_high_school_data, gender == 2)$income_before_tax
t = 7.0649, df = 846.02, p-value = 3.359e-12
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
20217.37 35772.52
sample estimates:
mean of x mean of y
51472.26 23477.31
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.
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.
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-scoreremove_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 removalzscore_threshold <-1.645# Adjust as needed# remove outliers based on z-score for the desired variabledf_filtered <-remove_outliers_zscore(df_gender_on_wealth, "wealth", zscore_threshold)df_filtered
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.↩︎
