<jupyter_start><jupyter_text>COVID-19 in Virginia Time series analysis of daily new cases The data for this analysis was retrieved from the Virginia Department of Health's (VDH) COVID-19 data page. Specifically, I am using the [Cases Summary Public Use Dataset](https://data.virginia.gov/Government/VDH-COVID-19-PublicUseDataset-Cases-Summary/3ves-mpfr), which is available for download with daily updated data.For this analysis, I will utilize the `pandas`, `numpy`, and `matplotlib` libraries.<jupyter_code># load libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# use pandas to load the data set downloaded from the VDH
covid = pd.read_csv(
    "VDH-COVID-19-PublicUseDataset-Cases-Summary.csv"
)

# examine the first 5 rows of the data set
covid.head(5)
# examine the last 5 rows of the data set
covid.tail(5)<jupyter_output><empty_output><jupyter_text>We can see that our data set includes three fields:- The report date, beginning on November 28, 2020 and ending on July 23, 2021- A summary indicator...we are interested in the "New Cases" indicator- A summary indicator value that gives the value for the specified indicator on the specified dateTo conduct our analysis, we need to filter our data set to include only the rows for which the Summary Indicator field contains "New Cases."<jupyter_code># create a new data frame for the new case data
new_cases = covid[covid["Summary Indicator"] == "New Cases"].copy()

# examine the first 5 rows of the new cases data set
new_cases.head(5)<jupyter_output><empty_output><jupyter_text>Before we being any time series analysis, we need to ensure that our date column is the proper data type.<jupyter_code># examine the data types
new_cases.dtypes<jupyter_output><empty_output><jupyter_text>Since the Report Date field is an object, in this case a string, we need to convert it to a proper `datetime` type. Additionally, the Summary Indicator Value field is a string; we need to remove the ',' character and convert it to a numeric data type.We also want to set the index of the data frame to the 'Report Date' field.<jupyter_code># convert 'Report Date' to datetime type
new_cases["Report Date"] = pd.to_datetime(new_cases["Report Date"])

# remove the ',' character and convert 'Summary Indicator Value' to numeric type
new_cases = new_cases.replace(',', '', regex = True)
new_cases["Summary Indicator Value"] = pd.to_numeric(new_cases["Summary Indicator Value"])

# ensure the dates are in the proper order
new_cases = new_cases.sort_values(by = ["Report Date"])

# set the index to the 'Report Date' field
new_cases = new_cases.set_index("Report Date")

# examine the data types
new_cases.dtypes<jupyter_output><empty_output><jupyter_text>Now let's visualize the daily new COVID-19 case data for Virginia using `matplotlib`. We start by defining a chart formatting function so we are not repeating a lot of formatting syntax throughout the analysis.<jupyter_code># chart formatting funcion
def ChartFormatter(title, legend = True):
    # set the x- and y-axis tick parameters
    ax.xaxis.set_tick_params(labelsize = 14)
    ax.yaxis.set_tick_params(labelsize = 14)

    # Use automatic StrMethodFormatter
    #ax.yaxis.set_major_formatter("{x:,.0f}")

    # hide the right and top spines
    ax.spines["right"].set_visible(False)
    ax.spines["top"].set_visible(False)

    if (legend):
        # add a legend
        ax.legend()
    
    ##### Figure object properties #####
    # give the chart a title
    fig.suptitle(title, fontsize = 20)

    # set the size of the chart
    fig.set_size_inches(12, 8)

# create variables for the Figure and Axes objects
fig, ax = plt.subplots()

##### Axes object properties #####
# create two line plots
ax.plot(new_cases["Summary Indicator Value"])

# set the x- and y-axis labels
ax.xaxis.set_label_text("Date", fontsize = 16)
ax.yaxis.set_label_text("New Cases", fontsize = 16)

ChartFormatter("Daily New Cases of COVID-19 in Virginia", legend = False)<jupyter_output><empty_output><jupyter_text>Analyzing the TrendThe interesting piece here is determining the optimal window for the moving average. Most media outlets report a 7-day moving average, but is that the most accurate? We will start there, and look at other values for our window size to determine whether we can find something better.<jupyter_code>new_cases["MovingAverage"] = new_cases.iloc[:, 1].rolling(7).mean()
new_cases["SquaredError"] = (new_cases["Summary Indicator Value"] - new_cases["MovingAverage"]) ** 2<jupyter_output><empty_output><jupyter_text>Now we can visualize the daily new case data along with our moving average.<jupyter_code># remove "na" rows
new_cases7 = new_cases.dropna(how = "any", inplace = True)

# create variables for the Figure and Axes objects
fig, ax = plt.subplots()

##### Axes object properties #####
# create two line plots
ax.plot(new_cases["Summary Indicator Value"], label = "Daily New Cases")
ax.plot(new_cases["MovingAverage"], label = "Moving Average")

# set the x- and y-axis labels
ax.xaxis.set_label_text("Date", fontsize = 16)
ax.yaxis.set_label_text("New Cases", fontsize = 16)

ChartFormatter("Daily New COVID-19 Cases in Virginia")
print(new_cases["SquaredError"].sum())<jupyter_output>98782913.3877551<jupyter_text>The chart above shows that the 7-day moving average does a very nice job of capturing the trend in the daily new cases data, so I suspect it will be tough to beat. Let's take a look at a 10-day moving average.<jupyter_code>new_cases["MovingAverage"] = new_cases.iloc[:, 1].rolling(10).mean()
new_cases["SquaredError"] = (new_cases["Summary Indicator Value"] - new_cases["MovingAverage"]) ** 2

# remove "na" rows
new_cases10 = new_cases.dropna(how = "any", inplace = True)

# create variables for the Figure and Axes objects
fig, ax = plt.subplots()

##### Axes object properties #####
# create two line plots
ax.plot(new_cases["Summary Indicator Value"], label = "Daily New Cases")
ax.plot(new_cases["MovingAverage"], label = "Moving Average")

# set the x- and y-axis labels
ax.xaxis.set_label_text("Date", fontsize = 16)
ax.yaxis.set_label_text("New Cases", fontsize = 16)

ChartFormatter("Daily New COVID-19 Cases in Virginia")
print(new_cases["SquaredError"].sum())<jupyter_output>103703149.95<jupyter_text>The 10-day moving average doesn't do as good of a job at capturing the trend and it's sum of squared errors is higher. Let's try 14 days for the window.<jupyter_code>new_cases["MovingAverage"] = new_cases.iloc[:, 1].rolling(14).mean()
new_cases["SquaredError"] = (new_cases["Summary Indicator Value"] - new_cases["MovingAverage"]) ** 2

# remove "na" rows
new_cases14 = new_cases.dropna(how = "any", inplace = True)

# create variables for the Figure and Axes objects
fig, ax = plt.subplots()

##### Axes object properties #####
# create two line plots
ax.plot(new_cases["Summary Indicator Value"], label = "Daily New Cases")
ax.plot(new_cases["MovingAverage"], label = "Moving Average")

# set the x- and y-axis labels
ax.xaxis.set_label_text("Date", fontsize = 16)
ax.yaxis.set_label_text("New Cases", fontsize = 16)

ChartFormatter("Daily New COVID-19 Cases in Virginia")
print(new_cases["SquaredError"].sum())<jupyter_output>113369660.24489795<jupyter_text>From the 10- and 14-day moving averages, it appears that going past 7 days produces greater error. Let's go in the opposite direction and try a 5-day moving average.<jupyter_code>new_cases["MovingAverage"] = new_cases.iloc[:, 1].rolling(5).mean()
new_cases["SquaredError"] = (new_cases["Summary Indicator Value"] - new_cases["MovingAverage"]) ** 2

# remove "na" rows
new_cases5 = new_cases.dropna(how = "any", inplace = True)

# create variables for the Figure and Axes objects
fig, ax = plt.subplots()

##### Axes object properties #####
# create two line plots
ax.plot(new_cases["Summary Indicator Value"], label = "Daily New Cases")
ax.plot(new_cases["MovingAverage"], label = "Moving Average")

# set the x- and y-axis labels
ax.xaxis.set_label_text("Date", fontsize = 16)
ax.yaxis.set_label_text("New Cases", fontsize = 16)

ChartFormatter("Daily New COVID-19 Cases in Virginia")
print(new_cases["SquaredError"].sum())<jupyter_output>78442693.79999998<jupyter_text>The 5-day moving average is our best try yet; it captures the trend nicely and has a lower sum of squared errors than the 7-day moving average. Now that we have a good window size, let's use it and get our first residuals.<jupyter_code># calculate the first residual and the first residual percentage
new_cases["R1"] = new_cases["Summary Indicator Value"] - new_cases["MovingAverage"]
new_cases["R1Error"] = abs((new_cases["Summary Indicator Value"] - new_cases["R1"]) / new_cases["Summary Indicator Value"])

new_cases
# create variables for the Figure and Axes objects
fig, ax = plt.subplots()

##### Axes object properties #####
# create two line plots
ax.plot(new_cases["R1"])

# set the x- and y-axis labels
ax.xaxis.set_label_text("Date", fontsize = 16)
ax.yaxis.set_label_text("New Cases", fontsize = 16)

ChartFormatter("Residual Daily New Cases Accounting for Trend", legend = False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10829417549994536, 'fraction_numerical': 0.01879576002622664, 'mean_word_length': 5.408963585434174, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 41, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '447629', 'license': 'CC0-1.0', 'n_tokens_mistral': 2984, 'n_tokens_neox': 2756, 'n_words': 1007, 'path': 'covid/covid.ipynb', 'repo_name': 'jsellis-wm/jsellis-wm.github.io', 'size': '9151'}