Stat 171 Final Report
Sumedha Chalasani, Swaraj Dash, Megan Nakajima, Rhea Verma
2026-03-15
Data Cleaning
# raw data
df <- read.csv("Stormwater_overflow_events.csv")
# data cleaning
df <- df %>%
drop_na() %>%
distinct()
df <- df %>%
dplyr::select(overflow_events, rain_mm, impervious_pct, combined_sewer, maintenance_overdue, upstream_construction, catchment_type, green_infra_index)Exploratory Data Analysis (EDA)
Summary Statistics of Response Variable:
## Range: 0 13##
## Mean: 2.616162##
## Variance: 4.392711Distribution: The response variable is count data, meaning a Poisson distribution would be a good fit. However, there is overdispersion present in the response variable, so it actually follows the Negative Binomial distribution. There are no rows of data with repeated or missing data, but rain_c was taken out as rain_mm was preferred as a more interpretable predictor.
# descriptive statistics for overflow_events
df %>%
summarise(
Mean = mean(df$overflow_events),
SD = sd(df$overflow_events),
Median = median(df$overflow_events),
IQR = IQR(df$overflow_events),
Min = min(df$overflow_events),
Max = max(df$overflow_events),
Variance = var(df$overflow_events)
)## Mean SD Median IQR Min Max Variance
## 1 2.616162 2.09588 2 3 0 13 4.392711# descriptive statistics for rain_mm
df %>%
summarise(
Mean = mean(df$rain_mm),
SD = sd(df$rain_mm),
Median = median(df$rain_mm),
IQR = IQR(df$rain_mm),
Min = min(df$rain_mm),
Max = max(df$rain_mm),
Variance = var(df$rain_mm)
)## Mean SD Median IQR Min Max Variance
## 1 27.9596 14.77619 24.8 19.4 3.3 87.6 218.3359# descriptive statistics for impervious_pct
df %>%
summarise(
Mean = mean(df$impervious_pct),
SD = sd(df$impervious_pct),
Median = median(df$impervious_pct),
IQR = IQR(df$impervious_pct),
Min = min(df$impervious_pct),
Max = max(df$impervious_pct),
Variance = var(df$impervious_pct)
)## Mean SD Median IQR Min Max Variance
## 1 40.63401 15.76052 40.3 21.3 5 89.9 248.3941# descriptive statistics for green_infra_index
df %>%
summarise(
Mean = mean(df$green_infra_index),
SD = sd(df$green_infra_index),
Median = median(df$green_infra_index),
IQR = IQR(df$green_infra_index),
Min = min(df$green_infra_index),
Max = max(df$green_infra_index),
Variance = var(df$green_infra_index)
)## Mean SD Median IQR Min Max Variance
## 1 49.76162 5.072855 49.6 7.1 36.6 64.2 25.73386# descriptive statistics for combined_sewer
prop.table(table(df$combined_sewer))##
## 0 1
## 0.5319865 0.4680135# descriptive statistics for maintenance_overdue
prop.table(table(df$maintenance_overdue))##
## 0 1
## 0.6498316 0.3501684# descriptive statistics for upstream_construction
prop.table(table(df$upstream_construction))##
## 0 1
## 0.7340067 0.2659933# descriptive statistics for catchment_types
prop.table(table(df$catchment_type))##
## Commercial Industrial Mixed Residential
## 0.2558923 0.1178451 0.2121212 0.4141414# find mean overflow for binary variables
df %>%
group_by(combined_sewer) %>%
summarise(mean_overflow = mean(overflow_events),
sd_overflow = sd(overflow_events),
n = n())## # A tibble: 2 × 4
## combined_sewer mean_overflow sd_overflow n
## <int> <dbl> <dbl> <int>
## 1 0 2.18 1.69 158
## 2 1 3.11 2.39 139df %>%
group_by(maintenance_overdue) %>%
summarise(mean_overflow = mean(overflow_events),
sd_overflow = sd(overflow_events),
n = n())## # A tibble: 2 × 4
## maintenance_overdue mean_overflow sd_overflow n
## <int> <dbl> <dbl> <int>
## 1 0 2.30 1.62 193
## 2 1 3.21 2.68 104df %>%
group_by(catchment_type) %>%
summarise(mean_overflow = mean(overflow_events),
sd_overflow = sd(overflow_events),
n = n())## # A tibble: 4 × 4
## catchment_type mean_overflow sd_overflow n
## <chr> <dbl> <dbl> <int>
## 1 Commercial 2.49 1.81 76
## 2 Industrial 2.97 2.51 35
## 3 Mixed 2.60 2.36 63
## 4 Residential 2.60 2.00 123# check correlations between numeric variables and overflow_events
cor(df$overflow_events, df$rain_mm)## [1] 0.3490611cor(df$overflow_events, df$impervious_pct)## [1] 0.2265991cor(df$overflow_events, df$green_infra_index)## [1] 0.01818317model <- glm.nb(overflow_events ~ rain_mm + impervious_pct + green_infra_index + combined_sewer + maintenance_overdue + upstream_construction + catchment_type,
data = df)
summary(model)##
## Call:
## glm.nb(formula = overflow_events ~ rain_mm + impervious_pct +
## green_infra_index + combined_sewer + maintenance_overdue +
## upstream_construction + catchment_type, data = df, init.theta = 16.66146599,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.105244 0.409053 -0.257 0.79696
## rain_mm 0.017673 0.002428 7.278 3.38e-13 ***
## impervious_pct 0.011296 0.002554 4.424 9.70e-06 ***
## green_infra_index -0.003929 0.007817 -0.503 0.61520
## combined_sewer 0.364351 0.079989 4.555 5.24e-06 ***
## maintenance_overdue 0.259273 0.080690 3.213 0.00131 **
## upstream_construction -0.038045 0.088794 -0.428 0.66831
## catchment_typeIndustrial 0.202087 0.135149 1.495 0.13484
## catchment_typeMixed -0.068650 0.117650 -0.584 0.55955
## catchment_typeResidential -0.070436 0.100747 -0.699 0.48447
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(16.6615) family taken to be 1)
##
## Null deviance: 431.19 on 296 degrees of freedom
## Residual deviance: 323.10 on 287 degrees of freedom
## AIC: 1122.3
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 16.66
## Std. Err.: 9.41
##
## 2 x log-likelihood: -1100.285# check for multicollinearity
vif(model)## GVIF Df GVIF^(1/(2*Df))
## rain_mm 1.033424 1 1.016575
## impervious_pct 1.049283 1 1.024345
## green_infra_index 1.021748 1 1.010815
## combined_sewer 1.045263 1 1.022381
## maintenance_overdue 1.040516 1 1.020057
## upstream_construction 1.055161 1 1.027211
## catchment_type 1.091805 3 1.014746numeric_variables <- df[, c("overflow_events", "rain_mm", "impervious_pct", "green_infra_index")]
ggpairs(numeric_variables)
Significant predictors:
- rain_mm: Highly significantly positive (more rainfall -> more overflow events).
- impervious_pct: Significantly positive (more impervious land -> more overflow events).
- combined_sewer: Catchments with combined sewers have significantly higher expected overflow counts.
- maintenance_overdue: Overdue maintenance increases expected overflows.
Non-significant predictors:
- green_infra_index
- upstream_construction
- catchment_type
There appears to be more overflow events when there was higher amounts of rain. The amount of overflow events seems to be somewhat evenly distributed amongst the catchment types.
There appears to be more data for residential catchment type. The distribution of count of impervious percentage per each catchment type appears to be pretty similar.
There appears to be a moderately even distribution of green infra index scores per each catchment type, with the distribution for each one being pretty similar regardless of maintenance overdue.
There appears to be less upstream construction overall, with a higher amount of no upstream construction if it is not a combined sewer.
There appears to be a higher amount of events where it was not a combined sewer and there was no upstream construction. There is a fairly even amount of this combination for there being and not being maintenance overdue. There appears to be a higher amount of events overall if there was maintenance overdue. Data seems spread evenly when there is less rain, but there is a high count of events when rain increases and maintenance is overdue. Outliers when this is true also do have a combined sewer and upstream construction.
Most of the data is between overflow count of 0 to 7. There is pretty even distribution of events per impervious percentage, except for a few outliers when it is 50%. Most of the rain amounts are low but the high overflow counts correspond with heavy rain.
Methodology
First Order Model
model_full <- glm.nb(overflow_events ~ ., data = training_data)
summary(model_full)##
## Call:
## glm.nb(formula = overflow_events ~ ., data = training_data, init.theta = 17.18612833,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4118147 0.4888816 -0.842 0.399586
## rain_mm 0.0209775 0.0030541 6.869 6.48e-12 ***
## impervious_pct 0.0115499 0.0029782 3.878 0.000105 ***
## combined_sewer 0.3767579 0.0965302 3.903 9.50e-05 ***
## maintenance_overdue1 0.3020182 0.0974090 3.101 0.001932 **
## upstream_construction1 -0.0663958 0.1090332 -0.609 0.542557
## catchment_typeIndustrial 0.1223575 0.1590551 0.769 0.441729
## catchment_typeMixed -0.1797131 0.1405208 -1.279 0.200930
## catchment_typeResidential -0.1132117 0.1211230 -0.935 0.349951
## green_infra_index -0.0004957 0.0094968 -0.052 0.958371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(17.1861) family taken to be 1)
##
## Null deviance: 310.85 on 206 degrees of freedom
## Residual deviance: 220.27 on 197 degrees of freedom
## AIC: 777.17
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 17.2
## Std. Err.: 12.0
##
## 2 x log-likelihood: -755.168# Perform wald test
wald.test(Sigma = vcov(model_full), b = coef(model_full), Terms = 1:9)## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 87.6, df = 9, P(> X2) = 5e-15# Evaluates the null hypothesis that all of the first 9 coefficients are equal to zero
par(mfrow=c(2,2))
plot(model_full)
Model with Added Interactions
model_int <- glm.nb(overflow_events ~ rain_mm*impervious_pct + rain_mm*combined_sewer + maintenance_overdue + upstream_construction + catchment_type + green_infra_index, data = training_data)
summary(model_int)##
## Call:
## glm.nb(formula = overflow_events ~ rain_mm * impervious_pct +
## rain_mm * combined_sewer + maintenance_overdue + upstream_construction +
## catchment_type + green_infra_index, data = training_data,
## init.theta = 20.24909611, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.0369885 0.5795251 0.064 0.94911
## rain_mm 0.0071365 0.0100476 0.710 0.47754
## impervious_pct 0.0029633 0.0069937 0.424 0.67177
## combined_sewer 0.2896734 0.2179430 1.329 0.18381
## maintenance_overdue1 0.2947343 0.0964392 3.056 0.00224 **
## upstream_construction1 -0.0779627 0.1083289 -0.720 0.47172
## catchment_typeIndustrial 0.1296796 0.1576730 0.822 0.41082
## catchment_typeMixed -0.1670980 0.1395593 -1.197 0.23118
## catchment_typeResidential -0.0997509 0.1202935 -0.829 0.40697
## green_infra_index -0.0011955 0.0095920 -0.125 0.90081
## rain_mm:impervious_pct 0.0002854 0.0002104 1.356 0.17495
## rain_mm:combined_sewer 0.0026983 0.0060476 0.446 0.65547
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(20.2491) family taken to be 1)
##
## Null deviance: 316.80 on 206 degrees of freedom
## Residual deviance: 222.16 on 195 degrees of freedom
## AIC: 779.16
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 20.2
## Std. Err.: 16.5
##
## 2 x log-likelihood: -753.156par(mfrow=c(2,2))
plot(model_int)
Model With Quadratic + Interaction Terms
model_quad <- glm.nb(overflow_events ~ poly(rain_mm, 2) + rain_mm*impervious_pct + rain_mm*combined_sewer + maintenance_overdue + upstream_construction + catchment_type + green_infra_index, data = training_data)
summary(model_quad)##
## Call:
## glm.nb(formula = overflow_events ~ poly(rain_mm, 2) + rain_mm *
## impervious_pct + rain_mm * combined_sewer + maintenance_overdue +
## upstream_construction + catchment_type + green_infra_index,
## data = training_data, init.theta = 20.02502277, link = log)
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2465092 0.4946702 0.498 0.61825
## poly(rain_mm, 2)1 1.8082376 2.2027456 0.821 0.41170
## poly(rain_mm, 2)2 -0.4790037 0.6674549 -0.718 0.47297
## rain_mm NA NA NA NA
## impervious_pct 0.0032095 0.0071042 0.452 0.65143
## combined_sewer 0.3033036 0.2219330 1.367 0.17174
## maintenance_overdue1 0.2985620 0.0965706 3.092 0.00199 **
## upstream_construction1 -0.0774548 0.1083617 -0.715 0.47474
## catchment_typeIndustrial 0.1377592 0.1581391 0.871 0.38368
## catchment_typeMixed -0.1715280 0.1399186 -1.226 0.22023
## catchment_typeResidential -0.1023701 0.1205642 -0.849 0.39583
## green_infra_index -0.0012642 0.0095658 -0.132 0.89486
## rain_mm:impervious_pct 0.0002722 0.0002140 1.272 0.20339
## rain_mm:combined_sewer 0.0022081 0.0061901 0.357 0.72131
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(20.025) family taken to be 1)
##
## Null deviance: 316.42 on 206 degrees of freedom
## Residual deviance: 221.39 on 194 degrees of freedom
## AIC: 780.63
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 20.0
## Std. Err.: 16.0
##
## 2 x log-likelihood: -752.634par(mfrow=c(2,2))
plot(model_quad)
Model with only Interactions
model_int_only <- glm.nb(overflow_events ~ rain_mm*impervious_pct + rain_mm*combined_sewer, data = training_data)
summary(model_int_only)##
## Call:
## glm.nb(formula = overflow_events ~ rain_mm * impervious_pct +
## rain_mm * combined_sewer, data = training_data, init.theta = 12.51779913,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.0693922 0.3435850 0.202 0.840
## rain_mm 0.0043287 0.0103112 0.420 0.675
## impervious_pct 0.0028085 0.0071276 0.394 0.694
## combined_sewer 0.2661240 0.2223116 1.197 0.231
## rain_mm:impervious_pct 0.0003169 0.0002155 1.471 0.141
## rain_mm:combined_sewer 0.0034290 0.0061956 0.553 0.580
##
## (Dispersion parameter for Negative Binomial(12.5178) family taken to be 1)
##
## Null deviance: 297.32 on 206 degrees of freedom
## Residual deviance: 221.35 on 201 degrees of freedom
## AIC: 779.46
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 12.52
## Std. Err.: 6.91
##
## 2 x log-likelihood: -765.463par(mfrow=c(2,2))
plot(model_int_only)
Reduced Model
model_reduced <- glm.nb(overflow_events ~ rain_mm + impervious_pct + combined_sewer, data = training_data)
summary(model_reduced)##
## Call:
## glm.nb(formula = overflow_events ~ rain_mm + impervious_pct +
## combined_sewer, data = training_data, init.theta = 11.16512861,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.403296 0.183172 -2.202 0.027684 *
## rain_mm 0.019973 0.003118 6.406 1.50e-10 ***
## impervious_pct 0.012357 0.002995 4.126 3.69e-05 ***
## combined_sewer 0.377227 0.098550 3.828 0.000129 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(11.1651) family taken to be 1)
##
## Null deviance: 291.75 on 206 degrees of freedom
## Residual deviance: 219.85 on 203 degrees of freedom
## AIC: 777.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 11.17
## Std. Err.: 5.59
##
## 2 x log-likelihood: -767.901par(mfrow=c(2,2))
plot(model_reduced)
Results
Final Model
# The final test model is a variation of the reduced model with the addition of maintenance_overdue variable.
final_model <- glm.nb(
overflow_events ~ rain_mm + impervious_pct + combined_sewer +
maintenance_overdue, data = training_data )
summary(final_model)##
## Call:
## glm.nb(formula = overflow_events ~ rain_mm + impervious_pct +
## combined_sewer + maintenance_overdue, data = training_data,
## init.theta = 15.94943497, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.480086 0.180369 -2.662 0.007775 **
## rain_mm 0.019945 0.003009 6.629 3.38e-11 ***
## impervious_pct 0.011828 0.002910 4.064 4.82e-05 ***
## combined_sewer 0.349652 0.095937 3.645 0.000268 ***
## maintenance_overdue1 0.283574 0.095959 2.955 0.003125 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(15.9494) family taken to be 1)
##
## Null deviance: 307.9 on 206 degrees of freedom
## Residual deviance: 222.8 on 202 degrees of freedom
## AIC: 771.66
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 15.9
## Std. Err.: 10.7
##
## 2 x log-likelihood: -759.658par(mfrow=c(2,2))
plot(final_model)
AIC Test
# All models' AIC values are outputted for comparison
AIC(model_full, model_int, model_quad, model_int_only, model_reduced, final_model)## df AIC
## model_full 11 777.1676
## model_int 13 779.1564
## model_quad 14 780.6336
## model_int_only 7 779.4635
## model_reduced 5 777.9010
## final_model 6 771.6579Likelihood Ratio Tests
# LRT for Full Model vs Full Model w/ Interaction Terms
full_int <- lrtest(model_full, model_int)[2, "Pr(>Chisq)"]
cat("Full vs Interactions:", full_int)## Full vs Interactions: 0.3658393# LRT for Full Model vs Full Model w/ Quadratic & Interaction Terms
full_quad <- lrtest(model_full, model_quad)[2, "Pr(>Chisq)"]
cat("\nFull vs Quadratic + Interactions:", full_quad)##
## Full vs Quadratic + Interactions: 0.4691799# LRT for Full Model vs Interaction Terms Only Model
full_int_only <- lrtest(model_full, model_int_only)[2, "Pr(>Chisq)"]
cat("\nFull vs Interactions Only:", full_int_only)##
## Full vs Interactions Only: 0.0357273# LRT for Full Model vs Reduced Model (w/o maintenance_overdue variable)
full_reduced <- lrtest(model_full, model_reduced)[2, "$Pr(>Chisq)"]
cat("\nFull vs Reduced:", full_reduced)##
## Full vs Reduced:# LRT for Interaction Terms Only Model vs Reduced Model (w/o maintenance_overdue variable)
int_only_reduced <- lrtest(model_int_only, model_reduced)[2, "Pr(>Chisq)"]
cat("\nInteractions Only vs Reduced:", int_only_reduced)##
## Interactions Only vs Reduced: 0.2955961# LRT for Final Model (Reduced Model w/ maintenance_overdue) vs Reduced Model (w/o maintenance_overdue)
final_reduced <- lrtest(final_model, model_reduced)[2, "Pr(>Chisq)"]
cat("\nFinal vs Reduced:", final_reduced)##
## Final vs Reduced: 0.004090648Testing Model
test_model <- glm.nb(
overflow_events ~ rain_mm + impervious_pct + combined_sewer +
maintenance_overdue, data = testing_data)## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reachedsummary(test_model)##
## Call:
## glm.nb(formula = overflow_events ~ rain_mm + impervious_pct +
## combined_sewer + maintenance_overdue, data = testing_data,
## init.theta = 21443.66227, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.197909 0.262670 -0.753 0.45118
## rain_mm 0.017334 0.004347 3.987 6.68e-05 ***
## impervious_pct 0.011097 0.004691 2.366 0.01799 *
## combined_sewer 0.441172 0.135650 3.252 0.00114 **
## maintenance_overdue1 0.226136 0.138066 1.638 0.10144
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(21443.66) family taken to be 1)
##
## Null deviance: 137.95 on 88 degrees of freedom
## Residual deviance: 103.65 on 84 degrees of freedom
## AIC: 341.19
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 21444
## Std. Err.: 492504
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -329.193# Checking dispersion in training/testing data
cat("Training Dispersion:", var(training_data$overflow_events) / mean(training_data$overflow_events))## Training Dispersion: 1.810762cat("\nTesting Dispersion:", var(testing_data$overflow_events) / mean(testing_data$overflow_events))##
## Testing Dispersion: 1.413655# ratio based off training data
cat("\nTheta vs SE:", summary(final_model)$theta / summary(final_model)$SE)##
## Theta vs SE: 1.489214- Convergence error in model run on testing data is due to only slight overdispersion present. Testing data sample is also small, which contributed to the error. Model works well on training data.
Evaluating model
# obtaining predicted values
pred_test <- predict(final_model, newdata = testing_data)
# calculating residual mean squared error
RMSE <- sqrt(mean((testing_data$overflow_events - pred_test)^2))
cat("RMSE:", RMSE)## RMSE: 2.685286# obtaining predicted values
df_pred <- testing_data %>%
mutate(predicted = predict(test_model, type = "response"),
actual = overflow_events)
# plot to compare actual vs predicted overflow events
ggplot(df_pred, aes(x = actual, y = predicted)) +
geom_point() +
geom_abline(color = "red", linetype = "dashed") +
labs(title = "Predicted vs Actual Overflow Events",
x = "Predicted", y = "Observed",
caption = paste("RMSE =", round(RMSE, 2)))