1 Introduction

The Oregon Department of Education tracks each incoming group of high school freshmen—called a cohort—over the course of their time in high school and calculates the percentage of those students who graduate within four years.² This graduation rate has been increasing steadily since the 2008-2009 school year.³

This data can be found in ODE’s Graduation Reports, which also break down each cohort by school, and then by several subcategories such as gender, race/ethnicity, and status as economically disadvantaged. However, it’s not possible to see, for example, how many Hispanic/Latino students who graduated are economically disadvantaged, which might help explain why they have a higher or lower graduation rate compared to other races/ethnicities. And while the reports do calculate graduation rates for students with/without disabilities, students who are/aren’t English learners, students who are/aren’t homeless, etc., we can’t determine the graduation rates for students in multiple such subgroups or for students in none of them.

On December 31, 2018, ODE released a cross-tabulated dataset for the 2016-17 cohort which contains counts for cohort and graduates from every combination of the following categories:

• Racial and ethnic group
• Gender
• Status as an English learner
• Status as a student with disabilities
• Status as homeless
• Status as economically disadvantaged

Using this new data we can attempt to disentangle these factors from one another. The questions we will focus on are:

• What are the differences in graduation rates between males and females…
• What are the effects of the four burdens at the end of the list above on graduation rates…
• How do the graduation rates of different races and ethnicities compare…

…when all other factors in the above list are equal?

# Libraries and helper functions we'll need.

library(viridisLite)
library(readxl)
library(loo)
library(coda)
library(rstan)

options(mc.cores = 4)  # number of cores in your processor

cols <- viridis(10)
turq <- cols[5]

# From the 'rethinking' package (https://github.com/rmcelreath/rethinking)
col.alpha <- function( acol , alpha=0.2 ) {
    acol <- col2rgb(acol)
    acol <- rgb(acol[1]/255, acol[2]/255, acol[3]/255, alpha)
    acol
}

2 The data

The cross-tabulated dataset can be downloaded from ODE’s School and District Profiles: Accountability Measures page under “Cohort Graduation Rates / Dropout Rates” with the title “Cross-Tabulated 2016-17 Four-year Cohort Graduation Rates”.

To protect the privacy of the students, some of the data in this file has been suppressed.⁴ If any group contains fewer than 10 students, an asterisk (*) is entered instead of the number of students in the group. So if there are fewer than 10 students of a certain group who graduate, there will be an * instead of the number of graduates. And if there are fewer than 10 students who do not graduate, there will be an * instead of the number of students in the cohort.

Rather than throw these rows away, we will show in section 3 how they can still be used in a Bayesian model—even the rows where both the size of the cohort and the number of graduates are suppressed. In fact we will argue that ignoring those double-suppressed rows can bias the inferred graduation rates toward the extremes of 0% and 100%.

Before loading the spreadsheet into R, a quick look reveals that whenever the number of graduates is suppressed, so is the size of the cohort. This quirk of this particular dataset will allow us to write a simpler model than if there were rows with unsuppressed cohort sizes but suppressed numbers of graduates.

Now let’s load the spreadsheet. If the cohort size is suppressed but the number of graduates isn’t, we’ll put a 1 for the cohort size. If both cohort size and number of graduates is suppressed, we’ll put a 0 for the cohort size. This will let us treat the two cases differently in the analysis.

ctab <- read_xlsx(
    "ESSA_STATE_Cohort_Graduation_Cross_Tabulated_1718_SP.xlsx",
    sheet = 2,
    na = "*"
)
ctab <- ctab[,c(-1, -10)]
for (i in seq_len(nrow(ctab))) {
    if (is.na(ctab$`Adjusted Cohort`[i]))
        if (is.na(ctab$`Number of Graduates`[i])) {
            ctab$`Adjusted Cohort`[i] <- 0
            ctab$`Number of Graduates`[i] <- 0
        } else
            ctab$`Adjusted Cohort`[i] <- 1
}
ctab

And format it for the models, giving it short column names and using 0 and 1 for the binary categories. We’ll also assign each race/ethnicity a unique ID.

df <- data.frame(
    female = as.integer(ifelse(ctab$Gender == "Female", 1, 0)),
    race_ethn = as.factor(ctab$`Race/Ethnicity`),
    eng_learn = as.integer(ifelse(ctab$`English Learners in High School` == "Y", 1, 0)),
    disability = as.integer(ifelse(ctab$`Students with Disabilities` == "Y", 1, 0)),
    homeless = as.integer(ifelse(ctab$Homeless == "Y", 1, 0)),
    econ_disad = as.integer(ifelse(ctab$`Economically Disadvantaged` == "Y", 1, 0)),
    cohort = as.integer(ctab$`Adjusted Cohort`),
    graduates = as.integer(ctab$`Number of Graduates`)
)
df$race_ethn_id <- as.integer(df$race_ethn)
df

3 The heart of the models

We will model the number of students who graduate as a binomial random variable,

\[ \textrm{graduates}_i \sim \operatorname{Binomial}(\textrm{cohort}_i, \theta). \]

Of course, we can only use this likelihood if we know both the size of the cohort and the number of graduates.

dbinom(16, 16, 1)

## [1] 1

plot(
    0,
    xlim = c(0, 1), ylim = c(0, 1),
    type = "n",
    xlab = "graduation rate theta", ylab = "density",
    main = "Posterior distributions p(theta | graduates=16, cohort=n)",
    yaxt = "n"
)
axis(2, 0)
for (i in 16:25)
    curve(dbinom(16, i, x), col = cols[i - 15], lwd = 2, n = 301, add = TRUE)
for (i in 1:10) {
    lines(c(0, 0.1), rep(1 - 0.05*i, 2), col = cols[i], lwd = 2)
    text(0.14, 1 - 0.05*i, labels = paste("n =", i + 15))
}