Water Quality in Central Valley
Yueyang Tang & Hyeyoon Cho
2026-05-10
setup
Import SAFER data
#Importing safer data
safer <- read.csv("SAFER_RA.csv")Hyeyoon’s part about k-12 schools
Then I wanted to find out how the systems are classified, so I put this code.
Water System Classifications
unique(safer$FEDERAL_CLASSIFICATION_TYPE)## [1] "COMMUNITY" "NON-TRANSIENT NON-COMMUNITY"
## [3] "TRANSIENT NON-COMMUNITY"The result it gave me – meaning there are 3 categories of public water systems.
And I wanted to find out how K-12 schools are classified according to these 3 categories, so I entered the code below.
The result was given below. I was able to find out K-12 schools are predominantly categorized as non-transient non-community.
Filtering K-12 Schools
schools <- safer %>%
filter(grepl("SCHOOL", SYSTEM_NAME, ignore.case = TRUE))
schools %>%
count(FEDERAL_CLASSIFICATION_TYPE)## FEDERAL_CLASSIFICATION_TYPE n
## 1 COMMUNITY 3
## 2 NON-TRANSIENT NON-COMMUNITY 333
## 3 TRANSIENT NON-COMMUNITY 4And I wanted to find out how many water systems are categorized as non-transient non-community as a whole.
Identifying At-Risk Schools
safer %>%
filter(FEDERAL_CLASSIFICATION_TYPE == "NON-TRANSIENT NON-COMMUNITY") %>%
nrow()## [1] 379schools_ntnc <- safer %>%
filter(FEDERAL_CLASSIFICATION_TYPE == "NON-TRANSIENT NON-COMMUNITY")So there are 333 K-12 schools, out of 379 non-transient non-community public water systems.
Since most non-transient non-community systems are K-12 schools, I wanted to find out which schools have the worst water quality. I used the “at risk” indicator.
schools_ntnc_atrisk <- schools_ntnc %>%
filter(RISK_ASSESSMENT_RESULT == "At-Risk")And I wanted to see this in a bar graph, by counties.
schools_ntnc_atrisk %>%
count(COUNTY) %>%
ggplot(aes(x = reorder(COUNTY, n), y = n)) +
geom_bar(stat = "identity", fill = "steelblue") +
coord_flip() +
labs(title = "NTNC Schools At-Risk by County",
x = "County",
y = "Number of Systems") +
theme_minimal()
### Correlation Analysis: Failing Schools Correlation between water
quality of already failing schools and median household income,
affordability, and operational capacity of a water system
schools_ntnc_failing <- schools_ntnc %>%
filter(FINAL_SAFER_STATUS == "Failing")
schools_ntnc_failing %>%
select(MHI, WATER_QUALITY_SCORE, WEIGHTED_WATER_QUALITY_SCORE,
AFFORDABILITY_SCORE, TMF_CAPACITY_SCORE) %>%
mutate(across(everything(), as.numeric)) %>%
cor(use = "pairwise.complete.obs") %>%
round(2)## MHI WATER_QUALITY_SCORE
## MHI 1.00 -0.33
## WATER_QUALITY_SCORE -0.33 1.00
## WEIGHTED_WATER_QUALITY_SCORE -0.33 1.00
## AFFORDABILITY_SCORE -0.58 -0.01
## TMF_CAPACITY_SCORE -0.01 -0.03
## WEIGHTED_WATER_QUALITY_SCORE AFFORDABILITY_SCORE
## MHI -0.33 -0.58
## WATER_QUALITY_SCORE 1.00 -0.01
## WEIGHTED_WATER_QUALITY_SCORE 1.00 -0.01
## AFFORDABILITY_SCORE -0.01 1.00
## TMF_CAPACITY_SCORE -0.03 -0.01
## TMF_CAPACITY_SCORE
## MHI -0.01
## WATER_QUALITY_SCORE -0.03
## WEIGHTED_WATER_QUALITY_SCORE -0.03
## AFFORDABILITY_SCORE -0.01
## TMF_CAPACITY_SCORE 1.00Findings
Relationship between MHI (median household income) and WATER_QUALITY_SCORE: -0.35 Among failing school water system, the lower the surrounding community’s income, the worse the water quality score. This is not a strong correlation but still meaningful.
The MHI and AFFORDABILITY_SCORE relationship at -0.60. Lower income communities are bearing a higher affordability burden, meaning they’re spending a larger share of their income on water even though the water is already failing.
MHI and TMF_CAPACITY_SCORE is basically zero (0.00), so income has no relationship with the technical and managerial capacity of the water system. This is interesting because it suggests the governance/management problems are spread across income levels equally
Correlation Analysis: All NTNC Schools
Correlation between water quality of all non-transient non-community schools and median household income, affordability, and operational capacity of a water system
schools_ntnc %>% select(MHI, WATER_QUALITY_SCORE, WEIGHTED_WATER_QUALITY_SCORE, AFFORDABILITY_SCORE, TMF_CAPACITY_SCORE) %>% mutate(across(everything(), as.numeric)) %>% cor(use = "pairwise.complete.obs") %>% round(2)## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `across(everything(), as.numeric)`.
## Caused by warning:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.## MHI WATER_QUALITY_SCORE
## MHI 1.00 -0.04
## WATER_QUALITY_SCORE -0.04 1.00
## WEIGHTED_WATER_QUALITY_SCORE -0.04 1.00
## AFFORDABILITY_SCORE -0.45 0.07
## TMF_CAPACITY_SCORE 0.02 0.05
## WEIGHTED_WATER_QUALITY_SCORE AFFORDABILITY_SCORE
## MHI -0.04 -0.45
## WATER_QUALITY_SCORE 1.00 0.07
## WEIGHTED_WATER_QUALITY_SCORE 1.00 0.07
## AFFORDABILITY_SCORE 0.07 1.00
## TMF_CAPACITY_SCORE 0.05 -0.03
## TMF_CAPACITY_SCORE
## MHI 0.02
## WATER_QUALITY_SCORE 0.05
## WEIGHTED_WATER_QUALITY_SCORE 0.05
## AFFORDABILITY_SCORE -0.03
## TMF_CAPACITY_SCORE 1.00Findings
MHI and water quality correlation is -0.05. The full population of NTNC schools income doesn’t really predict water quality. Thus, income doesn’t determine whether a school ends up with a bad water system in the first place.
But once a school is already failing, the poorer communities tend to have it worse within that failing group. Richer schools aren’t necessarily avoiding failure, but when they do fail they’re not as far down.
The affordability and water quality correlation is -0.47, which is noteworthy. Lower income communities are always paying more relative to what they earn, regardless of whether the system is failing or not.
Hyeyoon’s part about Failing Systems Across All Categories
Separately, I wanted to see, out of the 3 categories, which one has the most failing water system.
safer %>%
filter(FINAL_SAFER_STATUS == "Failing") %>%
count(FEDERAL_CLASSIFICATION_TYPE) %>%
ggplot(aes(x = reorder(FEDERAL_CLASSIFICATION_TYPE, n), y = n)) +
geom_bar(stat = "identity", fill = "tomato") +
coord_flip() +
labs(title = "Failing Water Systems by Federal Classification Type",
x = "Federal Classification Type",
y = "Number of Systems") +
theme_minimal()
#It was the community system. So this is when I thought we should pivot to the community system and analyze them, since they have far more numbers of failing water systems. Yaelle’s part
Data Overview
Before anything else, I want to know what fields are available andwhat the key categorical variables look like. How many system types are there, what are the possible risk statuses, etc.
names(safer)## [1] "TINWSYS_IS_NUMBER"
## [2] "WATER_SYSTEM_NUMBER"
## [3] "SYSTEM_NAME"
## [4] "REGULATING_AGENCY"
## [5] "COUNTY"
## [6] "FEDERAL_CLASSIFICATION_TYPE"
## [7] "SERVICE_CONNECTIONS"
## [8] "POPULATION"
## [9] "OWNER_TYPE"
## [10] "PL_ADDRESS"
## [11] "PL_ADDRESS_CITY_NAME"
## [12] "PL_ADDRESS_STATE_CODE"
## [13] "PL_ADDRESS_ZIP_CODE"
## [14] "LATITUDE_MEASURE"
## [15] "LONGITUDE_MEASURE"
## [16] "SERVICE_AREA_ECONOMIC_STATUS"
## [17] "MHI"
## [18] "CALENVIRO_SCREEN_SCORE"
## [19] "RISK_ASSESSMENT_RESULT"
## [20] "CURRENT_FAILING"
## [21] "FINAL_SAFER_STATUS"
## [22] "FAILING_START_DATE"
## [23] "PRIMARY_MCL_VIOLATION"
## [24] "PRIMARY_ANALYTES"
## [25] "SECONDARY_MCL_VIOLATION"
## [26] "SECONDARY_ANALYTES"
## [27] "E_COLI_VIOLATION"
## [28] "E_COLI_ANALYTES"
## [29] "TREATMENT_TECHNIQUE_VIOLATION"
## [30] "TT_ANALYTES"
## [31] "MONITORING_AND_REPORTING_VIOLATION"
## [32] "MONITORING_AND_REPORTING_ANALYTES"
## [33] "SOURCE_CAPACITY_VIOLATION"
## [34] "SOURCE_CAPACITY_ANALYTES"
## [35] "TOTAL_WEIGHTED_RISK_SCORE_BEFORE_DIVIDING_BY_CATEGORY_COUNT"
## [36] "WEIGHTED_WATER_QUALITY_SCORE"
## [37] "WATER_QUALITY_SCORE"
## [38] "WATER_QUALITY_PERCENTAGE_OF_TOTAL_RISK_SCORE"
## [39] "WATER_QUALITY_RISK_LEVEL"
## [40] "WEIGHTED_ACCESSIBILITY_SCORE"
## [41] "ACCESSIBILITY_SCORE"
## [42] "ACCESSIBILITY_PERCENTAGE_OF_TOTAL_RISK_SCORE"
## [43] "ASSESSIBILITY_RISK_LEVEL"
## [44] "WEIGHTED_AFFORDABILITY_SCORE"
## [45] "AFFORDABILITY_SCORE"
## [46] "AFFORDABILITY_PERCENTAGE_OF_TOTAL_RISK_SCORE"
## [47] "AFFORDABILITY_RISK_LEVEL"
## [48] "WEIGHTED_TMF_CAPACITY_SCORE"
## [49] "TMF_CAPACITY_SCORE"
## [50] "TMF_CAPACITY_PERCENTAGE_OF_TOTAL_RISK_SCORE"
## [51] "TMF_CAPACITY_RISK_LEVEL"
## [52] "HISTORY_OF_E_COLI_PRESENCE_RISK_LEVEL"
## [53] "HISTORY_OF_E_COLI_PRESENCE_THRESHOLD_MET"
## [54] "HISTORY_OF_E_COLI_PRESENCE_RAW_SCORE"
## [55] "HISTORY_OF_E_COLI_PRESENCE_WEIGHTED_SCORE"
## [56] "INCREASING_PRESENCE_OF_WATER_QUALITY_TRENDS_TOWARD_MCL_RISK_LEVEL"
## [57] "INCREASING_PRESENCE_OF_WATER_QUALITY_TRENDS_TOWARD_MCL_THRESHOLD_MET"
## [58] "INCREASING_PRESENCE_OF_WATER_QUALITY_TRENDS_TOWARD_MCL_RAW_SCORE"
## [59] "INCREASING_PRESENCE_OF_WATER_QUALITY_TRENDS_TOWARD_MCL_WEIGHTED_SCORE"
## [60] "TREATMENT_TECHNIUQE_VIOLATIONS_RISK_LEVEL"
## [61] "TREATMENT_TECHNIUQE_VIOLATIONS_THRESHOLD_MET"
## [62] "TREATMENT_TECHNIUQE_VIOLATIONS_RAW_SCORE"
## [63] "TREATMENT_TECHNIUQE_VIOLATIONS_WEIGHTED_SCORE"
## [64] "PAST_PRESENCE_ON_THE_FAILING_LIST_RISK_LEVEL"
## [65] "PAST_PRESENCE_ON_THE_FAILING_LIST_THRESHOLD_MET"
## [66] "PAST_PRESENCE_ON_THE_FAILING_LIST_RAW_SCORE"
## [67] "PAST_PRESENCE_ON_THE_FAILING_LIST_WEIGHTED_SCORE"
## [68] "CONSTITUENTS_OF_EMERGING_CONCERN_RISK_LEVEL"
## [69] "CONSTITUENTS_OF_EMERGING_CONCERN_THRESHOLD_MET"
## [70] "CONSTITUENTS_OF_EMERGING_CONCERN_RAW_SCORE"
## [71] "CONSTITUENTS_OF_EMERGING_CONCERN_WEIGHTED_SCORE"
## [72] "PERCENTAGE_OF_SOURCES_EXCEEDING_AN_MCL_RISK_LEVEL"
## [73] "PERCENTAGE_OF_SOURCES_EXCEEDING_AN_MCL_THRESHOLD_MET"
## [74] "PERCENTAGE_OF_SOURCES_EXCEEDING_AN_MCL_RAW_SCORE"
## [75] "PERCENTAGE_OF_SOURCES_EXCEEDING_AN_MCL_WEIGHTED_SCORE"
## [76] "NUMBER_OF_WATER_SOURCES_RISK_LEVEL"
## [77] "NUMBER_OF_WATER_SOURCES_THRESHOLD_MET"
## [78] "NUMBER_OF_WATER_SOURCES_RAW_SCORE"
## [79] "NUMBER_OF_WATER_SOURCES_WEIGHTED_RAW_SCORE"
## [80] "ABESENCE_OF_INTERTIES_RISK_LEVEL"
## [81] "ABESENCE_OF_INTERTIES_THRESHOLD_MET"
## [82] "ABESENCE_OF_INTERTIES_RAW_SCORE"
## [83] "ABESENCE_OF_INTERTIES_WEIGHTED_SCORE"
## [84] "SOURCE_CAPACITY_VIOLATION_RISK_LEVEL"
## [85] "SOURCE_CAPACITY_VIOLATIONS_THRESHOLD_MET"
## [86] "SOURCE_CAPACITY_VIOLATIONS_RAW_SCORE"
## [87] "SOURCE_CAPACITY_VIOLATIONS_WEIGHTED_SCORE"
## [88] "BOTTLED_WATER_OR_HAULED_WATER_RELIANCE_RISK_LEVEL"
## [89] "BOTTLED_WATER_OR_HAULED_WATER_RELIANCE_THRESHOLD_MET"
## [90] "BOTTLED_WATER_OR_HAULED_WATER_RELIANCE_RAW_SCORE"
## [91] "BOTTLED_WATER_OR_HAULED_WATER_RELIANCE_WEIGHTED_SCORE"
## [92] "DWR_DROUGHT_AND_WATER_SHORTAGE_RISK_ASSESSMENT_PERCENTILE_RISK_LEVEL"
## [93] "DWR_DROUGHT_AND_WATER_SHORTAGE_RISK_ASSESSMENT_PERCENTILE_THRESHOLD_MET"
## [94] "DWR_DROUGHT_AND_WATER_SHORTAGE_RISK_ASSESSMENT_PERCENTILE_RAW_SCORE"
## [95] "DWR_DROUGHT_AND_WATER_SHORTAGE_RISK_ASSESSMENT_PERCENTILE_WEIGHTED_SCORE"
## [96] "CRITICALLY_OVERDRAFTED_GROUNDWATER_BASIN_RISK_LEVEL"
## [97] "CRITICALLY_OVERDRAFTED_GROUNDWATER_BASIN_THRESHOLD_MET"
## [98] "CRITICALLY_OVERDRAFTED_GROUNDWATER_BASIN_RAW_SCORE"
## [99] "CRITICALLY_OVERDRAFTED_GROUNDWATER_BASIN_WEIGHTED_SCORE"
## [100] "PERCENT_OF_MEDIAN_HOUSEHOLD_INCOME_MHI_RISK_LEVEL"
## [101] "PERCENT_OF_MEDIAN_HOUSEHOLD_INCOME_MHI_THRESHOLD_MET"
## [102] "PERCENT_OF_MEDIAN_HOUSEHOLD_INCOME_MHI_RAW_SCORE"
## [103] "PERCENT_OF_MEDIAN_HOUSEHOLD_INCOME_MHI_WEIGHTED_SCORE"
## [104] "EXTREME_WATER_BILL_RISK_LEVEL"
## [105] "EXTREME_WATER_BILL_THRESHOLD_MET"
## [106] "EXTREME_WATER_BILL_RAW_SCORE"
## [107] "EXTREME_WATER_BILL_WEIGHTED_SCORE"
## [108] "HOUSEHOLD_SOCIOECONOMIC_BURDEN_RISK_LEVEL"
## [109] "HOUSEHOLD_SOCIOECONOMIC_BURDEN_THRESHOLD_MET"
## [110] "HOUSEHOLD_SOCIOECONOMIC_BURDEN_RAW_SCORE"
## [111] "HOUSEHOLD_SOCIOECONOMIC_BURDEN_WEIGHTED_SCORE"
## [112] "TOTAL_NET_ANNUAL_INCOME_RISK_LEVEL"
## [113] "TOTAL_NET_ANNUAL_INCOME_THRESHOLD_MET"
## [114] "TOTAL_NET_ANNUAL_INCOME_RAW_SCORE"
## [115] "TOTAL_NET_ANNUAL_INCOME_WEIGHTED_SCORE"
## [116] "OPERATING_RATIO_RISK_LEVEL"
## [117] "OPERATING_RATIO_THRESHOLD_MET"
## [118] "OPERATING_RATIO_RAW_SCORE"
## [119] "OPERATING_RATIO_WEIGHTED_SCORE"
## [120] "DAYS_CASH_ON_HAND_RISK_LEVEL"
## [121] "DAYS_CASH_ON_HAND_THRESHOLD_MET"
## [122] "DAYS_CASH_ON_HAND_RAW_SCORE"
## [123] "DAYS_CASH_ON_HAND_WEIGHTED_SCORE"
## [124] "OPERATOR_CERTIFICATION_VIOLATIONS_RISK_LEVEL"
## [125] "OPERATOR_CERTIFICATION_VIOLATIONS_THRESHOLD_MET"
## [126] "OPERATOR_CERTIFICATION_VIOLATIONS_RAW_SCORE"
## [127] "OPERATOR_CERTIFICATION_VIOLATIONS_WEIGHTED_SCORE"
## [128] "MONITORING_AND_REPORTING_VIOLATIONS_RISK_LEVEL"
## [129] "MONITORING_AND_REPORTING_VIOLATIONS_THRESHOLD_MET"
## [130] "MONITORING_AND_REPORTING_VIOLATIONS_SCORE"
## [131] "MONITORING_AND_REPORTING_VIOLATIONS_WEIGHTED_SCORE"
## [132] "SIGNIFICANT_DEFICIENCIES_RISK_LEVEL"
## [133] "SIGNIFICANT_DEFICIENCIES_THRESHOLD_MET"
## [134] "SIGNIFICANT_DEFICIENCIES_RAW_SCORE"
## [135] "SIGNIFICANT_DEFICIENCIES_WEIGHTED_SCORE"
## [136] "FUNDING_RECEIVED_SINCE_2017"
## [137] "TECHNICAL_ASSISTANCE_FUNDING_SINCE_2017"
## [138] "REGIONAL_BOARD"
## [139] "AUTOMATICALLY_AT_RISK_REASON"
## [140] "CREATED_DATE"unique(safer$FEDERAL_CLASSIFICATION_TYPE)## [1] "COMMUNITY" "NON-TRANSIENT NON-COMMUNITY"
## [3] "TRANSIENT NON-COMMUNITY"unique(safer$FINAL_SAFER_STATUS)## [1] "Not At-Risk" "At-Risk" "Not Assessed"
## [4] "Failing" "Potentially At-Risk"unique(safer$CURRENT_FAILING)## [1] "Not Failing" "Failing"unique(safer$RISK_ASSESSMENT_RESULT)## [1] "Not At-Risk" "At-Risk" "Not Assessed"
## [4] "Potentially At-Risk"Adding Binary Indicator Columns
To make filtering easier later, I added three 0/1 columns, one for each system type.
safer$ntnc <- as.integer(safer$FEDERAL_CLASSIFICATION_TYPE == "NON-TRANSIENT NON-COMMUNITY")
safer$community <- as.integer(safer$FEDERAL_CLASSIFICATION_TYPE == "COMMUNITY")
safer$tnc <- as.integer(safer$FEDERAL_CLASSIFICATION_TYPE == "TRANSIENT NON-COMMUNITY")Creating the Analysis Subset: safer_sub
The full dataset has a lot of columns I don’t need. I pulled out just the ones relevant to this analysis to keep things cleaner.
safer_sub <- safer[, c("TINWSYS_IS_NUMBER", "WATER_SYSTEM_NUMBER", "SYSTEM_NAME",
"REGULATING_AGENCY", "COUNTY", "FEDERAL_CLASSIFICATION_TYPE",
"SERVICE_CONNECTIONS", "POPULATION", "OWNER_TYPE",
"PL_ADDRESS", "PL_ADDRESS_CITY_NAME",
"LATITUDE_MEASURE", "LONGITUDE_MEASURE",
"SERVICE_AREA_ECONOMIC_STATUS", "MHI", "CALENVIRO_SCREEN_SCORE",
"RISK_ASSESSMENT_RESULT", "CURRENT_FAILING",
"WATER_QUALITY_SCORE", "WATER_QUALITY_PERCENTAGE_OF_TOTAL_RISK_SCORE",
"WATER_QUALITY_RISK_LEVEL",
"ntnc", "community", "tnc")]Which Type of Water System Is Failing the Most?
Now that I know the data structure, the first question is: are all three system types equally at risk, or is one type driving most of the failures?
Failure Rates by Water System Type
I looked at raw counts first, then row proportions to compare failure rates across types.
#raw counts, how many systems of each type are failing vs. not
current_fail <- table(safer$FEDERAL_CLASSIFICATION_TYPE, safer$CURRENT_FAILING)
current_fail##
## Failing Not Failing
## COMMUNITY 379 2434
## NON-TRANSIENT NON-COMMUNITY 46 333
## TRANSIENT NON-COMMUNITY 0 12prop.table(current_fail, margin = 1)##
## Failing Not Failing
## COMMUNITY 0.1347316 0.8652684
## NON-TRANSIENT NON-COMMUNITY 0.1213720 0.8786280
## TRANSIENT NON-COMMUNITY 0.0000000 1.0000000# Convert the counts into row percentages, compare failure RATES across system types
final_status <- table(safer$FINAL_SAFER_STATUS)
barplot(final_status)
Community systems make up the vast majority of the dataset (2,813 total vs. 379 NTNC and just 12 TNC). Of all currently failing systems, 379 out of 425 (about 89%) are community systems. That matters because community systems serve people where they live. Only 12 TNC systems appear in the data, which is a very small sample.
SAFER Status × System Type: Cross-Tab and Stacked Bar Chart
But raw counts can be misleading if one type just has more systems. So I built a cross-tab to check whether the composition stays the same across risk levels.
#I am trying to answer are Failing or At-Risk systems disproportionately one type, or does the breakdown stay roughly the same across risk levels?
counts <- table(safer$FINAL_SAFER_STATUS, safer$FEDERAL_CLASSIFICATION_TYPE)
props <- round(prop.table(counts, margin = 1), 2)
combined <- matrix(paste0(counts, " (", props, ")"),
nrow = nrow(counts),
dimnames = dimnames(counts))
combined <- cbind(combined, Total = rowSums(counts))
combined## COMMUNITY NON-TRANSIENT NON-COMMUNITY
## At-Risk "571 (0.87)" "87 (0.13)"
## Failing "379 (0.89)" "46 (0.11)"
## Not Assessed "143 (0.81)" "22 (0.12)"
## Not At-Risk "1331 (0.88)" "182 (0.12)"
## Potentially At-Risk "389 (0.9)" "42 (0.1)"
## TRANSIENT NON-COMMUNITY Total
## At-Risk "0 (0)" "658"
## Failing "0 (0)" "425"
## Not Assessed "12 (0.07)" "177"
## Not At-Risk "0 (0)" "1513"
## Potentially At-Risk "0 (0)" "431"counts <- table(safer$FEDERAL_CLASSIFICATION_TYPE, safer$FINAL_SAFER_STATUS)
bp <- barplot(counts,
col = c("steelblue", "orange", "lightgray"),
legend = TRUE,
args.legend = list(x = "topleft", bty = "n"))
# add labels to the bars
for (i in seq_len(ncol(counts))) {
ypos <- cumsum(counts[, i]) - counts[, i] / 2
labels <- ifelse(counts[, i] > 0, counts[, i], "")
text(bp[i], ypos, labels, col = "white", cex = 0.9)
}
Community systems make up 87-90% of every SAFER category, and NTNC systems hover around 10-13%. The composition looks the same whether the system is failing or safe. Transient Non-Community systems are effectively absent from the risk analysis.
TNC systems only show up in the “Not Assessed” row. Maybe they aren’t being evaluated in the SAFER framework.
Who Is Running the Failing Systems?
From the earlier analysis we know community systems make up 87-90% of every SAFER category. That means the risk story lives mostly inside community systems. Then I wanted to know if certain ownership types are more associated with failure.
- Who runs these community systems (private, public, etc.)?
# count owner types among community systems only
community <- table(safer$OWNER_TYPE[safer$FEDERAL_CLASSIFICATION_TYPE == "COMMUNITY"])
community##
## FEDERAL GOVERNMENT LOCAL MIXED (PUBLIC/PRIVATE)
## 35 1004 2
## PRIVATE STATE GOVERNMENT
## 1725 47prop_community <- round(prop.table(community), 2)
prop_community##
## FEDERAL GOVERNMENT LOCAL MIXED (PUBLIC/PRIVATE)
## 0.01 0.36 0.00
## PRIVATE STATE GOVERNMENT
## 0.61 0.0261% of community water systems in California are privately owned, with local government second at 36%. So I wanted to see whether private systems also fail at a disproportionate rate.
- Do some owner types fail or hit “At-Risk” status more than others?
# subset to community systems only
comm <- safer[safer$FEDERAL_CLASSIFICATION_TYPE == "COMMUNITY", ]
# cross-tab
tab <- table(comm$OWNER_TYPE, comm$FINAL_SAFER_STATUS)
tab##
## At-Risk Failing Not Assessed Not At-Risk
## FEDERAL GOVERNMENT 8 0 0 21
## LOCAL 168 96 116 507
## MIXED (PUBLIC/PRIVATE) 0 0 0 2
## PRIVATE 377 282 24 785
## STATE GOVERNMENT 18 1 3 16
##
## Potentially At-Risk
## FEDERAL GOVERNMENT 6
## LOCAL 117
## MIXED (PUBLIC/PRIVATE) 0
## PRIVATE 257
## STATE GOVERNMENT 9Of the 379 failing community systems, 282 (74%) are privately owned
Where Are the Failures Concentrated?
After I know it’s mostly private community systems failing, I wanted to see where geographically. Both by number of failing systems per county, and by how many people are actually affected.
Total Population Served by Failing or At-Risk Systems
# filter
subset_data <- safer[safer$FEDERAL_CLASSIFICATION_TYPE == "COMMUNITY" &
safer$OWNER_TYPE == "PRIVATE" &
safer$FINAL_SAFER_STATUS == "Failing", ]
# total population served by these systems
sum(subset_data$POPULATION, na.rm = TRUE)## [1] 138211# count by county, sorted high to low
sort_county <- sort(table(subset_data$COUNTY), decreasing = TRUE)
sort_county##
## KERN SONOMA SAN DIEGO MONTEREY SAN BERNARDINO
## 40 26 21 18 18
## MADERA LOS ANGELES FRESNO TULARE SAN LUIS OBISPO
## 15 14 9 9 8
## STANISLAUS VENTURA SAN JOAQUIN SANTA CLARA RIVERSIDE
## 8 8 7 7 6
## MENDOCINO SANTA BARBARA TUOLUMNE SAN BENITO SAN MATEO
## 5 5 5 4 4
## SANTA CRUZ CONTRA COSTA HUMBOLDT IMPERIAL INYO
## 4 3 3 3 3
## SACRAMENTO AMADOR BUTTE DEL NORTE LAKE
## 3 2 2 2 2
## SHASTA SUTTER TEHAMA TRINITY YUBA
## 2 2 2 2 2
## EL DORADO LASSEN MARIPOSA MERCED MONO
## 1 1 1 1 1
## NAPA NEVADA YOLO
## 1 1 1# bar chart: most-failing counties on top
sort_county <- rev(sort_county)
par(mar = c(4, 10, 4, 2))
bp <- barplot(sort_county,
horiz = TRUE,
las = 1,
cex.names = 0.2,
col = "steelblue")
Kern County has by far the most failing private community systems (40),
followed by Sonoma (26) and San Diego (21). But system count alone
doesn’t tell the full story. A county with fewer systems could still
affect far more people. So I also looked at population impact.
Counties with the Largest Affected Populations
# all failing or at-risk systems
risk_data <- safer[safer$FINAL_SAFER_STATUS %in% c("Failing", "At-Risk"), ]
# total population served
sum(risk_data$POPULATION, na.rm = TRUE)## [1] 2527963# break down by system type and status
aggregate(POPULATION ~ FEDERAL_CLASSIFICATION_TYPE + FINAL_SAFER_STATUS,
data = risk_data, FUN = sum, na.rm = TRUE)## FEDERAL_CLASSIFICATION_TYPE FINAL_SAFER_STATUS POPULATION
## 1 COMMUNITY At-Risk 1941950
## 2 NON-TRANSIENT NON-COMMUNITY At-Risk 26096
## 3 COMMUNITY Failing 543769
## 4 NON-TRANSIENT NON-COMMUNITY Failing 16148# population by county, failing or at-risk
county_pop <- aggregate(POPULATION ~ COUNTY, data = risk_data,
FUN = sum, na.rm = TRUE)
county_pop[order(-county_pop$POPULATION), ]## COUNTY POPULATION
## 19 LOS ANGELES 382287
## 36 SAN BERNARDINO 328986
## 15 KERN 240519
## 10 FRESNO 203912
## 33 RIVERSIDE 145607
## 41 SANTA BARBARA 121355
## 38 SAN JOAQUIN 100822
## 53 TULARE 92873
## 24 MERCED 88365
## 55 VENTURA 87307
## 13 IMPERIAL 87107
## 16 KINGS 82656
## 1 ALAMEDA 76739
## 49 STANISLAUS 74148
## 56 YOLO 56033
## 37 SAN DIEGO 50937
## 35 SAN BENITO 28930
## 28 NAPA 28920
## 40 SAN MATEO 26596
## 48 SONOMA 23164
## 27 MONTEREY 22783
## 29 NEVADA 20252
## 20 MADERA 19780
## 34 SACRAMENTO 15420
## 39 SAN LUIS OBISPO 11993
## 50 SUTTER 10320
## 23 MENDOCINO 8850
## 57 YUBA 7486
## 6 COLUSA 7209
## 26 MONO 6717
## 12 HUMBOLDT 5717
## 17 LAKE 5088
## 52 TRINITY 5080
## 54 TUOLUMNE 4877
## 5 CALAVERAS 4611
## 43 SANTA CRUZ 4187
## 21 MARIN 4079
## 3 AMADOR 4032
## 22 MARIPOSA 3900
## 44 SHASTA 3747
## 51 TEHAMA 3445
## 46 SISKIYOU 2914
## 18 LASSEN 2681
## 7 CONTRA COSTA 2460
## 32 PLUMAS 2144
## 14 INYO 1974
## 4 BUTTE 1834
## 42 SANTA CLARA 1759
## 45 SIERRA 1108
## 9 EL DORADO 910
## 8 DEL NORTE 645
## 25 MODOC 629
## 2 ALPINE 625
## 11 GLENN 500
## 30 ORANGE 372
## 47 SOLANO 322
## 31 PLACER 250# narrowing to failing only
failing_data <- safer[safer$FINAL_SAFER_STATUS == "Failing", ]
county_pop2 <- aggregate(POPULATION ~ COUNTY, data = failing_data,
FUN = sum, na.rm = TRUE)
county_pop2[order(-county_pop2$POPULATION), ]## COUNTY POPULATION
## 1 ALAMEDA 76689
## 8 FRESNO 75398
## 12 KERN 55845
## 44 TULARE 45413
## 13 KINGS 42123
## 20 MERCED 41099
## 26 RIVERSIDE 37710
## 16 LOS ANGELES 34540
## 30 SAN DIEGO 34071
## 29 SAN BERNARDINO 30194
## 40 STANISLAUS 11807
## 17 MADERA 8233
## 4 COLUSA 6959
## 23 NAPA 5530
## 39 SONOMA 5407
## 22 MONTEREY 5213
## 43 TRINITY 4300
## 32 SAN LUIS OBISPO 4032
## 10 IMPERIAL 3446
## 9 HUMBOLDT 3346
## 33 SAN MATEO 3186
## 46 VENTURA 3128
## 21 MONO 2894
## 14 LAKE 2449
## 31 SAN JOAQUIN 2095
## 5 CONTRA COSTA 1576
## 36 SANTA CRUZ 1576
## 25 PLUMAS 1514
## 28 SAN BENITO 1433
## 35 SANTA CLARA 1300
## 3 BUTTE 1187
## 34 SANTA BARBARA 1017
## 45 TUOLUMNE 1017
## 19 MENDOCINO 808
## 27 SACRAMENTO 420
## 48 YUBA 420
## 38 SISKIYOU 390
## 47 YOLO 360
## 42 TEHAMA 282
## 41 SUTTER 280
## 37 SHASTA 221
## 11 INYO 181
## 2 AMADOR 180
## 7 EL DORADO 150
## 15 LASSEN 140
## 18 MARIPOSA 130
## 6 DEL NORTE 118
## 24 NEVADA 110Counties with the most people served by failing systems are not necessarily the ones with the most failing systems. Alameda County ranks first in affected population, but when I went back to check, there is actually only one failing water system there. It just serves a very large number of people. This is a good reminder that system count and population impact tell different stories.
Fresno and Kern follow, and unlike Alameda, their high numbers reflect a concentration of many smaller failing systems.
So, this lead me to think whether system size matters. Maybe smaller water systems will more likely to fail compared to larger ones. If so, what size is most likely to fail.
System size
# convert to numeric first
safer <- safer %>%
mutate(WATER_QUALITY_SCORE = as.numeric(WATER_QUALITY_SCORE))## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `WATER_QUALITY_SCORE = as.numeric(WATER_QUALITY_SCORE)`.
## Caused by warning:
## ! NAs introduced by coercion# basic summary
summary(safer$WATER_QUALITY_SCORE)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0000 0.0000 0.0800 0.2738 0.5000 2.0000 178sum(safer$WATER_QUALITY_SCORE == 0, na.rm = TRUE)## [1] 1380sum(safer$WATER_QUALITY_SCORE > 0, na.rm = TRUE)## [1] 1646The dataset contains 3,026 systems with a water quality score. 1,380 of them (46%) score exactly 0, meaning they have no recorded water quality violations at all.
# histogram
ggplot(safer, aes(x = WATER_QUALITY_SCORE)) +
geom_histogram(binwidth = 0.05, fill = "#2c7bb6", color = "white") +
labs(
title = "Distribution of Water Quality Score",
x = "Water Quality Score",
y = "Count"
) +
theme_minimal()## Warning: Removed 178 rows containing non-finite outside the scale range
## (`stat_bin()`).
# log-transformed, easier to see the spread
hist(log10(safer$WATER_QUALITY_SCORE))
hist(log10(safer$WATER_QUALITY_SCORE+1))
there isn’t a smooth spectrum of water quality problems. Systems tend to
either have minor issues or significantly worse ones.
System Size Distribution
hist(safer$SERVICE_CONNECTIONS)
hist(log10(safer$SERVICE_CONNECTIONS))
safer %>%
filter(FINAL_SAFER_STATUS %in% c("Failing", "At-Risk", "Potentially At-Risk", "Not At-Risk")) %>%
ggplot(aes(x = log10(SERVICE_CONNECTIONS), fill = FINAL_SAFER_STATUS)) +
geom_density(alpha = 0.2) +
scale_fill_manual(values = c(
"Failing" = "red",
"At-Risk" = "yellow",
"Potentially At-Risk" = "blue",
"Not At-Risk" = "green"
))
A raw histogram of SERVICE_CONNECTIONS was difficult to read: most
systems are very small while a handful are extremely large, so small
systems get visually crushed. So I use a log10 transformation fixed
this.
The distribution peaked around log10 ≈ 1.5, meaning the most common water system in California serves only about 30 connections. These are very small, often serving a single neighborhood or rural community.
Then I overlaid density curves for all four risk categories. The peak at log10 ≈ 1.5 is tallest for Failing systems and progressively shorter as risk level decreases, meaning small systems are far more likely to be failing.
Who Bears the Risk most
correlation
# change to numeric
safer <- safer %>%
mutate(
MHI = as.numeric(MHI),
SERVICE_CONNECTIONS = as.numeric(SERVICE_CONNECTIONS),
WATER_QUALITY_SCORE = as.numeric(WATER_QUALITY_SCORE)
)## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `MHI = as.numeric(MHI)`.
## Caused by warning:
## ! NAs introduced by coercion# using log10(score + 1) to retain zero-score systems
# correlation
cor(safer$MHI, log10(safer$WATER_QUALITY_SCORE+1), use="pairwise.complete.obs")## [1] -0.1093724cor(safer$SERVICE_CONNECTIONS, log10(safer$WATER_QUALITY_SCORE+1), use="pairwise.complete.obs")## [1] -0.01915681Both correlations are weak…… So I decided to give up on this angel… No!! I wasn’t ready to drop the income angle entirely.
A low linear correlation doesn’t rule out a distributional pattern. So I looked at SERVICE_AREA_ECONOMIC_STATUS instead, which classifies communities into DAC tiers rather than treating income as a continuous variable.
safer$risk_flag <- ifelse(safer$FINAL_SAFER_STATUS %in% c("Failing", "At-Risk"),
"Failing/At-Risk", "Other")
table(safer$SERVICE_AREA_ECONOMIC_STATUS)##
## DAC Missing Non-DAC SDAC
## 12 938 3 1062 1189sum(safer$SERVICE_AREA_ECONOMIC_STATUS == "")## [1] 12#clean the blank data, put into "missing" category
safer$SERVICE_AREA_ECONOMIC_STATUS[safer$SERVICE_AREA_ECONOMIC_STATUS == ""] <- "Missing"
table(safer$SERVICE_AREA_ECONOMIC_STATUS)##
## DAC Missing Non-DAC SDAC
## 938 15 1062 1189#write csv for visualization in datawrapper
economic_status <- table(safer$risk_flag, safer$SERVICE_AREA_ECONOMIC_STATUS)
write.csv(economic_status, "economic_status.csv")Nearly half (49%) of Failing or At-Risk water systems serve severely disadvantaged communities (SDAC), compared to just 31% among systems that are not struggling. Meanwhile, systems serving wealthier Non-DAC communities make up 39% of the “Other” group but only 21% of Failing/At-Risk systems.
So, the poorer the community, the more likely its water system is in trouble.
MHI Comparison
Communities served by Failing or At-Risk systems have an average household income of $62,426, nearly $14,640 less than the $77,066 average for systems that are doing fine.
I ran a t-test to confirm this gap is statistically significant. The p-value is essentially zero, meaning the difference is almost certainly not due to chance.
safer$MHI <- as.numeric(safer$MHI)
aggregate(MHI ~ risk_flag, data = safer,
FUN = function(x) c(mean = mean(x, na.rm = TRUE),
median = median(x, na.rm = TRUE),
n = sum(!is.na(x))))## risk_flag MHI.mean MHI.median MHI.n
## 1 Failing/At-Risk 62425.69 58710.50 1078.00
## 2 Other 77065.84 70910.00 2111.00t.test(MHI ~ risk_flag, data = safer)##
## Welch Two Sample t-test
##
## data: MHI by risk_flag
## t = -12.022, df = 2516.1, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Failing/At-Risk and group Other is not equal to 0
## 95 percent confidence interval:
## -17028.14 -12252.15
## sample estimates:
## mean in group Failing/At-Risk mean in group Other
## 62425.69 77065.84Next, I look into Central Valley counties (Kern, Fresno, Merced, Tulare, and Madera), which I mentioned before that suffer the highest failing/at risk water systems in California. Failing or At-Risk systems serve communities earning on average $52,391, compared to $60,415 for systems that are not struggling. The gap is about $8,000.
central_valley <- safer %>%
filter(COUNTY %in% c("KERN", "FRESNO", "MERCED", "TULARE", "MADERA"))
t.test(MHI ~ risk_flag, data = central_valley)##
## Welch Two Sample t-test
##
## data: MHI by risk_flag
## t = -4.2331, df = 501.71, p-value = 2.742e-05
## alternative hypothesis: true difference in means between group Failing/At-Risk and group Other is not equal to 0
## 95 percent confidence interval:
## -11749.24 -4300.19
## sample estimates:
## mean in group Failing/At-Risk mean in group Other
## 52390.77 60415.49aggregate(MHI ~ risk_flag, data = central_valley,
FUN = function(x) c(mean = mean(x, na.rm = TRUE),
median = median(x, na.rm = TRUE),
n = sum(!is.na(x))))## risk_flag MHI.mean MHI.median MHI.n
## 1 Failing/At-Risk 52390.77 52299.00 313.00
## 2 Other 60415.49 59118.00 236.00What chemical drives most failing water systems?
We know that poorer communities are more exposed to failing water systems. The next question is: what chemical drives most failing water systems?
I started with the full statewide picture, then drilled down by county.
fail_all <- safer[safer$FINAL_SAFER_STATUS == "Failing" &
safer$PRIMARY_MCL_VIOLATION == "YES",]
analytes_all <- unlist(strsplit(fail_all$PRIMARY_ANALYTES, ";\\s*"))
analytes_all <- trimws(analytes_all)
violation_all <- sort(table(analytes_all), decreasing = TRUE)
violation_all## analytes_all
## NITRATE ARSENIC
## 86 76
## 1,2,3-TRICHLOROPROPANE COMBINED URANIUM
## 65 32
## TOTAL HALOACETIC ACIDS (HAA5) TTHM
## 25 23
## FLUORIDE GROSS ALPHA PARTICLE ACTIVITY
## 14 8
## NITRATE-NITRITE SELENIUM
## 6 6
## PERCHLORATE 1,2-DIBROMO-3-CHLOROPROPANE
## 3 2
## 1,1-DICHLOROETHYLENE ALUMINUM, DISSOLVED
## 1 1
## ANTIMONY, TOTAL CHROMIUM, HEX
## 1 1length(unique(analytes_all))## [1] 16#write csv for visualization in datawrapper
write.csv(violation_all, "violation_all.csv")Among California’s Failing water systems with a Primary MCL violation, Nitrate, Arsenic, and TCP are the contaminants most frequently exceeding safe drinking water standard
county_analytes_long <- safer %>%
filter(FINAL_SAFER_STATUS == "Failing" & PRIMARY_MCL_VIOLATION == "YES") %>%
select(COUNTY, PRIMARY_ANALYTES) %>%
mutate(analyte = strsplit(PRIMARY_ANALYTES, ";\\s*")) %>%
unnest(analyte) %>%
mutate(analyte = trimws(analyte)) %>%
count(COUNTY, analyte, name = "frequent") %>%
arrange(COUNTY, desc(frequent))
write_csv(county_analytes_long, "county_analytes_long.csv")Two patterns I want to mention:
First, 1,2,3-Trichloropropane (known as TCP) is heavily concentrated in Central Valley counties, particularly Kern, Fresno, Merced, and Tulare. TCP is a byproduct of pesticides that were banned decades ago, but because it binds to groundwater and breaks down extremely slowly, it continues to contaminate drinking water sources today.
Second, nitrate violations cluster in the same region, driven by decades of intensive agricultural fertilizer use and livestock waste that leach into groundwater.