Exploring Diversity in California Schools

Importing the Data

data <- readr::read_tsv("filesenr.asp.txt")
head(data) %>% kable()

longlat_raw <- read_tsv("pubschls.txt")
longlat <- longlat_raw %>%
  select(CDSCode, Latitude, Longitude, School, City)
head(longlat) %>% kable()

url <- 'https://en.wikipedia.org/wiki/List_of_California_locations_by_income'
income_data <- read_html(url) %>%
  html_nodes("table") %>%
  .[2] %>%
  html_table() %>%
  .[[1]] %>%
  set_names(c("county", "population", "population_density", 
              "per_capita_income", "median_household_income",
              "median_family_income"))
head(income_data) %>% kable()

income_data$population <- as.character(income_data$population)
income_data$population <- parse_number(income_data$population, locale = locale(grouping_mark = ","))

income_data$population_density <- as.character(income_data$population_density)
income_data$population_density <- parse_number(income_data$population_density, locale = locale(grouping_mark = ","))

income_data$per_capita_income <- as.character(income_data$per_capita_income)
income_data$per_capita_income <- parse_number(income_data$per_capita_income, locale = locale(grouping_mark = ","))

income_data$median_household_income <- as.character(income_data$median_household_income)
income_data$median_household_income <- parse_number(income_data$median_household_income, locale = locale(grouping_mark = ","))

income_data$median_family_income <- as.character(income_data$median_family_income)
income_data$median_family_income <- parse_number(income_data$median_family_income, locale = locale(grouping_mark = ","))

head(income_data) %>% kable()

The Diversity Index

To be able to compare the ethnic diversity between schools we need a measure that describes diversity to begin with. I’ll use a diversity index developed in 1991 by a researcher at the University of North Carolina at Chapel Hill. The index simply asks: “What is the probability that two people in a population picked at random will be from different ethnicity?”

To make sense of this measure, we will look at a plot:

Above, we have a population of 10 people split into two ethnicities, with 5 in each. The colored squares represent random picks with the same ethnicity, light grey squares are picks with different ethnicities, and dark grey indicates impossible picks (the same person being picked twice). Now the diversity index can be calculated by dividing the number of light grey squares with the sum of light grey and colored squares.

First, lets go back to the enrollemnt data. We would like to have it in a different format where each row represents a school. Another thing we want to do is look across each ethnic category within each school. The new dataset that we will be working with is:

small_data <- data %>%
  select(CDS_CODE, ETHNIC, ENR_TOTAL) %>%
  group_by(CDS_CODE, ETHNIC) %>%
  summarize(ENR_TOTAL = sum(ENR_TOTAL)) %>%
  spread(ETHNIC, ENR_TOTAL, fill = 0) %>%
  ungroup()
head(small_data) %>% kable()

Before we move on, lets reference the data documentation to see what each of the ethnic numbers mean. According to the File Structure:

• Code 0 = Not reported
• Code 1 = American Indian or Alaska Native, Not Hispanic
• Code 2 = Asian, Not Hispanic
• Code 3 = Pacific Islander, Not Hispanic
• Code 4 = Filipino, Not Hispanic
• Code 5 = Hispanic or Latino
• Code 6 = African American, not Hispanic
• Code 7 = White, not Hispanic
• Code 9 = Two or More Races, Not Hispanic

Now we have two decisions to make before moving on. How shoule we take care of the “Not reported” cases, and “Two or More Races, Not Hispanic”? (These are the types of decisions you make during the EDA phase, which then set up your analysis later on)

Lets start with the case of code 0. We notice that the category generally is quite small compared to the other groups, so we need to reason about what would happen if we drop them. If we treat them as a separate group, it would inflate the diversity index so we will drop them.

Dealing with multi-ethnicity is slighlty harder. We have four options:

• Assume everyone in “Two or More Races, Not Hispanic” is the same ethnicity
• Assume everyone in “Two or More Races, Not Hispanic” is all different ethnicities
• Ignore the group “Two or More Races, Not Hispanic”
• Distribute the group “Two or More Races, Not Hispanic” evenly into the other groups

Options 3 and 4 will both yield low diversity indices considering that the choice of “Two or More Races, Not Hispanic” is because these students are different enough not to be put into one of the more precisely defined groups. The first option, while appealling from a computational standpoint, would treat a black-white mixture and an Asian-American Indian as members of the same group even though they have no racial identity in common. Therefore, we will work with the second option which might overestimate the diversity of any given population as it ignores people with “Two or More Races, Not Hispanic” with identical mixes (siblings would be a prime example).

After getting that out of the way, we can thus create a function to calculate the diversity index by using the geometric interpretation from before.

diversity <- function(..., y) {
 x <- sapply(list(...), cbind)
 total <- cbind(x, y)
 1 - (rowSums(x ^ 2) - rowSums(x)) / (rowSums(total) ^ 2 - rowSums(total))
}

data_diversity <- small_data %>%
  mutate(diversity = diversity(`1`, `2`, `3`, `4`, `5`, `6`, `7`, y = `9`))

Now we are ready to bring back the data and start to calculate diversity indices for the schools. Lets run a summary to see the ranges we get and make sure everything is working correctly.

summary(data_diversity$diversity)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.0000  0.3140  0.5205  0.4711  0.6425  1.0000     108

All the indices are between 0 and 1 which is good! This means that the diversity index worked as intended. However there are 108 schools that have NA valued diversity. Lets look at those schools:

filter(data_diversity, is.na(diversity)) %>% head() %>% kable()

It seems that these schools do not have missing values but rather have such a low number of kids that the calculations break down. We will exclude these schools for now.

data_diversity <- data_diversity %>%
  filter(!is.na(diversity))

The first thing we should do with the data is look at the distribution of the diversity indices.

ggplot(data_diversity, aes(x=diversity)) +
  geom_histogram(binwidth = 0.01) +
  xlab("Diversity Index") +
  ylab(NULL) +
  ggtitle("Diversity at California Schools") +
  theme(plot.title = element_text(hjust = 0.5))

We notice that most of the schools follow a nice shape with a median above 0.5. The only weird thing are the spikes at 0 and 1. Lets look at the 0’s first:

data_diversity %>%
  filter(diversity == 0) %>%
  head(20) %>% kable()

It seems that all of the schools with a 0 diversity index is because they have students from one ethnicity, which is predominantly category 5 which is “Hispanic or Latino”. Now lets look at the 1’s:

data_diversity %>%
  filter(diversity == 1) %>%
  head(20) %>% kable()

Intuitively, the only way a school can have a diveristy index of 1 is if the maximum number of students from each ethnicity is 1. Therefore, these schools are rather small.

Lets take a look at the total enrollment in the schools. We have seen a couple of instances with low enrollment and we wouldn’t want those to distort our view.

total_students <- data %>%
  select(CDS_CODE:SCHOOL, ENR_TOTAL) %>%
  group_by(CDS_CODE) %>%
  summarise(ENR_TOTAL = sum(ENR_TOTAL)) %>%
  ungroup()

meta <- data %>%
  select(CDS_CODE:SCHOOL) %>%
  distinct() 

ggplot(total_students, aes(ENR_TOTAL)) +
  geom_histogram(bins = 1000) +
  xlab("Total Enrollment") +
  ylab(NULL) +
  ggtitle("Total Enrollment at California Schools") +
  theme(plot.title = element_text(hjust = 0.5))

We see a hump at the 600-700 mark, but a very large spike closer to 0. We should take a closer look at these small values to see what is happening.

total_students %>%
  filter(ENR_TOTAL < 250) %>%
  ggplot(aes(ENR_TOTAL)) +
  geom_histogram(bins = 250) +
  xlab("Total Enrollment") +
  ylab(NULL) +
  ggtitle("Total Enrollment at Smaller California Schools") +
  theme(plot.title = element_text(hjust = 0.5))

Here, we will make another decision that will change the results of the analysis, but we will restrict the investigation to only include schools with total enrollments of 50 or larger. We do this so that these smaller schools do not skew the distribution of the diversity index.

data_big <- total_students %>%
  left_join(meta, by = "CDS_CODE") %>%
  filter(ENR_TOTAL > 50)
data_enrollment <- data_diversity %>%
  right_join(data_big, by = "CDS_CODE")

One thing we can look at is the relationship between enrollment and diversity. We will fit a simple linear model to start with. I expect larger schools to have more diversity, but lets see:

data_enrollment %>%
  ggplot(aes(ENR_TOTAL, diversity)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = "lm") +
  theme_minimal() +
  ylim(0,1) +
  labs(x="Enrollment", y="Diversity Index", title = "Relationship between Enrollment and Diversity") +
  theme(plot.title = element_text(hjust = 0.5))

The slope seems quite small but still about a 5% increase in diversity per 2000 students. Lets check the null hypothesis to make sure for statistical significance.

data_enrollment %>%
  mutate(ENR_TOTAL = ENR_TOTAL / 1000) %>%
  lm(diversity ~ ENR_TOTAL, data = .) %>%
  summary()

## 
## Call:
## lm(formula = diversity ~ ENR_TOTAL, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.51078 -0.15751  0.05026  0.17040  0.39595 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.449955   0.003429 131.223  < 2e-16 ***
## ENR_TOTAL   0.024282   0.004063   5.976 2.37e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2136 on 9472 degrees of freedom
## Multiple R-squared:  0.003756,   Adjusted R-squared:  0.003651 
## F-statistic: 35.71 on 1 and 9472 DF,  p-value: 2.371e-09

The results show that we can reject the null hypothesis of no correlation and we get a better estimate of the slope of 2.3%. Its important to keep in mind that we have a miniscule R-squared, so we dont think there is much of a relationship here.

Does this have anything to do with money?

I will end this post by bringing in the income data and checking if there is a correlation between median household income in the county in which the school is located and the schools’ diversity indices.

data_enrollment %>%
  left_join(income_data, by = c("COUNTY" = "county")) %>%
  ggplot(aes(x=median_household_income, diversity)) +
  geom_jitter(alpha = 0.1, width = 500) +
  geom_smooth(method = "lm") +
  ylim(0,1) +
  labs(x="Median Household Income", y = "Diversity Index", title = "Correlation between Diversity Index and Median Household Income \n(in CA counties)") +
  theme(plot.title = element_text(hjust = 0.5))

As always, if you have a question or a suggestion related to the topic covered in this article, please feel free to contact me!