Getting Started With Exploratory Data Analysis (EDA)

This notebook serves as a starter guide or template for exploratory data analysis. It will go over the topics mentioned in the EDA guide.

# let's start off by importing the libraries we will need for eda
import pandas as pd
import numpy as np 

# for visualizations : 
import seaborn as sns
import matplotlib.pyplot as plt

The dataset we will be using in this tutorial is from Analyze Boston. Analyze Boston is the City of Boston’s data hub and is a great resource for data sets regarding the city.

We will be working with the 2022 311 Service Requests dataset. The dataset consists of service requests from all channels of engagement. 311 allows you to report non-emergency issues or request non-emergency City services.

Link to dataset

# to run in colab, run the following lines 
# from google.colab import drive
# drive.mount('/content/drive')

Mounted at /content/drive

# read in dataset
df = pd.read_csv('311-requests.csv') 
pd.set_option('display.max_columns', 6)

# let's look at the first five rows of the dataset
df.head()

How many observations/rows are there?

How many variables/columns are there?

What kinds of variables are there? Qualitative? Quantitative? Both?

# number of observations 
df.shape[0]

146373

# to see column name, count, and dtype of each column 
df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 146373 entries, 0 to 146372
Data columns (total 29 columns):
 #   Column                          Non-Null Count   Dtype  
---  ------                          --------------   -----  
 0   case_enquiry_id                 146373 non-null  int64  
 1   open_dt                         146373 non-null  object 
 2   target_dt                       129475 non-null  object 
 3   closed_dt                       125848 non-null  object 
 4   ontime                          146373 non-null  object 
 5   case_status                     146373 non-null  object 
 6   closure_reason                  146373 non-null  object 
 7   case_title                      146371 non-null  object 
 8   subject                         146373 non-null  object 
 9   reason                          146373 non-null  object 
 10  type                            146373 non-null  object 
 11  queue                           146373 non-null  object 
 12  department                      146373 non-null  object 
 13  submittedphoto                  55740 non-null   object 
 14  closedphoto                     0 non-null       float64
 15  location                        146373 non-null  object 
 16  fire_district                   146149 non-null  object 
 17  pwd_district                    146306 non-null  object 
 18  city_council_district           146358 non-null  object 
 19  police_district                 146306 non-null  object 
 20  neighborhood                    146210 non-null  object 
 21  neighborhood_services_district  146358 non-null  object 
 22  ward                            146373 non-null  object 
 23  precinct                        146270 non-null  object 
 24  location_street_name            145030 non-null  object 
 25  location_zipcode                110519 non-null  float64
 26  latitude                        146373 non-null  float64
 27  longitude                       146373 non-null  float64
 28  source                          146373 non-null  object 
dtypes: float64(4), int64(1), object(24)
memory usage: 32.4+ MB

There are 146373 rows (observations).

There are 29 columns (variables).

There are both categorical and numerical variables. At quick glance there seems to be more categorical variables than numerical variables.

Categorical Variables: case_status, neighborhood, source, etc.

Numerical Variables: … maybe not?

The case_enquiry_id is a unique identifier for each row, closedphoto has 0 non-null values so it might be worth it to drop this column since there is no additional information we can gather, columns such as location_zipcode, latitude, longitude not exactly numeric varaibles, since they are numbers that represent different codes.

Cleaning

Let’s convert the three time variables (open_dt, target_dt, and closed_dt) from objects to pandas datetime objects. Let’s focus on service requests for a set period of time in 2022. We will start by filtering for service requests that were opened from January 2022 to March 2022.

# changing the three columns with dates and times to pandas datetime object 
df['open_dt'] = pd.to_datetime(df['open_dt'])
df['target_dt'] = pd.to_datetime(df['target_dt'])
df['closed_dt'] = pd.to_datetime(df['closed_dt'])

# output is long, but run the line below to check the type of the three columns 
#df.dtypes

# filter data for 311 requests from january 2022 to march 2022 
df_filtered = df.loc[(df['open_dt'] >= '2022-01-01') &
                  (df['open_dt'] < '2022-03-31')]
df_filtered.head()

From our previous observation, since closedphoto column does not contain any non-null values, let’s drop it.

# drop closedphoto column
df_filtered = df_filtered.drop(columns=['closedphoto'])

After filtering the service requests, let’s see how many observations we are left with.

# how many requests were opened from Jan 2022 to March 2022
df_filtered.shape[0]

66520

From a quick preview of the dataframe, we can see that some of the requests are still open. Let’s see how many observations are open vs. closed and then how many are ontime vs. overdue from the set of requests from January 2022 to March 2022.

# checking how many open vs. closed cases
df_filtered['case_status'].value_counts()

Closed    59420
Open       7100
Name: case_status, dtype: int64

# visualize case_status in pie chart, set color palette 
colors = sns.color_palette('muted')[0:5]
ax = df_filtered['case_status'].value_counts().plot.pie(colors=colors)
ax.yaxis.set_visible(False)

# checking how many ontime vs. overdue cases 
df_filtered['ontime'].value_counts()

ONTIME     55089
OVERDUE    11431
Name: ontime, dtype: int64

# visualize ontime in pie chart, set color palette 
colors = sns.color_palette('bright')[0:5]
ax = df_filtered['ontime'].value_counts().plot.pie(colors=colors)
ax.yaxis.set_visible(False)

Descriptive Statistics

Pandas makes this easy! We can use describe() to get the descriptive statistics of the numerical columns.

df_filtered.describe()

As mentioned before, the case_enquiry_id, location_zipcode, latitude, and longitude columns are not numeric variables. The descriptive statistics are not very useful in this situation.

What would be a useful numeric variable is the duration of a request. Let’s calculate the duration of each of the requests from January 2022 to March 2022 and add it as a new column in our dataframe.

# calculating case duration and adding a new column (case_duration) to the dataframe 
duration = df_filtered['closed_dt'] - df_filtered['open_dt']
df_filtered = df_filtered.assign(case_duration=duration)
df_filtered.head()

Now we can see the new case_duration column. Some values are NaT, which means there is a missing date. This makes sense because the case_status is OPEN.

Let’s filter out the open cases and focus on analyzing the duration of the closed cases.

# filter out the open cases
df_closed = df_filtered.loc[(df_filtered['case_status'] == "Closed")]
df_closed.head()

With the closed cases, let’s calculate the descriptive statistics of the new case_duration column.

# let's calculate the descriptive statistics again 
# using double brackets to display in a *fancy* table format
df_closed[['case_duration']].describe()

From the table, we can see that the average case duration is ~4.5 days.
The standard deviation for the case duration is ~15.4 days.
The minimum time a case takes to close is 4 minutes.
The maximum time a case takes to close is ~181.6 days.
The inter-quartile range (IQR) is the difference between the 25% and 75% quantiles.

We can also calculate the mode and median.

df_closed['case_duration'].mode()

0   0 days 00:00:54
1   0 days 00:00:57
2   0 days 00:01:03
Name: case_duration, dtype: timedelta64[ns]

df_closed['case_duration'].median()

Timedelta('0 days 09:01:45')

The descriptive statistics summary in table form is nice, but it would be nice to visualize the data in a histogram. Simply trying to plot using the values in the case_duration column will case an error.

Currently, the values in case_duration are of type timedelta64[ns], df_closed['case_duration'] is a Timedelta Series. We will need to apply what is called a frequency conversion to the values.

“Timedelta Series, TimedeltaIndex, and Timedelta scalars can be converted to other ‘frequences’ by dividing by another timedelta, or by astyping to a specific timedelta type.” (See the link below for more information and code examples!)

https://pandas.pydata.org/pandas-docs/stable/user_guide/timedeltas.html

# dividing the case_duration values by Timedelta of 1 day 
duration_days = ( df_closed['case_duration'] / pd.Timedelta(days=1))

# adding calculation to dataframe under duration_in_days column 
df_closed = df_closed.assign(duration_in_days=duration_days)

# display descriptive statistics summary with new column addition 
df_closed[['duration_in_days']].describe()

# using seaborn library for visualizations 
sns.set_theme() # use this if you dont want the visualizations to be default matplotlibstyle
sns.displot(df_closed, x="duration_in_days", binwidth=1)

<seaborn.axisgrid.FacetGrid at 0x16ae584c0>

From the plot above, the data seems to be skewed right meaning the right tail is much longer than the left. Let’s try playing with different bin widths.

# trying different bin sizes 
sns.displot(df_closed, x="duration_in_days", binwidth=5)

<seaborn.axisgrid.FacetGrid at 0x16b031790>

# trying different bin sizes
sns.displot(df_closed, x="duration_in_days", binwidth=25)

<seaborn.axisgrid.FacetGrid at 0x16b06b340>

Since the data is heavily skewed, let’s apply log transformation to the data. The log transformation will hopefully reduce or remove the skewness of the original data. The assumption is that the original data follows a log-normal distribution.

# log-scale transformation since the data is heavliy skewed 
# add bin_width parameter to change bin sizes
sns.displot(df_closed, x="duration_in_days", log_scale=True)

<seaborn.axisgrid.FacetGrid at 0x16b5577f0>

Which neighborhoods had the most requests from January 2022 - March 2022?

To answer this question, we will take a look at the neighborhood column.

# has 25 unique values so a pie chart probably is not the best option
len(df_closed['neighborhood'].unique()) 

25

# plot neighborhood counts 
sns.countplot(x="neighborhood", data=df_closed).set_title('Number of Requests by Neighborhood')

Text(0.5, 1.0, 'Number of Requests by Neighborhood')

Yikes! The x-axis labels are pretty hard to read. Let’s fix that by plotting the bars horizontally.

# fixing orientation of the labels 
sns.countplot(y="neighborhood", data=df_closed).set_title('Number of Requests by Neighborhood')

Text(0.5, 1.0, 'Number of Requests by Neighborhood')

From the plot we can see that Dorchester has the most requests, followed by South Boston/South Boston Waterfront, then Roxbury. There’s a bar that doesn’t have a name…that’s strange. Let’s display the exact counts for each neighborhood.

# displaying number of requests by neighborhood in table form 
df_closed['neighborhood'].value_counts()

Dorchester                                      9148
South Boston / South Boston Waterfront          5608
Roxbury                                         5097
East Boston                                     4420
Allston / Brighton                              3945
Jamaica Plain                                   3696
South End                                       3666
Downtown / Financial District                   3419
Back Bay                                        2740
Greater Mattapan                                2429
Hyde Park                                       2308
Charlestown                                     2096
Roslindale                                      2083
Boston                                          1803
West Roxbury                                    1698
Beacon Hill                                     1595
Fenway / Kenmore / Audubon Circle / Longwood    1034
Mission Hill                                     990
South Boston                                     600
                                                 476
Brighton                                         293
Allston                                          141
Mattapan                                          72
Chestnut Hill                                      4
Name: neighborhood, dtype: int64

There are 476 requests without a neighborhood label.

# uncomment and run the line below to check for the empty neighborhood label 
# print(df_closed['neighborhood'].unique())

# gather the rows where neighborhood == ' ' 
df_no_neighborhood = df_closed.loc[(df_closed['neighborhood'] == ' ')]
df_no_neighborhood.head(15) # display first 15 rows

print(df_no_neighborhood['latitude'].unique())
print(df_no_neighborhood['longitude'].unique())

[42.3594]
[-71.0587]

The latitude and longitude values are the same for all of the rows without a neighborhood value. We can use the Geopy module to convert the latitude and longitude coordinates to a place or location address - also referred to as reverse geocoding.

# import geopy 
from geopy.geocoders import Nominatim 

# make a Nominatim object and initialize, specify a user_agent 
# Nominatim requires this value to be set to your application name, to be able to limit the number of requests per application
# Nominatim is a free service but provides low request limits: https://operations.osmfoundation.org/policies/nominatim/
geolocator = Nominatim(user_agent="eda_geotest")

# set latitude and longitude and convert to string 
lat = str(df_no_neighborhood['latitude'].unique()[0])
long = str(df_no_neighborhood['longitude'].unique()[0])

# get the location information
location = geolocator.reverse(lat + "," +long)

# display location information, add .raw for more details
print(location.raw)

{'place_id': 264213803, 'licence': 'Data © OpenStreetMap contributors, ODbL 1.0. https://osm.org/copyright', 'osm_type': 'way', 'osm_id': 816277585, 'lat': '42.3594696', 'lon': '-71.05880376899256', 'display_name': "Sears' Crescent and Sears' Block Building, Franklin Avenue, Downtown Crossing, Downtown Boston, Boston, Suffolk County, Massachusetts, 02201, United States", 'address': {'building': "Sears' Crescent and Sears' Block Building", 'road': 'Franklin Avenue', 'neighbourhood': 'Downtown Crossing', 'suburb': 'Downtown Boston', 'city': 'Boston', 'county': 'Suffolk County', 'state': 'Massachusetts', 'ISO3166-2-lvl4': 'US-MA', 'postcode': '02201', 'country': 'United States', 'country_code': 'us'}, 'boundingbox': ['42.3593149', '42.3596061', '-71.0592779', '-71.0584887']}

Quick Google Maps search of the location confirms that (42.3594, -71.0587) is Government Center. The output from geopy is Sear’s Crescent and Sears’ Block which are a pair of buildings adjacent to City Hall and City Hall Plaza, Government Center.

Another quick look at the output from geopy shows that the lat and lon values are similar but different from the latitude and longitude values in the dataset.

The requests without a neighborhood value have a general location of Government Center. At least we can confirm that requests without a neighborhood value are not outside of Boston or erroneous.

During January 2022 - March 2022, where did the most case requests come from?

To answer this question, we will take a look at the source column.

# has only 5 unique values so in this case we can use a pie chart 
len(df_closed['source'].unique())

5

# displaying the number of requests by each source type 
df_closed['source'].value_counts()

Citizens Connect App    32066
Constituent Call        21051
City Worker App          3795
Self Service             1632
Employee Generated        876
Name: source, dtype: int64

# visualizing the breakdown of where case requests come from 
# seaborn doesn't have a default pie chart but you can add seaborn color palettes to matplotlib plots

colors = sns.color_palette('pastel')[0:5]
ax = df_closed['source'].value_counts().plot.pie(colors=colors)

# label each slice with the percentage of requests per source 
ax = df_closed['source'].value_counts().plot.pie(colors=colors,autopct='%1.1f%%')

# run the following to remove the default column name label *source*
#ax.yaxis.set_visible(False)

From the pie chart, 54% of the requests from January 2022 - March 2022 came from the Citizens Connect App, 35.4% came from a Constituent Call, followed by 6.4% from the City Worker App.

How many different types of requests were there from January 2022 - March 2022?

To answer this question, we will take a look at the reason column.

# how many different reasons are there 
len(df_closed['reason'].unique())

38

# number of requests by reason 
df_closed['reason'].value_counts()

Enforcement & Abandoned Vehicles     14908
Code Enforcement                     10437
Street Cleaning                       8477
Sanitation                            5993
Highway Maintenance                   5032
Signs & Signals                       2202
Street Lights                         1774
Recycling                             1690
Housing                               1529
Needle Program                        1298
Building                              1293
Park Maintenance & Safety             1001
Trees                                  762
Animal Issues                          580
Environmental Services                 560
Employee & General Comments            366
Health                                 344
Graffiti                               297
Administrative & General Requests      261
Notification                           141
Traffic Management & Engineering       113
Abandoned Bicycle                      108
Sidewalk Cover / Manhole                53
Catchbasin                              40
Fire Hydrant                            26
Noise Disturbance                       24
Programs                                23
Pothole                                 22
Air Pollution Control                   13
Operations                              11
Neighborhood Services Issues             9
Weights and Measures                     8
Cemetery                                 7
Generic Noise Disturbance                7
Parking Complaints                       5
Fire Department                          3
Office of The Parking Clerk              2
Billing                                  1
Name: reason, dtype: int64

There were 38 different types of requests from January 2022 - March 2022, the top three with most requests being Enforcement & Abandoned Vehicles with 14,908 requests, Code Enforcement with 10,437 requests, then Street Cleaning with 8,477 requests.

# top case request reason by neighborhood 
df_closed.groupby(['neighborhood'])['reason'].describe()

# get counts for each request reason by neighborhood 
reason_by_neighborhood = df_closed.groupby(['neighborhood', 'reason'])['duration_in_days'].describe()[['count']]
reason_by_neighborhood

# run this cell to write the reason by neighborhood to a csv to see all rows of data 
reason_by_neighborhood.to_csv('reasons_by_neighborhood.csv')

# let's take a look at the South End neighborhood specifically 
south_end_df = df_closed.loc[(df_closed['neighborhood'] == 'South End')]
south_end_df.groupby(['reason'])['duration_in_days'].describe()[['count']]

What types of cases typically take the longest to resolve?

To answer this question, let’s take a look at the duration_in_days and reason columns.

# what types of cases typically take the longest 
# case_duration by reason 

sns.catplot(x="reason", y="duration_in_days", kind="box", data=df_closed,)

<seaborn.axisgrid.FacetGrid at 0x16ba4bf10>

# The chart is kind of difficult to read... 
# Let's fix the size of the chart and flip the labels on the x-axis 

sns.catplot(y="reason", x="duration_in_days", kind="box", data=df_closed,
            height = 8, aspect = 1.25)

<seaborn.axisgrid.FacetGrid at 0x16bfc04c0>

Box plots display the five-number-summary, which includes: the minimum, the maximum, the sample median, and the first and third quartiles.

The box plot shows the distribution duration_in_days in a way that allows comparisions between case reasons. Box plots show the distribution of a numerical variable broken down by a categorical variable.

The box shows the quartiles of the duration_in_days and the whiskers extend to show the rest of the distribution (minimum and maximum). Points that are shown outside of the whiskers are determined to be outliers. The line inside the box is the median.

# descriptive statistics for duration_in_days by case reason 
# box plot in table form 
pd.set_option('display.max_columns', None)
df_closed.groupby(['reason'])['duration_in_days'].describe()

Graffiti cases take on average take the longests time to resolve, 60.796 days.

Do cases typically take longer in one neighborhood over another?

# do cases typically take longer in one neighborhood over another?

sns.catplot(y="neighborhood", x="duration_in_days", kind="box", data=df_closed,
            height = 8, aspect = 1.25)

<seaborn.axisgrid.FacetGrid at 0x16c38f880>

The box plot above shows several outliers for each category (neighborhood) making it difficult to read and quite overwhelming.

Let’s display the information in table form.

# in table form 
df_closed.groupby(['neighborhood'])['duration_in_days'].describe()

In January 2022 - March 2022, cases took the longest in Chestnut Hill. Cases typically lasted on average 19.075 days but there were only 4 cases located in Chestnut Hill during this time. Smaller sample sizes could mean more variability (look at standard deviation to explain the spread of observations).

We can further look at the population of Chestnut Hill versus the other neighborhoods to try and make sense of this low case count. Additionally, we can broaden the time period of the cases to see if Chestnut Hill still has a low case count.

From the table above we can see how long cases take by each neighborhood, it would be interesting to further breakdown by case reason for each neighborhood.

Wrap Up, Next Steps

Further analysis could be done using the 311 dataset. Using the 311 data from previous years, we can see how number of requests have changed over the years, or how case duration may have changed over the years.

Since most requests have latitude and longitude coordinates it could be interesting to plot each case request on a map to see if there are clusters of requests in certain locations.

Next steps could include gathering demographic data to overlay on top of the 311 dataset for further analysis. Another possible next step would be to build a model to predict how long a request could take given the request reason, subject, location, source, etc.

Click here to download this Jupyter Notebook (make sure you are signed in with your BU email)!