Examining COVID-19 Vaccination Rates in California
How are race/ethnicity and income related to vaccination rates?
Manas Khatore, Sahana Noru, Nikki Trueblood, and Sean Yang
from datascience import *
import pandas as pd
import seaborn as sns
sns.set_theme(style="ticks", color_codes=True)
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import array
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn import metrics
from sklearn.metrics import mean_squared_error
Association between race/ethnicity and COVID-19 vaccination rates
raw_table = Table().read_table('covid19vaccinesbycountybydemographic.csv')
raw_table
race_table = raw_table.where('demographic_category', 'Race/Ethnicity').drop('county_type', 'demographic_category')
race_table
The different race/ethnicity categories in "race_table" are American Indian or Alaska Native, Asian, Black or African American, Latino, Multiracial, Native Hawaiian or Other Pacific Islander, Other Race, Unknown and White. Since the table is sorted according to the administration date of the vaccine, we decided to clean the table by only keeping the rows of the most recent administration date (11/12/2021).
race_and_county_vax = race_table.where("administered_date", "2021-11-12")
race_and_county_vax
# adding the percent of population fully vaccinated and with one dose of the vaccine to the table
race_and_county_vax = race_and_county_vax.relabeled('est_age_5plus_pop', 'eligible population')
fully_vax_percent = 100*race_and_county_vax.column('cumulative_fully_vaccinated')/race_and_county_vax.column('eligible population')
one_dose_percent = 100*race_and_county_vax.column('cumulative_unvax_total_pop')/race_and_county_vax.column('eligible population')
race_and_county_vax = race_and_county_vax.with_columns('percent_fully_vaccinated', fully_vax_percent, 'percent_one_dose', one_dose_percent)
race_and_county_vax
Looking at the table, we noticed that many values in the "Other Race" and "Unknown" categories are nan. Checking the percentages in the table that are over 100:
above_100 = race_and_county_vax.where('percent_fully_vaccinated', are.above(100))
above_100.show()
Since the "Native Hawaiian or Other Pacific Islander" demographic largely contains percentages that are over 100, we decided to fully remove it from the table, as well as removing "Other Race" and "Unknown."
# removing demographics and percentages higher than 100
# converting to a DataFrame for visualization purposes
final_race_and_county_vax = race_and_county_vax.where('demographic_value', are.not_equal_to('Native Hawaiian or Other Pacific Islander')).where('demographic_value', are.not_equal_to('Other Race')).where('demographic_value', are.not_equal_to('Unknown')).where('percent_fully_vaccinated', are.below_or_equal_to(100)).where('percent_one_dose', are.below_or_equal_to(100))
final_race_and_county_vax = final_race_and_county_vax.select("county", "demographic_value", "est_population", "eligible population", "cumulative_fully_vaccinated", "cumulative_at_least_one_dose", "percent_fully_vaccinated", "percent_one_dose")
final_race_and_county_vax_df = final_race_and_county_vax.to_df()
final_race_and_county_vax_df
%matplotlib inline
sns.catplot(x="percent_fully_vaccinated", y="demographic_value", data=final_race_and_county_vax_df)
# black dots/diamonts in the boxplot represent outliers
sns.catplot(x='percent_fully_vaccinated', y='demographic_value', kind='box', data=final_race_and_county_vax_df)
Seeing the differences in fully vaccinated rates along the different demographic groups, we decided to perform an A/B test to determine whether certain groups were getting vaccinated at statistically significantly higher rates than others. Instead of analyzing the differences in each combination of groups, we decided to split the demographics into "traditionally represented" and "traditionally underrepresented" racial groups, using sources about presence in higher education.
represented = final_race_and_county_vax.where('demographic_value', are.not_equal_to("American Indian or Alaska Native")).where('demographic_value', are.not_equal_to("Black or African American")).where('demographic_value', are.not_equal_to("Latino"))
represented = represented.group('county', sum)
new_percent_full = represented.column('cumulative_fully_vaccinated sum') / represented.column('eligible population sum') * 100
new_percent_one_dose = represented.column('cumulative_at_least_one_dose sum') / represented.column('eligible population sum') * 100
represented = represented.drop('percent_fully_vaccinated sum', 'percent_one_dose sum').with_columns('percent_fully_vaccinated', new_percent_full, 'percent_at_least_one_dose', new_percent_one_dose).drop('demographic_value sum').where('county', are.not_equal_to('Statewide'))
labels = []
for i in range(represented.num_rows):
labels = np.append(labels, 'represented')
represented = represented.with_column('label', labels)
represented
underrepresented = final_race_and_county_vax.where('demographic_value', are.not_equal_to("Asian")).where('demographic_value', are.not_equal_to("Multiracial")).where('demographic_value', are.not_equal_to("White"))
underrepresented = underrepresented.group('county', np.sum)
new_percent_full = underrepresented.column('cumulative_fully_vaccinated sum') / underrepresented.column('eligible population sum') * 100
new_percent_one_dose = underrepresented.column('cumulative_at_least_one_dose sum') / underrepresented.column('eligible population sum') * 100
underrepresented = underrepresented.drop('percent_fully_vaccinated sum', 'percent_one_dose sum').with_columns('percent_fully_vaccinated', new_percent_full, 'percent_at_least_one_dose', new_percent_one_dose).drop('demographic_value sum').where('county', are.not_equal_to('Statewide'))
labels = make_array()
for i in np.arange(underrepresented.num_rows):
labels = np.append(labels, "underrepresented")
underrepresented = underrepresented.with_column('label', labels)
underrepresented
A/B Testing
# we shuffle the labels and define the function test_stat() that will calculate test statistics
# we calculate the observed test statistic from the data
represented_df = represented.to_df()
underrepresented_df = underrepresented.to_df()
frames = [represented_df, underrepresented_df]
rep_vs_underrep = pd.concat(frames)
shuffled_labels = rep_vs_underrep.sample(n = 58)['label']
def test_stat(tbl):
represented_values = tbl.loc[(tbl['label'] == 'represented')]
underrepresented_values = tbl.loc[(tbl['label'] == 'underrepresented')]
return represented_values['percent_fully_vaccinated'].mean() - underrepresented_values['percent_fully_vaccinated'].mean()
observed_difference = test_stat(rep_vs_underrep)
observed_difference
# one simulation of the A/B test
def one_simulated_difference():
shuffled_labels = rep_vs_underrep.sample(n = 116)['label']
shuffled_table = rep_vs_underrep.drop(columns = ['label'])
shuffled_table['label'] = shuffled_labels.tolist()
return test_stat(shuffled_table)
one_simulated_difference()
# completing 10000 simulations and storing the results in an array
differences = []
for i in range(10000):
new_difference = one_simulated_difference()
differences = np.append(differences, new_difference)
differences_df = pd.DataFrame(differences, columns = ['Differences'])
differences_df
# plotting the results of the simulations
# the red dot represents the observed test statistic
differences_df.hist()
plt.scatter(observed_difference, -0.002, color='red', s=40);
above_observed = differences_df.loc[differences_df['Differences'] >= observed_difference]
above_observed.astype(bool).sum(axis=0)/10000
Our calculated p-value is 0.0011. Using a standard p-value cutoff of 0.05, we would reject the null hypothesis. From our analysis, we would conclude that traditionally underrepresented groups are being vaccinated a lower rates than other groups in California.
sns.catplot(x="percent_one_dose", y="demographic_value", data=final_race_and_county_vax_df)
# black dots/diamonts represent outliers
sns.catplot(x='percent_one_dose', y='demographic_value', kind='box', data=final_race_and_county_vax_df)
The above visualizations show that some groups that had low rates of full vaccination have high rates of one dose vaccination, such as the American Indian/Alaska Native and Multiracial demographics. This can potentially be explained by these groups getting vaccinated in larger numbers more recently, or that members in these groups only received one does and did not return for the second.
Linear Regression between COVID-19 vaccination rates and median income per county
# importing a table that contains the number of people who are fully vaccinated or received at least one dose per county
vaccines = Table().read_table('covid19vaccinesbycounty.csv')
new_vaccines = vaccines.where('county', are.not_equal_to('All CA Counties')).where('county', are.not_equal_to('All CA and Non-CA Counties')).where('county', are.not_equal_to('Outside California')).where('county', are.not_equal_to('Unknown'))
new_vaccines = new_vaccines.group('county', sum).select('county', 'fully_vaccinated sum', 'at_least_one_dose sum')
new_vaccines
# importing a table that contains the median income per county in 2018
# adding a column with the total population in each county
# using the total population and the new_vaccines table to calculate the percent of fully vaccinated and at least one dose per county
income_table = Table().read_table('incomebycounty.csv')
int_values = income_table.apply(lambda x: int((x[0:x.index(',')] + x[x.index(',') + 1:])), 'Value')
income_table = income_table.drop('Value').with_column("Median income", int_values)
population_array = make_array(1682253, 1204, 40474, 211632,45292,21839,1165927, 27743,191185,1008654,28917,136463, 179702, 19016,909235,152486,68163,32730,10014009,156255,262321,17131,91601,281202,8700,13195,439035,138019,102241,3186989, 404739,19790,2418185,1585055,64209,2181654,3298634,873965,779233,282424,764442,448229,1936259,270861,182155,3236,44076,453491,488863,552878,99633,65829,16112,473117,55620,843843,216403,81575)
income_table = income_table.with_column('Total population', population_array)
full_vax = 100*new_vaccines.column('fully_vaccinated sum')/income_table.column('Total population')
one_dose = 100*new_vaccines.column('at_least_one_dose sum')/income_table.column('Total population')
income_table = income_table.with_columns('Fully vaccinated percentage', full_vax, 'At least one dose percentage', one_dose)
income_table_df = income_table.to_df()
income_table_df
# plotting fully vaccinated percentage vs. median income
fig, ax = plt.subplots(figsize = (20,10))
plt.scatter(data = income_table_df, x = "Median income", y = "Fully vaccinated percentage")
plt.xlabel("Median income")
plt.ylabel("Fully vaccinated percentage")
income_and_fully_vax = income_table_df[["Median income", "Fully vaccinated percentage"]]
sns.regplot(income_and_fully_vax['Median income'], income_and_fully_vax['Fully vaccinated percentage'])
print("The correlation coefficient is: " + str(income_and_fully_vax.corr(method ='pearson').iloc[0, 1]))
# plotting at least one dose percentage vs. median income
fig, ax = plt.subplots(figsize = (20,10))
plt.scatter(data = income_table_df, x = "Median income", y = "At least one dose percentage")
plt.xlabel("Median income") # You need to manually add these in!
plt.ylabel("At least one dose percentage")
income_and_one_dose = income_table_df[["Median income", "At least one dose percentage"]]
sns.regplot(income_and_one_dose['Median income'], income_and_one_dose['At least one dose percentage'])
print("The correlation coefficient is: " + str(income_and_one_dose.corr(method ='pearson').iloc[0, 1]))
The graphs above show that there is a strong, positive correlation between median income and the percentage of people who are fully vaccinated or have at least one dose in each county. Seeing that there is a correlation, we decided to create linear regression models that would be able to predict the percentage of vaccinated people given the median income of an area.
# Linear model predicting fully vaccinated percentage from median income
train_data, test_data = train_test_split(income_table_df, test_size = 0.20, random_state = 50)
X_train = train_data[['Median income']]
y_train = train_data['Fully vaccinated percentage']
X_test = test_data[['Median income']]
y_test = test_data['Fully vaccinated percentage']
model = LinearRegression(fit_intercept=True)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
print("The coefficient of determination is " + str(metrics.r2_score(y_test, y_pred)))
print("The mean squared error is " + str(mean_squared_error(y_test, y_pred)))
# Linear model predicting one dose percentage from median income
X_train = train_data[['Median income']]
y_train = train_data['At least one dose percentage']
X_test = test_data[['Median income']]
y_test = test_data['At least one dose percentage']
model = LinearRegression(fit_intercept=True)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
print("The coefficient of determination is " + str(metrics.r2_score(y_test, y_pred)))
print("The mean squared error is " + str(mean_squared_error(y_test, y_pred)))
The coefficient of determination shows the proportion of outcomes that can be explained by the model. In order to reduce the mean squared error for the models, we decided to standardize the units by converting each measure to the number of standard deviations it is away from its mean.
standardized_income = (income_table_df["Median income"]-income_table_df["Median income"].mean())/income_table_df["Median income"].std()
standardized_full_vax = (income_table_df["Fully vaccinated percentage"]-income_table_df["Fully vaccinated percentage"].mean())/income_table_df["Fully vaccinated percentage"].std()
standardized_one_dose = (income_table_df["At least one dose percentage"]-income_table_df["At least one dose percentage"].mean())/income_table_df["At least one dose percentage"].std()
standards = {"Median income": standardized_income, "Fully vaccinated percentage": standardized_full_vax, "At least one dose percentage": standardized_one_dose}
standardized_income_table_df = pd.DataFrame(data=standards)
train_data, test_data = train_test_split(standardized_income_table_df, test_size = 0.20, random_state = 50)
X_train = train_data[['Median income']]
y_train = train_data['Fully vaccinated percentage']
X_test = test_data[['Median income']]
y_test = test_data['Fully vaccinated percentage']
model = LinearRegression(fit_intercept=True)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
print("The coefficient of determination is " + str(metrics.r2_score(y_test, y_pred)))
print("The mean squared error is " + str(mean_squared_error(y_test, y_pred)))
X_train = train_data[['Median income']]
y_train = train_data['At least one dose percentage']
X_test = test_data[['Median income']]
y_test = test_data['At least one dose percentage']
model = LinearRegression(fit_intercept=True)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
print("The coefficient of determination is " + str(metrics.r2_score(y_test, y_pred)))
print("The mean squared error is " + str(mean_squared_error(y_test, y_pred)))