<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Guided Walkthrough of Crash Modeling Tools
First, perform the necessary imports.
Step1:  Load in the initial data and see what predictors are available
Here we will look at crash data from I-90 in Washington State. The purpose of this demo is to show how the tool works, not to provide specific model results with deep meaning.
Step2:  Observe the summary statistics for the variables in our dataset
Here we consider two types of variables, namely continuous (cont) and categorical (cat). Below, we provide descriptions of each variable as well as its type. 
Note
Step3:  Show summary stats for the continuous variables
Step4:  Develop a crash frequency model
Here, we develop a model to predict crash frequency as a function of several characteristics of the roadway (here, we consider the annual average daily traffic, lane width, average grade, the presence of horizontal curvature, and surface type (e.g., asphalt concrete, portland cement concrete etc.) in our regression model). The model is a negative binomial (NB) regression model as is often the standard in highway safety modeling. The NB model is used to model count data and is preferrable to the Poisson model as it does not assume the mean and variance of the data are equal and is able to account for such overdispersion (variance>mean) by allowing the Poisson parameter to itself follow a gamma distribution.  It is critical we include an offset term since the segments are not the same length.
Mean 3-year crash count in dataset = 2.46
Varance of 3-year crash count in dataset = 26.41 
Variance-to-mean ratio = 10.72 --> Overdispersion
First compute the offset term
The offset term equals the natural logarithm of the segment length times the duration of the study period (here, we consider 3 years
Step5:  Now estimate the regression model and output the summary
Step6:  We will now compute the value of the safety performance function (spf) on a data set also from I-90
The spf simply involves computing the value of the regression model for a set of sites (i.e., road segments). The results returned are the predicted crash counts for the segments as estimated via the previously developed NB regression model.
Step7:  Implementation of the Emprical Bayes Method
Now, we will apply the Emprical Bayes (EB) method to computed the expected safety for each road segment. The EB method combines the expected crash count for each site as computed from the NB regression model and the observed crash count at each site in a weighted sum, where the weight factor is based on the dispersion parameter associated with the NB model (i.e., 1/scale parameter).
The EB method has many benefits including
Step8:  Calculation of Accident Reduction Potential
Accident Reduction Potential (ARP) is a commonly used method to rank/prioritize sites for safety improvement/treatment implementation. If $K_i$ (i.e., the observed crash count) is greater than $SPF_i$ (i.e., the predicted crash count), the ARP increases in value (compared to the converse) and the given site will be a high priority for treatment. The closer the value of $w_i$ is to 1, the higher the variance associated with the SPF, thus corresponding to lower ARP (low priority for treatment).
ARP is calculated as follows
Step9:  Confidence and Prediction Intervals for Negative Binomial Regression Model
The negative binomial model (and all other mixed-Poisson) regression models arise when an error term following a given distibution (known as the mixture distribution) enters the functional form of the Poisson parameter (denoted m). In the case of the NB model, the error term is gamma distributed which further implies that m itself follows a gamma distribution. Ultimately, there are three levels of variation in a mixed-Poisson model. Below, we describe these levels with respect to the NB model. WE then calculate and plot formulas for the confidence intervals (CIs) for $\mu$ and the prediction intervals (PIs) for m and y, as outlined in Wood (2005). Plotting of CIs and PIs gives a more complete picture for crash modeling as the user is presented with a range of values instead of just a point estimate as is obtained directly from a regression model.
Three Levels of Hierarchy in a Mixed-Poisson Model
Step10:  Calculate the values of the Poisson mean
Step11:  Calculate confidence intervals the Poisson means (mu)
Step12:  Calculate prediction intervals for the Poisson parameters (m)
LB = lower bound of interval, UB = upper bound of interval
Step13:  Calculate prediction intervals for the predicted responses (y)
Step14:  Plotting of the Intervals
Below, we generate a plot of the aforementioned confidence and prediction intervals (CIs and PIs, respectively) associated with the NB (aka Poisson-gamma) regression model.
<ASSISTANT_TASK:>
Python Code:
from crash_modeling_tools import *
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
Explanation: Guided Walkthrough of Crash Modeling Tools
First, perform the necessary imports.
End of explanation
crash_data = pd.read_csv('../data/crash_modeling_tools_demo_data/crash_data_final_90.csv')
crash_data = crash_data.dropna()
crash_data.head()
Explanation: Load in the initial data and see what predictors are available
Here we will look at crash data from I-90 in Washington State. The purpose of this demo is to show how the tool works, not to provide specific model results with deep meaning.
End of explanation
show_summary_stats(crash_data,[0,9])
Explanation: Observe the summary statistics for the variables in our dataset
Here we consider two types of variables, namely continuous (cont) and categorical (cat). Below, we provide descriptions of each variable as well as its type. 
Note: Not all variables will be used in the modeling.
surf_typ (cat): surface type
road_inv (n/a): road inventory number (i.e., route number), not a variable
lanewid (cont): lane width in feet
no_lanes (cont): number of lanes
lshldwid (cont): left shoulder width in feet
rshldwid (cont): right shoulder width in feet
medwid (cont): median width in feet
seg_lng (cont): segment length in feet
avg_grad (cont): average grade of the road segment
curve (cat): indicator for presence of curvature on the road segment
aadt_09 (cont): annual average daily traffic count for road segment in 2009
acc_ct_09 (cont): accident count for road segment in 2009
aadt_10 (cont): annual average daily traffic count for road segment in 2010
acc_ct_10 (cont): accident count for road segment in 2010
aadt_1 (cont): annual average daily traffic count for road segment in 2011
acc_ct_11 (cont): accident count for road segment in 2011
avg_addt (cont): average value of annual average daily traffic count for road segment from 2009-2011
tot_acc_count (cont): total accident count for road segment from 2009-2011
log_avg_aadt(cont): natural logarithm of annual average daily traffic count from 2009-2011
Show summary stats for the categorical variables we will consider
End of explanation
show_summary_stats(crash_data)
Explanation: Show summary stats for the continuous variables
End of explanation
offset_term = np.log(crash_data['seg_lng'] * 3)
Explanation: Develop a crash frequency model
Here, we develop a model to predict crash frequency as a function of several characteristics of the roadway (here, we consider the annual average daily traffic, lane width, average grade, the presence of horizontal curvature, and surface type (e.g., asphalt concrete, portland cement concrete etc.) in our regression model). The model is a negative binomial (NB) regression model as is often the standard in highway safety modeling. The NB model is used to model count data and is preferrable to the Poisson model as it does not assume the mean and variance of the data are equal and is able to account for such overdispersion (variance>mean) by allowing the Poisson parameter to itself follow a gamma distribution.  It is critical we include an offset term since the segments are not the same length.
Mean 3-year crash count in dataset = 2.46
Varance of 3-year crash count in dataset = 26.41 
Variance-to-mean ratio = 10.72 --> Overdispersion
First compute the offset term
The offset term equals the natural logarithm of the segment length times the duration of the study period (here, we consider 3 years: 2009, 2010, and 2011)
End of explanation
# need to cast log_avg_aadt to float since statsmodels thinks (incorrectly) that it is categorical
crash_data['log_aadt'] = crash_data.log_avg_aadt.astype(np.float)
mod_nb = smf.glm('tot_acc_ct~log_aadt+lanewid+avg_grad+C(curve)+C(surf_typ)',data=crash_data,offset=offset_term,
                 family=sm.families.NegativeBinomial()).fit()
mod_nb.summary()
Explanation: Now estimate the regression model and output the summary
End of explanation
data_eb = pd.read_csv('../data/crash_modeling_tools_demo_data/crash_data_eb.csv')
data_eb = data_eb.dropna()
data_eb['log_aadt'] = data_eb.log_avg_aadt.astype(np.float)
compute_spf(mod_nb,data_eb)
Explanation: We will now compute the value of the safety performance function (spf) on a data set also from I-90
The spf simply involves computing the value of the regression model for a set of sites (i.e., road segments). The results returned are the predicted crash counts for the segments as estimated via the previously developed NB regression model.
End of explanation
eb_safety_estimates = estimate_empirical_bayes(mod_nb,data_eb,data_eb['seg_lng'],data_eb['tot_acc_ct'])
eb_safety_estimates.head()
Explanation: Implementation of the Emprical Bayes Method
Now, we will apply the Emprical Bayes (EB) method to computed the expected safety for each road segment. The EB method combines the expected crash count for each site as computed from the NB regression model and the observed crash count at each site in a weighted sum, where the weight factor is based on the dispersion parameter associated with the NB model (i.e., 1/scale parameter).
The EB method has many benefits including:
- Easier to implement than full Bayes method
- Takes into account variance of the SPF (i.e., regression model) estimate
- Provides more robust estimate of site safety by considering observed crash count in conjunction with predicted crash count from similar sites
The EB method is formulated as follows:
$${\pi}_i = w_iSPF_i + (1-w_i)K_i$$
Where,
- ${\pi}_i$ = expected safety for site i
- $w_i$ = weight factor for site i
- $SPF_i$ = safety performance function for site i
- $K_i$ = observed crash count for site i
The weight factor $w_i$ is calculated as follows:
$$w_i = 1 / (1 + SPF_i/(kL_i))$$
Where,
- $k = 1/(NB scale parameter)$
- $L_i$ = segment length
- All other variables as aforementioned
Below, we calculate the EB esimates for each site and display the first few rows (not ranked in any specific order):
End of explanation
arp = calc_accid_reduc_potential(mod_nb,data_eb,data_eb['seg_lng'],data_eb['tot_acc_ct'])
arp.head()
Explanation: Calculation of Accident Reduction Potential
Accident Reduction Potential (ARP) is a commonly used method to rank/prioritize sites for safety improvement/treatment implementation. If $K_i$ (i.e., the observed crash count) is greater than $SPF_i$ (i.e., the predicted crash count), the ARP increases in value (compared to the converse) and the given site will be a high priority for treatment. The closer the value of $w_i$ is to 1, the higher the variance associated with the SPF, thus corresponding to lower ARP (low priority for treatment).
ARP is calculated as follows:
$$ARP_i = (1-w_i)(K_i-SPF_i)$$
Where,
- $ARP_i$ = accident reduction potential for site i
- $w_i$ = weight factor for site i
- $SPF_i$ = safety performance function for site i
- $K_i$ = observed crash count for site i
Below, we calculate ARP and show a few values (not in any particular order):
End of explanation
data_design = pd.read_csv('../data/crash_modeling_tools_demo_data/data_design.csv')
var_eta_hat = calc_var_eta_hat(mod_nb,data_design)
Explanation: Confidence and Prediction Intervals for Negative Binomial Regression Model
The negative binomial model (and all other mixed-Poisson) regression models arise when an error term following a given distibution (known as the mixture distribution) enters the functional form of the Poisson parameter (denoted m). In the case of the NB model, the error term is gamma distributed which further implies that m itself follows a gamma distribution. Ultimately, there are three levels of variation in a mixed-Poisson model. Below, we describe these levels with respect to the NB model. WE then calculate and plot formulas for the confidence intervals (CIs) for $\mu$ and the prediction intervals (PIs) for m and y, as outlined in Wood (2005). Plotting of CIs and PIs gives a more complete picture for crash modeling as the user is presented with a range of values instead of just a point estimate as is obtained directly from a regression model.
Three Levels of Hierarchy in a Mixed-Poisson Model:
- Lowest level: mean response ($\mu$), also known as the Poisson mean, which itself follows a lognormal distribution (that can be approximated by a normal distribution)
- One level up: Poisson parameter (m), alternately known as the safety, which when conditioned on the Poisson mean ($\mu$), follows the mixture distribution (here, gamma)
- Finally: predicted response (y), i.e., the crash frequency at a new site which when conditioned on the Poisson parameter (m), follows a Poisson distribution
Interested users can find full specification of the formulas for the CIs and PIs in the source code or in the reference below.
Reference:
Wood, G. (2005). "Confidence and prediction intervals for generalised linear accident models" Accident Analysis & Prevention, 37(2), 267-273.
Load in the data and compute the variance of the linear predictor
End of explanation
mu_hat = calc_mu_hat_nb(mod_nb,data_design)
mu_hat
Explanation: Calculate the values of the Poisson mean
End of explanation
ci_mu_nb = calc_ci_mu_nb(mu_hat,var_eta_hat)
ci_mu_nb.head()
Explanation: Calculate confidence intervals the Poisson means (mu)
End of explanation
pi_m_nb = calc_pi_m_nb(mod_nb,mu_hat,var_eta_hat)
pi_m_nb.head() 
Explanation: Calculate prediction intervals for the Poisson parameters (m)
LB = lower bound of interval, UB = upper bound of interval
End of explanation
pi_y_nb = calc_pi_y_nb(mod_nb,mu_hat,var_eta_hat)
pi_y_nb.head()
Explanation: Calculate prediction intervals for the predicted responses (y)
End of explanation
% matplotlib inline
aadt_range = np.arange(9700,148400,100)
plot_and_save_nb_cis_and_pis(data_design,mod_nb,mu_hat,var_eta_hat,aadt_range,'AADT')
Explanation: Plotting of the Intervals
Below, we generate a plot of the aforementioned confidence and prediction intervals (CIs and PIs, respectively) associated with the NB (aka Poisson-gamma) regression model.
End of explanation
<END_TASK>