Lecture 5 – Data 100, Spring 2023¶

Lisa Yan
Content by Lisa Yan, Will Fithian, Joseph Gonzalez, Deborah Nolan, Sam Lau

In [1]:
import numpy as np
import pandas as pd
In [2]:
import matplotlib.pyplot as plt
import seaborn as sns
#%matplotlib inline
plt.rcParams['figure.figsize'] = (12, 9)

sns.set()
sns.set_context('talk')
np.set_printoptions(threshold=20, precision=2, suppress=True)
pd.set_option('display.max_rows', 30)
pd.set_option('display.max_columns', None)
pd.set_option('display.precision', 2)
# This option stops scientific notation for pandas
pd.set_option('display.float_format', '{:.2f}'.format)




Different File Types¶

There are many file types for storing structured data: CSV, TSV, JSON, XML, ASCII, SAS...

  • Documentation will be your best friend to understand how to process many of these file types.
  • In lecture, we will cover TSV and JSON since pandas supports them out-of-box.


TSV¶

TSV (Tab-Separated Values) files are very similar to CSVs, but now items are delimited by tabs.

Let's check out cdc_tuberculosis.tsv, which is the same data but now in a TSV.

  1. To the Jupyter view!

  2. To the Python view!

    Quick Python reminders:

    • Python's print() prints each string (including the newline), and an additional newline on top of that.
    • We use the repr() function to return the raw string with all special characters:
In [3]:
with open("data/cdc_tuberculosis.tsv", "r") as f:
    i = 0
    for row in f:
        print(repr(row)) # print raw strings
        i += 1
        if i > 3:
            break

'\tNo. of TB cases\t\t\tTB incidence\t\t\n'
'U.S. jurisdiction\t2019\t2020\t2021\t2019\t2020\t2021\n'
'Total\t"8,900"\t"7,173"\t"7,860"\t2.71\t2.16\t2.37\n'
'Alabama\t87\t72\t92\t1.77\t1.43\t1.83\n'

The pd.read_csv function also reads in TSVs if we specify the delimiter with parameter sep='\t' (documentation).

In [4]:
tuberculosis_df_tsv = pd.read_csv("data/cdc_tuberculosis.tsv", sep='\t')
tuberculosis_df_tsv.head()
Out[4]:
Unnamed: 0 No. of TB cases Unnamed: 2 Unnamed: 3 TB incidence Unnamed: 5 Unnamed: 6
0 U.S. jurisdiction 2019 2020 2021 2019.00 2020.00 2021.00
1 Total 8,900 7,173 7,860 2.71 2.16 2.37
2 Alabama 87 72 92 1.77 1.43 1.83
3 Alaska 58 58 58 7.91 7.92 7.92
4 Arizona 183 136 129 2.51 1.89 1.77

Side note: there was a question last time on how pandas differentiates a comma delimiter vs. a comma within the field itself, e.g., 8,900. Check out the documentation for the quotechar parameter.


JSON¶

The City of Berkeley Open Data website has a dataset with COVID-19 Confirmed Cases among Berkeley residents by date.

Let's first check out this website.

Next, let's download this file, saving it as a JSON (note the source URL file type):

Reproducible Data Science¶

In the interest of reproducible data science we will download the data programatically. We have defined some helper functions in the ds100_utils.py file. I can then reuse these helper functions in many different notebooks.

In [5]:
# just run this cell
from ds100_utils import fetch_and_cache

Occasionally, you will want to modify code that you have imported from a local Python library. To reimport those modifications you can either use the python importlib library:

from importlib import reload
reload(utils)

or use iPython magic which will intelligently import code when files change:

%load_ext autoreload
%autoreload 2
In [6]:
covid_file = fetch_and_cache(
    "https://data.cityofberkeley.info/api/views/xn6j-b766/rows.json?accessType=DOWNLOAD",
    "confirmed-cases.json",
    force=False)
covid_file          # a file path wrapper object

Using cached version that was downloaded (UTC): Tue Jan 31 14:33:04 2023
Out[6]:
PosixPath('data/confirmed-cases.json')


File size¶

Often, I like to start my analysis by getting a rough estimate of the size of the data. This will help inform the tools I use and how I view the data. If it is relatively small I might use a text editor or a spreadsheet to look at the data. If it is larger, I might jump to more programmatic exploration or even used distributed computing tools.

However here we will use Python tools to probe the file.

Since these seem to be text files I might also want to investigate the number of lines, which often corresponds to the number of records.

In [7]:
import os

print(covid_file, "is", os.path.getsize(covid_file) / 1e6, "MB")

with open(covid_file, "r") as f:
    print(covid_file, "is", sum(1 for l in f), "lines.")

data/confirmed-cases.json is 0.183341 MB
data/confirmed-cases.json is 1559 lines.


File contents¶

Because we have a text file in a visual IDE like Jupyter/DataHub, I'm going to visually explore the data via the built-in file explorer.

  1. To the Python view...?
In [8]:
with open(covid_file, "r") as f:
    i = 0
    for row in f:
        print(repr(row)) # print raw strings
        i += 1
        if i > 5:
            break

'{\n'
'  "meta" : {\n'
'    "view" : {\n'
'      "id" : "xn6j-b766",\n'
'      "name" : "COVID-19 Confirmed Cases",\n'
'      "assetType" : "dataset",\n'
  1. To the text editor view???

  1. Back to the Python view.

    In order to load the JSON file into pandas, Let's first do some EDA with the Python json package to understand the particular structure of this JSON file so that we can decide what (if anything) to load into Pandas.


EDA: Digging into JSON¶

Python has relatively good support for JSON data since it closely matches the internal python object model. In the following cell we import the entire JSON datafile into a python dictionary using the json package.

In [9]:
import json

with open(covid_file, "rb") as f:
    covid_json = json.load(f)

The covid_json variable is now a dictionary encoding the data in the file:

In [10]:
type(covid_json)
Out[10]:
dict

Examine what keys are in the top level json object¶

We can list the keys to determine what data is stored in the object.

In [11]:
covid_json.keys()
Out[11]:
dict_keys(['meta', 'data'])

Observation: The JSON dictionary contains a meta key which likely refers to meta data (data about the data). Meta data often maintained with the data and can be a good source of additional information.


We can investigate the meta data further by examining the keys associated with the metadata.

In [12]:
covid_json['meta'].keys()
Out[12]:
dict_keys(['view'])

The meta key contains another dictionary called view. This likely refers to meta-data about a particular "view" of some underlying database. We will learn more about views when we study SQL later in the class.

In [13]:
covid_json['meta']['view'].keys()
Out[13]:
dict_keys(['id', 'name', 'assetType', 'attribution', 'averageRating', 'category', 'createdAt', 'description', 'displayType', 'downloadCount', 'hideFromCatalog', 'hideFromDataJson', 'newBackend', 'numberOfComments', 'oid', 'provenance', 'publicationAppendEnabled', 'publicationDate', 'publicationGroup', 'publicationStage', 'rowsUpdatedAt', 'rowsUpdatedBy', 'tableId', 'totalTimesRated', 'viewCount', 'viewLastModified', 'viewType', 'approvals', 'clientContext', 'columns', 'grants', 'metadata', 'owner', 'query', 'rights', 'tableAuthor', 'tags', 'flags'])

Notice that this a nested/recursive data structure. As we dig deeper we reveal more and more keys and the corresponding data:

meta
|-> data
    | ... (haven't explored yet)
|-> view
    | -> id
    | -> name
    | -> attribution 
    ...
    | -> description
    ...
    | -> columns
    ...

There is a key called description in the view sub dictionary. This likely contains a description of the data:

In [14]:
print(covid_json['meta']['view']['description'])

Counts of confirmed COVID-19 cases among Berkeley residents by date. As of 6/21/22, this dataset will be updated weekly instead of daily. As of 11/14/22, this dataset only includes PCR cases.

Examining the Data Field for Records¶

We can look at a few entries in the data field. This is what we'll load into Pandas.

In [15]:
for i in range(3):
    print(f"{i:03} | {covid_json['data'][i]}")

000 | ['row-2yyp~r8a3~phgq', '00000000-0000-0000-F944-963C2BD56F87', 0, 1674521616, None, 1674521616, None, '{ }', '2019-12-01T00:00:00', '0', '0']
001 | ['row-svsf_gzh2~cz9t', '00000000-0000-0000-BA8E-0D8297E66451', 0, 1674521616, None, 1674521616, None, '{ }', '2019-12-02T00:00:00', '0', '0']
002 | ['row-w244_ivf6-rdcu', '00000000-0000-0000-A608-DB4DF9DB1B16', 0, 1674521616, None, 1674521616, None, '{ }', '2019-12-03T00:00:00', '0', '0']

Observations:

  • These look like equal-length records, so maybe data is a table!
  • But what do each of values in the record mean? Where can we find column headers?

Back to the metadata.


Columns Metadata¶

Another potentially useful key in the metadata dictionary is the columns. This returns a list:

In [16]:
type(covid_json['meta']['view']['columns'])
Out[16]:
list

Let's go back to the file explorer.

Based on the contents of this key, what are reasonable names for each column in the data table?

Summary of exploring the JSON file¶

  1. The above metadata tells us a lot about the columns in the data including column names, potential data anomalies, and a basic statistic.
  2. Because of its non-tabular structure, JSON makes it easier (than CSV) to create self-documenting data, meaning that information about the data is stored in the same file as the data.
  3. Self documenting data can be helpful since it maintains its own description and these descriptions are more likely to be updated as data changes.


3. Finally, read data into a pandas DataFrame¶

After our above EDA, let's finally go about loading the data (not the metadata) into a pandas dataframe.

In the following block of code we:

  1. Translate the JSON records into a dataframe:

    • fields: covid_json['meta']['view']['columns']
    • records: covid_json['data']

  2. Remove columns that have no metadata description. This would be a bad idea in general but here we remove these columns since the above analysis suggests that they are unlikely to contain useful information.

  3. Examine the tail of the table.
In [17]:
# Load the data from JSON and assign column titles
covid = pd.DataFrame(
    covid_json['data'],
    columns=[c['name'] for c in covid_json['meta']['view']['columns']])

covid.tail()
Out[17]:
sid id position created_at created_meta updated_at updated_meta meta Date New Cases Cumulative Cases
1144 row-vvwh-m4y7.mr6x 00000000-0000-0000-8DE9-36E51F431D12 0 1674521616 None 1674521616 None { } 2023-01-18T00:00:00 6 22182
1145 row-zex5.4t7g-ehhz 00000000-0000-0000-7A55-165F57F862F6 0 1674521616 None 1674521616 None { } 2023-01-19T00:00:00 0 22182
1146 row-gh6h_y9cb_6knb 00000000-0000-0000-E039-77F6B2766B82 0 1674521616 None 1674521616 None { } 2023-01-20T00:00:00 0 22182
1147 row-cb4r~rtkn.xg9x 00000000-0000-0000-CCC4-D0C4B5D12422 0 1674521616 None 1674521616 None { } 2023-01-21T00:00:00 0 22182
1148 row-4en9-p4vi.fq5m 00000000-0000-0000-7F40-89071F462EE8 0 1674521616 None 1674521616 None { } 2023-01-22T00:00:00 0 22182



Data Faithfulness: Mauna Loa CO2 data¶

CO2 concentrations have been monitored at Mauna Loa Observatory since 1958 (website link).

In [18]:
co2_file = "data/co2_mm_mlo.txt"

Let's do some EDA!!

How do we read the file into Pandas?¶

Let's instead check out this file with JupyterLab.

  • Note it's a .txt file.
  • Do we trust this file extension?
  • What structure is it?






Looking at the first few lines of the data, we spot some relevant characteristics:

  • The values are separated by white space, possibly tabs.
  • The data line up down the rows. For example, the month appears in 7th to 8th position of each line.
  • The 71st and 72nd lines in the file contain column headings split over two lines.

We can use read_csv to read the data into a Pandas data frame, and we provide several arguments to specify that the separators are white space, there is no header (we will set our own column names), and to skip the first 72 rows of the file.

In [19]:
co2 = pd.read_csv(
    co2_file, header = None, skiprows = 72,
    sep = r'\s+'       #delimiter for continuous whitespace (stay tuned for regex next lecture))
)
co2.head()
Out[19]:
0 1 2 3 4 5 6
0 1958 3 1958.21 315.71 315.71 314.62 -1
1 1958 4 1958.29 317.45 317.45 315.29 -1
2 1958 5 1958.38 317.50 317.50 314.71 -1
3 1958 6 1958.46 -99.99 317.10 314.85 -1
4 1958 7 1958.54 315.86 315.86 314.98 -1

Congratulations! You've wrangled your first set of real world data!


...But our columns aren't named. We need to do more EDA.



¶

Real World Example: Wrangling CO2 Measurements¶

Let's continue looking at the CO2 concentrations at Mauna Loa Observatory:

In [20]:
co2.head()
Out[20]:
0 1 2 3 4 5 6
0 1958 3 1958.21 315.71 315.71 314.62 -1
1 1958 4 1958.29 317.45 317.45 315.29 -1
2 1958 5 1958.38 317.50 317.50 314.71 -1
3 1958 6 1958.46 -99.99 317.10 314.85 -1
4 1958 7 1958.54 315.86 315.86 314.98 -1

We need column names!

Exploring Variable Feature Types¶

The NOAA webpage might have some useful tidbits (in this case it doesn't). Let's go back to the raw data file to identify each feature.



We'll rerun pd.read_csv, but this time with some custom column names.

In [21]:
co2 = pd.read_csv(
    co2_file, header = None, skiprows = 72,
    sep = '\s+', #regex for continuous whitespace (next lecture)
    names = ['Yr', 'Mo', 'DecDate', 'Avg', 'Int', 'Trend', 'Days']
)
co2.head()
Out[21]:
Yr Mo DecDate Avg Int Trend Days
0 1958 3 1958.21 315.71 315.71 314.62 -1
1 1958 4 1958.29 317.45 317.45 315.29 -1
2 1958 5 1958.38 317.50 317.50 314.71 -1
3 1958 6 1958.46 -99.99 317.10 314.85 -1
4 1958 7 1958.54 315.86 315.86 314.98 -1



Let's start exploring!!¶

Scientific studies tend to have very clean data, right...? Let's jump right in and make a time series plot of CO2 monthly averages.

In [22]:
sns.lineplot(x='DecDate', y='Avg', data=co2);

The code above uses the seaborn plotting library (abbreviated sns). We won't cover this library in detail until next week, so focus on the plots themselves, not the code used to create them.

Yikes! Plotting the data uncovered a problem. It looks like we have some missing values. What happened here?

In [23]:
co2.head()
Out[23]:
Yr Mo DecDate Avg Int Trend Days
0 1958 3 1958.21 315.71 315.71 314.62 -1
1 1958 4 1958.29 317.45 317.45 315.29 -1
2 1958 5 1958.38 317.50 317.50 314.71 -1
3 1958 6 1958.46 -99.99 317.10 314.85 -1
4 1958 7 1958.54 315.86 315.86 314.98 -1
In [24]:
co2.tail()
Out[24]:
Yr Mo DecDate Avg Int Trend Days
733 2019 4 2019.29 413.32 413.32 410.49 26
734 2019 5 2019.38 414.66 414.66 411.20 28
735 2019 6 2019.46 413.92 413.92 411.58 27
736 2019 7 2019.54 411.77 411.77 411.43 23
737 2019 8 2019.62 409.95 409.95 411.84 29

Some data have unusual values like -1 and -99.99.

Let's check the description at the top of the file again.

  1. -1 signifies a missing value for the number of days Days the equipment was in operation that month.
  2. -99.99 denotes a missing monthly average Avg

How can we fix this? First, let's explore other aspects of our data. Understanding our data will help us decide what to do with the missing values.


Quality Checks: Reasoning about the data¶

First, we consider the shape of the data. How many rows should we have?

  • If chronological order, we should have one record per month.
  • Data from March 1958 to August 2019.
  • We should have $ 12 \times (2019-1957) - 2 - 4 = 738 $ records.
In [25]:
co2.shape
Out[25]:
(738, 7)

Nice!! The number of rows (i.e. records) match our expectations.



Let's now check the quality of each feature.

Understanding Missing Value 1: Days¶

Days is a time field, so let's analyze other time fields to see if there is an explanation for missing values of days of operation.

Let's start with months Mo.

Are we missing any records? The number of months should have 62 or 61 instances (March 1957-August 2019).

In [26]:
co2["Mo"].value_counts().sort_index()
Out[26]:
1     61
2     61
3     62
4     62
5     62
6     62
7     62
8     62
9     61
10    61
11    61
12    61
Name: Mo, dtype: int64

As expected Jan, Feb, Sep, Oct, Nov, and Dec have 61 occurrences and the rest 62.



Next let's explore days Days itself, which is the number of days that the measurement equipment worked.

In [27]:
sns.displot(co2['Days']);
plt.title("Distribution of days feature")
plt.show() # suppresses unneeded plotting output

In terms of data quality, a handful of months have averages based on measurements taken on fewer than half the days. In addition, there are nearly 200 missing values--that's about 27% of the data!



Finally, let's check the last time feature, year Yr.

Let's check to see if there is any connection between missingness and the year of the recording.

In [28]:
sns.scatterplot(x="Yr", y="Days", data=co2);
plt.title("Day field by Year"); # the ; suppresses output

Observations:

  • All of the missing data are in the early years of operation.
  • It appears there may have been problems with equipment in the mid to late 80s.

Potential Next Steps:

  • Confirm these explanations through documentation about the historical readings.
  • Maybe drop earliest recordings? However, we would want to delay such action until after we have examined the time trends and assess whether there are any potential problems.



Understanding Missing Value 2: Avg¶

Next, let's return to the -99.99 values in Avg to analyze the overall quality of the CO2 measurements.

In [29]:
# Histograms of average CO2 measurements
sns.displot(co2['Avg']);

The non-missing values are in the 300-400 range (a regular range of CO2 levels).

We also see that there are only a few missing Avg values (<1% of values). Let's examine all of them:

In [30]:
co2[co2["Avg"] < 0]
Out[30]:
Yr Mo DecDate Avg Int Trend Days
3 1958 6 1958.46 -99.99 317.10 314.85 -1
7 1958 10 1958.79 -99.99 312.66 315.61 -1
71 1964 2 1964.12 -99.99 320.07 319.61 -1
72 1964 3 1964.21 -99.99 320.73 319.55 -1
73 1964 4 1964.29 -99.99 321.77 319.48 -1
213 1975 12 1975.96 -99.99 330.59 331.60 0
313 1984 4 1984.29 -99.99 346.84 344.27 2

There doesn't seem to be a pattern to these values, other than that most records also were missing Days data.

Drop, NaN, or Impute Missing Avg Data?¶

How should we address the invalid Avg data?

A. Drop records

B. Set to NaN

C. Impute using some strategy

Remember we want to fix the following plot:

In [31]:
sns.lineplot(x='DecDate', y='Avg', data=co2)
plt.title("CO2 Average By Month");

Since we are plotting Avg vs DecDate, we should just focus on dealing with missing values for Avg.

Let's consider a few options:

  1. Drop those records
  2. Replace -99.99 with NaN
  3. Substitute it with a likely value for the average CO2?

What do you think are the pros and cons of each possible action?



Let's examine each of these three options.

In [32]:
# 1. Drop missing values
co2_drop = co2[co2['Avg'] > 0]

# 2. Replace NaN with -99.99
co2_NA = co2.replace(-99.99, np.NaN)

We'll also use a third version of the data. First, we note that the dataset already comes with a substitute value for the -99.99.

From the file description:

The interpolated column includes average values from the preceding column (average) and interpolated values where data are missing. Interpolated values are computed in two steps...

The Int feature has values that exactly match those in Avg, except when Avg is -99.99, and then a reasonable estimate is used instead. So, the third version of our data will use the Int feature instead of Avg.

In [33]:
# 3. Use interpolated column which estimates missing Avg values
co2_impute = co2.copy()
co2_impute['Avg'] = co2['Int']


What's a reasonable estimate?

To answer this question, let's zoom in on a short time period, say the measurements in 1958 (where we know we have two missing values).

In [34]:
# results of plotting data in 1958

def line_and_points(data, ax, title):
    # assumes single year, hence Mo
    ax.plot('Mo', 'Avg', data=data)
    ax.scatter('Mo', 'Avg', data=data)
    ax.set_xlim(2, 13)
    ax.set_title(title)
    ax.set_xticks(np.arange(3, 13))

def data_year(data, year):
    return data[data["Yr"] == 1958]
    
# uses matplotlib subplots
# you may see more next week; focus on output for now
fig, axes = plt.subplots(ncols = 3, figsize=(12, 4), sharey=True)

year = 1958
line_and_points(data_year(co2_drop, year), axes[0], title="1. Drop Missing")
line_and_points(data_year(co2_NA, year), axes[1], title="2. Missing Set to NaN")
line_and_points(data_year(co2_impute, year), axes[2], title="3. Missing Interpolated")

fig.suptitle(f"Monthly Averages for {year}")
plt.tight_layout()

In the big picture since there are only 7 Avg values missing (<1% of 738 months), any of these approaches would work.

However there is some appeal to option C: Imputing:

  • Shows seasonal trends for CO2
  • We are plotting all months in our data as a line plot


Let's replot our original figure with option 3:

In [35]:
sns.lineplot(x='DecDate', y='Avg', data=co2_impute)
plt.title("CO2 Average By Month, Imputed");

Looks pretty close to what we see on the NOAA website!

Presenting the data: A Discussion on Data Granularity¶

From the description:

  • monthly measurements are averages of average day measurements.
  • The NOAA GML website has datasets for daily/hourly measurements too.

The data you present depends on your research question.

How do CO2 levels vary by season?

  • You might want to keep average monthly data.

Are CO2 levels rising over the past 50+ years, consistent with global warming predictions?

  • You might be happier with a coarser granularity of average year data!
In [36]:
co2_year = co2_impute.groupby('Yr').mean()
sns.lineplot(x='Yr', y='Avg', data=co2_year)
plt.title("CO2 Average By Year");

Indeed, we see a rise by nearly 100 ppm of CO2 since Mauna Loa began recording in 1958.