Data in proprietary formats

• The PISA data are provided in proprietary formats (SAS and SPSS).
• This means that the data are stored in a particular encoding scheme, designed so that decoding and reading the data is accomplished by particular software or hardware.
• In R, you can use the haven package to import the PISA data.

library(tidyverse)
library(haven)
pisa2018 <- read_sav(here::here("data", "CY07_MSU_STU_QQQ.sav")) %>% 
                as_factor()  # swap code and labels for labelled factors
dim(pisa2018)

## [1] 612004   1118

• Since the data is big, it will take a while to read the data in.
• Every row corresponds to a student.

Are girls worse in maths than boys?

• The gender gap in mathematics is a common discussion, with the concern being that girls tend to score lower than boys on average in standardized math tests.
• The PISA data provides an opportunity to explore the gender gap across numerous countries.
• In the pisa2018 data, the sex of the student is in variable ST004D01T and the country/region is in variable CNT.
• Let's rename these to sensible names, e.g. sex and country.
• We'll also modify some country names so that it can be joined with the map data later.
• We will focus on using PV1MATH and will not cover any analysis that require us to use all 10 plausible values in this course.

Bootstrap using R: Part 1

There are a number of ways of doing this in R but we will use sample_n function in dplyr 📦.

set.seed(2020) # for reproducibility
boot_sample1 <- pisa2018c %>% 
  group_by(country, sex) %>% 
  sample_n(size = n(), replace = TRUE)

We can then treat boot_sample1 as we did before to obtain another set of estimates for the gender gap for mathematics score by country.

boot_sample1 %>%
  summarise(avg = weighted.mean(PV1MATH, SENWT)) %>%
  ungroup() %>%
  pivot_wider(country, names_from = sex, values_from = avg) %>%
  mutate(diff = Female - Male, country = fct_reorder(country, diff))

## # A tibble: 76 x 4
##    country                Female  Male   diff
##    <fct>                   <dbl> <dbl>  <dbl>
##  1 Albania                  438.  434.   3.74
##  2 Argentina                373.  389. -15.4 
##  3 Australia                490.  495.  -4.57
##  4 Austria                  497.  508. -11.1 
##  5 Baku                     415.  424.  -9.02
##  6 Belarus                  470.  472.  -1.83
##  7 Belgium                  502.  517. -14.9 
##  8 Bosnia and Herzegovina   402.  410.  -7.43
##  9 Brazil                   376.  391. -14.6 
## 10 Brunei                   435.  425.   9.41
## # … with 66 more rows

Bootstrap using R: Part 2

• We need to repeat this process a reasonable number of times.
• We will do 100 times.
• To make this process easier, we will use the map_dfr function from purrr 📦.

boot_ests <- map_dfr(1:100, ~{
    pisa2018c %>%
      group_by(country, sex) %>%
      sample_n(size = n(), replace = TRUE) %>%
      summarise(avg = weighted.mean(PV1MATH, SENWT)) %>%
      ungroup() %>%
      pivot_wider(country, names_from = sex, values_from = avg) %>%
      mutate(diff = Female - Male, country = fct_reorder(country, diff)) %>%
      mutate(boot_id = .x)
  })

• The .x is substituted from an element from the first argument in map_dfr.

Data Generating Process: Simple Linear Model

• Unlike typical data, the data generating process is known for simulated data.
• For example, we can generate a set of observation that is only linearly dependent on an independent continuous variable.
• Mathematically, the data generating process is a simple linear model: $$\color{blue}{y_{i}} = \color{red}{\beta_0} + \color{red}{\beta_1}\color{blue}{x_i} + e_{i},$$ where for $i = 1, ..., n$, $y_i$ is the $i$-th response, $x_i$ is the corresponding covariate, $e_i$ is a random error, $\beta_0$ is the intercept and $\beta_1$ is the slope.
• Typically we assume that $e_i$s are independent and $e_i \sim N(0, \sigma^2)$.
• In practice, we only observe those colored in blue and we estimate those colored in red (for simplicity assuming the data generating process is known but the model parameters are unknown).

Simulated Data in R: Simple Linear Model

Suppose that $n = 200$, $\beta_0 = 3$, $\beta_1 = -2$ and $\sigma^2 = 1$.

set.seed(2020) # for reproducibility
n <- 200 # sample size
b0 <- 3  # intercept
b1 <- -2 # slope
sim_df <- tibble(id = 1:n) %>%  # initialise data set
    mutate(x = runif(n(), 0, 10), # draw x from Uniform[0,10]
           y = b0 + b1 * x + rnorm(n(), 0, 1))

Obtain least squares estimates (or maximum likelihood estimate) for $\beta_0$ and $\beta_1$:

fit1 <- lm(y ~ x, data = sim_df)
coef(fit1)

## (Intercept)           x 
##    2.975177   -2.020419

The estimates are pretty close to the true values!

Simulated Data in R: Simple Linear Model

• Suppose now that $e_i \sim \text{Gamma}(2, 1) - 2$.
• How good are the estimates under least squares when the error is not normally distributed?

sim2_df <- tibble(id = 1:n) %>%  
    mutate(x = runif(n(), 0, 10), 
           y = b0 + b1 * x + rgamma(n(), 2, 1) - 2)
fit2 <- lm(y ~ x, data = sim2_df)
coef(fit2)

## (Intercept)           x 
##    2.818092   -1.997902

• The estimate of the slope is still good but the estimate of the intercept is not as good as before.
• In the tutorial, you will also learn to simulate data when the dependent variable is a categorical variable 🔧

Null data

$$H_0: e_i \sim N(0, \hat{\sigma}^2)\qquad \text{vs.} \qquad H_1: \text{not }H_0$$

• $\hat{\sigma}$ is the estimate of $\sigma$ from the model fit.

c(summary(fit1)$sigma, summary(fit2)$sigma)

## [1] 1.140534 1.229640

• We simulate observations of size $n = 200$ from $N(0, \hat{\sigma}^2)$. These are called null data as it is generated under the null hypothesis.
• We draw a scatter plot with these observations against the original $x$-values.
• We repeat this process $K = 19$ times then produce a lineup of the plots with the original residual plot randomly included into it.