%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

"""
Explanation: Idomatic Pandas

Q: How do I make my pandas code faster with parallelism?
A: You don’t need parallelism, you can use Pandas better.
-- Matthew Rocklin

Now that we have been exposed to the basic functionality of pandas, lets explore some more advanced features that will be useful when addressing more complex data management tasks.
As most statisticians/data analysts will admit, often the lion's share of the time spent implementing an analysis is devoted to preparing the data itself, rather than to coding or running a particular model that uses the data. This is where Pandas and  Python's standard library are beneficial, providing high-level, flexible, and efficient tools for manipulating your data as needed.
As you may already have noticed, there are sometimes mutliple ways to achieve the same goal using pandas. Importantly, some approaches are better than others, in terms of performance, readability and ease of use. We will cover some important ways of maximizing your pandas efficiency.
End of explanation
"""


cdystonia = pd.read_csv("../data/cdystonia.csv", index_col=None)
cdystonia.head()

"""
Explanation: Reshaping DataFrame objects
In the context of a single DataFrame, we are often interested in re-arranging the layout of our data. 
This dataset in from Table 6.9 of Statistical Methods for the Analysis of Repeated Measurements by Charles S. Davis, pp. 161-163 (Springer, 2002). These data are from a multicenter, randomized controlled trial of botulinum toxin type B (BotB) in patients with cervical dystonia (spasmodic torticollis) from nine U.S. sites.

Randomized to placebo (N=36), 5000 units of BotB (N=36), 10,000 units of BotB (N=37)
Response variable: total score on Toronto Western Spasmodic Torticollis Rating Scale (TWSTRS), measuring severity, pain, and disability of cervical dystonia (high scores mean more impairment)
TWSTRS measured at baseline (week 0) and weeks 2, 4, 8, 12, 16 after treatment began
End of explanation
"""


stacked = cdystonia.stack()
stacked

"""
Explanation: This dataset includes repeated measurements of the same individuals (longitudinal data). Its possible to present such information in (at least) two ways: showing each repeated measurement in their own row, or in multiple columns representing multiple measurements.
The stack method rotates the data frame so that columns are represented in rows:
End of explanation
"""


stacked.unstack().head()

"""
Explanation: Have a peek at the structure of the index of the stacked data (and the data itself).
To complement this, unstack pivots from rows back to columns.
End of explanation
"""


# Write your answer here

"""
Explanation: Exercise
Which columns uniquely define a row? Create a DataFrame called cdystonia2 with a hierarchical index based on these columns.
End of explanation
"""


twstrs_wide = cdystonia2['twstrs'].unstack('obs')
twstrs_wide.head()

"""
Explanation: If we want to transform this data so that repeated measurements are in columns, we can unstack the twstrs measurements according to obs.
End of explanation
"""


cdystonia_wide = (cdystonia[['patient','site','id','treat','age','sex']]
                  .drop_duplicates()
                  .merge(twstrs_wide, right_index=True, left_on='patient', how='inner'))
cdystonia_wide.head()

"""
Explanation: We can now merge these reshaped outcomes data with the other variables to create a wide format DataFrame that consists of one row for each patient.
End of explanation
"""


(cdystonia.set_index(['patient','site','id','treat','age','sex','week'])['twstrs']
     .unstack('week').head())

"""
Explanation: A slightly cleaner way of doing this is to set the patient-level information as an index before unstacking:
End of explanation
"""


pd.melt(cdystonia_wide, id_vars=['patient','site','id','treat','age','sex'], 
        var_name='obs', value_name='twsters').head()

"""
Explanation: To convert our "wide" format back to long, we can use the melt function, appropriately parameterized. This function is useful for DataFrames where one
or more columns are identifier variables (id_vars), with the remaining columns being measured variables (value_vars). The measured variables are "unpivoted" to
the row axis, leaving just two non-identifier columns, a variable and its corresponding value, which can both be renamed using optional arguments.
End of explanation
"""


(cdystonia[['patient','site','id','treat','age','sex']]
                  .drop_duplicates()
                  .merge(twstrs_wide, right_index=True, left_on='patient', how='inner')
                  .head())

"""
Explanation: This illustrates the two formats for longitudinal data: long and wide formats. Its typically better to store data in long format because additional data can be included as additional rows in the database, while wide format requires that the entire database schema be altered by adding columns to every row as data are collected.
The preferable format for analysis depends entirely on what is planned for the data, so it is imporant to be able to move easily between them.
Method chaining
In the DataFrame reshaping section above, you probably noticed how several methods were strung together to produce a wide format table:
End of explanation
"""


cdystonia_subset = cdystonia[['patient','site','id','treat','age','sex']]
cdystonia_complete = cdystonia_subset.drop_duplicates()
cdystonia_merged = cdystonia_complete.merge(twstrs_wide, right_index=True, left_on='patient', how='inner')
cdystonia_merged.head()

"""
Explanation: This approach of seqentially calling methods is called method chaining, and despite the fact that it creates very long lines of code that must be properly justified, it allows for the writing of rather concise and readable code. Method chaining is possible because of the pandas convention of returning copies of the results of operations, rather than in-place operations. This allows methods from the returned object to be immediately called, as needed, rather than assigning the output to a variable that might not otherwise be used. For example, without method chaining we would have done the following:
End of explanation
"""


measles = pd.read_csv("../data/measles.csv", index_col=0, encoding='latin-1', parse_dates=['ONSET'])
measles.head()

"""
Explanation: This necessitates the creation of a slew of intermediate variables that we really don't need.
Let's transform another dataset using method chaining. The measles.csv file contains de-identified cases of measles from an outbreak in Sao Paulo, Brazil in 1997. The file contains rows of individual records:
End of explanation
"""


(measles.assign(AGE_GROUP=pd.cut(measles.YEAR_AGE, [0,5,10,15,20,25,30,35,40,100], right=False))
                        .groupby(['ONSET', 'AGE_GROUP']))

"""
Explanation: The goal is to summarize this data by age groups and bi-weekly period, so that we can see how the outbreak affected different ages over the course of the outbreak.
The best approach is to build up the chain incrementally. We can begin by generating the age groups (using cut) and grouping by age group and the date (ONSET):
End of explanation
"""


(measles.assign(AGE_GROUP=pd.cut(measles.YEAR_AGE, [0,5,10,15,20,25,30,35,40,100], right=False))
                        .groupby(['ONSET', 'AGE_GROUP'])
                        .size()).head(10)

"""
Explanation: What we then want is the number of occurences in each combination, which we can obtain by checking the size of each grouping:
End of explanation
"""


(measles.assign(AGE_GROUP=pd.cut(measles.YEAR_AGE, [0,5,10,15,20,25,30,35,40,100], right=False))
                        .groupby(['ONSET', 'AGE_GROUP'])
                        .size()
                        .unstack()).head(5)

"""
Explanation: This results in a hierarchically-indexed Series, which we can pivot into a DataFrame by simply unstacking:
End of explanation
"""


(measles.assign(AGE_GROUP=pd.cut(measles.YEAR_AGE, [0,5,10,15,20,25,30,35,40,100], right=False))
                        .groupby(['ONSET', 'AGE_GROUP'])
                        .size()
                        .unstack()
                        .fillna(0)).head(5)

"""
Explanation: Now, fill replace the missing values with zeros:
End of explanation
"""


case_counts_2w = (measles.assign(AGE_GROUP=pd.cut(measles.YEAR_AGE, [0,5,10,15,20,25,30,35,40,100], right=False))
                        .groupby(['ONSET', 'AGE_GROUP'])
                        .size()
                        .unstack()
                        .fillna(0)
                        .resample('2W')
                        .sum())

case_counts_2w

"""
Explanation: Finally, we want the counts in 2-week intervals, rather than as irregularly-reported days, which yields our the table of interest:
End of explanation
"""


case_counts_2w.plot(cmap='hot')

"""
Explanation: From this, it is easy to create meaningful plots and conduct analyses:
End of explanation
"""


cdystonia.pivot(index='patient', columns='obs', values='twstrs').head()

"""
Explanation: Pivoting
The pivot method allows a DataFrame to be transformed easily between long and wide formats in the same way as a pivot table is created in a spreadsheet. It takes three arguments: index, columns and values, corresponding to the DataFrame index (the row headers), columns and cell values, respectively.
For example, we may want the twstrs variable (the response variable) in wide format according to patient, as we saw with the unstacking method above:
End of explanation
"""


# Write your answer here

"""
Explanation: Exercise
Try pivoting the cdystonia DataFrame without specifying a variable for the cell values:
End of explanation
"""


cdystonia.pivot_table(index=['site', 'treat'], columns='week', values='twstrs', 
                      aggfunc=max).head(20)

"""
Explanation: A related method, pivot_table, creates a spreadsheet-like table with a hierarchical index, and allows the values of the table to be populated using an arbitrary aggregation function.
End of explanation
"""


pd.crosstab(cdystonia.sex, cdystonia.site)

"""
Explanation: For a simple cross-tabulation of group frequencies, the crosstab function (not a method) aggregates counts of data according to factors in rows and columns. The factors may be hierarchical if desired.
End of explanation
"""


vessels = pd.read_csv('../data/AIS/vessel_information.csv')
vessels.tail(10)

vessels.duplicated(subset='names').tail(10)

"""
Explanation: Data transformation
There are a slew of additional operations for DataFrames that we would collectively refer to as transformations which include tasks such as:

removing duplicate values
replacing values
grouping values.

Dealing with duplicates
We can easily identify and remove duplicate values from DataFrame objects. For example, say we want to remove ships from our vessels dataset that have the same name:
End of explanation
"""


vessels.drop_duplicates(['names']).tail(10)

"""
Explanation: These rows can be removed using drop_duplicates
End of explanation
"""


cdystonia.treat.value_counts()

"""
Explanation: Value replacement
Frequently, we get data columns that are encoded as strings that we wish to represent numerically for the purposes of including it in a quantitative analysis. For example, consider the treatment variable in the cervical dystonia dataset:
End of explanation
"""


treatment_map = {'Placebo': 0, '5000U': 1, '10000U': 2}

cdystonia['treatment'] = cdystonia.treat.map(treatment_map)
cdystonia.treatment

"""
Explanation: A logical way to specify these numerically is to change them to integer values, perhaps using "Placebo" as a baseline value. If we create a dict with the original values as keys and the replacements as values, we can pass it to the map method to implement the changes.
End of explanation
"""


scores = pd.Series([99, 76, 85, -999, 84, 95])

"""
Explanation: Alternately, if we simply want to replace particular values in a Series or DataFrame, we can use the replace method. 
An example where replacement is useful is replacing sentinel values with an appropriate numeric value prior to analysis. A large negative number is sometimes used in otherwise positive-valued data to denote missing values.
End of explanation
"""


scores.replace(-999, np.nan)

"""
Explanation: In such situations, we can use replace to substitute nan where the sentinel values occur.
End of explanation
"""


cdystonia2.treat.replace({'Placebo': 0, '5000U': 1, '10000U': 2})

"""
Explanation: We can also perform the same replacement that we used map for with replace:
End of explanation
"""


# Write your answer here

"""
Explanation: Inidcator variables
For some statistical analyses (e.g. regression models or analyses of variance), categorical or group variables need to be converted into columns of indicators--zeros and ones--to create a so-called design matrix. The Pandas function get_dummies (indicator variables are also known as dummy variables) makes this transformation straightforward.
Let's consider the DataFrame containing the ships corresponding to the transit segments on the eastern seaboard. The type variable denotes the class of vessel; we can create a matrix of indicators for this. For simplicity, lets filter out the 5 most common types of ships.
Exercise
Create a subset of the vessels DataFrame called vessels5 that only contains the 5 most common types of vessels, based on their prevalence in the dataset.
End of explanation
"""


pd.get_dummies(vessels5.type).head(10)

"""
Explanation: We can now apply get_dummies to the vessel type to create 5 indicator variables.
End of explanation
"""


cdystonia.age.describe()

"""
Explanation: Discretization
Pandas' cut function can be used to group continuous or countable data in to bins. Discretization is generally a very bad idea for statistical analysis, so use this function responsibly!
Lets say we want to bin the ages of the cervical dystonia patients into a smaller number of groups:
End of explanation
"""


pd.cut(cdystonia.age, [20,30,40,50,60,70,80,90])[:30]

"""
Explanation: Let's transform these data into decades, beginnnig with individuals in their 20's and ending with those in their 80's:
End of explanation
"""


pd.cut(cdystonia.age, [20,30,40,50,60,70,80,90], right=False)[:30]

"""
Explanation: The parentheses indicate an open interval, meaning that the interval includes values up to but not including the endpoint, whereas the square bracket is a closed interval, where the endpoint is included in the interval. We can switch the closure to the left side by setting the right flag to False:
End of explanation
"""


pd.cut(cdystonia.age, [20,40,60,80,90], labels=['young','middle-aged','old','really old'])[:30]

"""
Explanation: Since the data are now ordinal, rather than numeric, we can give them labels:
End of explanation
"""


pd.qcut(cdystonia.age, 4)[:30]

"""
Explanation: A related function qcut uses empirical quantiles to divide the data. If, for example, we want the quartiles -- (0-25%], (25-50%], (50-70%], (75-100%] -- we can just specify 4 intervals, which will be equally-spaced by default:
End of explanation
"""


quantiles = pd.qcut(vessels.max_loa, [0, 0.01, 0.05, 0.95, 0.99, 1])
quantiles[:30]

"""
Explanation: Alternatively, one can specify custom quantiles to act as cut points:
End of explanation
"""


# Write your answer here

"""
Explanation: Exercise
Use the discretized segment lengths as the input for get_dummies to create 5 indicator variables for segment length:
End of explanation
"""


cdystonia_cat = cdystonia.assign(treatment=cdystonia.treat.astype('category')).drop('treat', axis=1)
cdystonia_cat.dtypes

cdystonia_cat.treatment.head()

cdystonia_cat.treatment.cat.codes

"""
Explanation: Categorical Variables
One of the keys to maximizing performance in pandas is to use the appropriate types for your data wherever possible. In the case of categorical data--either the ordered categories as we have just created, or unordered categories like race, gender or country--the use of the categorical to encode string variables as numeric quantities can dramatically improve performance and simplify subsequent analyses.
When text data are imported into a DataFrame, they are endowed with an object dtype. This will result in relatively slow computation because this dtype runs at Python speeds, rather than as Cython code that gives much of pandas its speed. We can ameliorate this by employing the categorical dtype on such data.
End of explanation
"""


cdystonia.treat.astype('category', ordered=True).head()

"""
Explanation: This creates an unordered categorical variable. To create an ordinal variable, we can specify order=True as an argument to astype:
End of explanation
"""


cdystonia.treat.astype('category').cat.set_categories(['Placebo', '5000U', '10000U'], ordered=True).head()

"""
Explanation: However, this is not the correct order; by default, the categories will be sorted alphabetically, which here gives exactly the reverse order that we need. 
To specify an arbitrary order, we can used the set_categories method, as follows:
End of explanation
"""


cdystonia_cat.treatment.cat

"""
Explanation: Notice that we obtained set_categories from the cat attribute of the categorical variable. This is known as the category accessor, and is a device for gaining access to Categorical variables' categories, analogous to the string accessor that we have seen previously from text variables.
End of explanation
"""


cdystonia_cat['treatment'] = (cdystonia.treat.astype('category').cat
                              .set_categories(['Placebo', '5000U', '10000U', '20000U'], ordered=True))

"""
Explanation: Additional categoried can be added, even if they do not currently exist in the DataFrame, but are part of the set of possible categories:
End of explanation
"""


cdystonia_cat.treatment.cat.remove_categories('20000U').head()

"""
Explanation: To complement this, we can remove categories that we do not wish to retain:
End of explanation
"""


cdystonia_cat.treatment.cat.remove_unused_categories().head()

"""
Explanation: Or, even more simply:
End of explanation
"""


vessels_merged = (pd.read_csv('../data/AIS/vessel_information.csv', index_col=0)
                          .merge(pd.read_csv('../data/AIS/transit_segments.csv'), left_index=True, right_on='mmsi'))

vessels_merged['registered'] = vessels_merged.flag.astype('category')

%timeit vessels_merged.groupby('flag').avg_sog.mean().sort_values()

%timeit vessels_merged.groupby('registered').avg_sog.mean().sort_values()

vessels_merged[['flag','registered']].memory_usage()

"""
Explanation: For larger datasets, there is an appreciable gain in performance, both in terms of speed and memory usage.
End of explanation
"""


cdystonia_grouped = cdystonia.groupby(cdystonia.patient)

"""
Explanation: Data aggregation and GroupBy operations
One of the most powerful features of Pandas is its GroupBy functionality. On occasion we may want to perform operations on groups of observations within a dataset. For exmaple:

aggregation, such as computing the sum of mean of each group, which involves applying a function to each group and returning the aggregated results
slicing the DataFrame into groups and then doing something with the resulting slices (e.g. plotting)
group-wise transformation, such as standardization/normalization
End of explanation
"""


cdystonia_grouped

"""
Explanation: This grouped dataset is hard to visualize
End of explanation
"""


for patient, group in cdystonia_grouped:
    print('patient', patient)
    print('group', group)

"""
Explanation: However, the grouping is only an intermediate step; for example, we may want to iterate over each of the patient groups:
End of explanation
"""


cdystonia_grouped.agg(np.mean).head()

"""
Explanation: A common data analysis procedure is the split-apply-combine operation, which groups subsets of data together, applies a function to each of the groups, then recombines them into a new data table.
For example, we may want to aggregate our data with with some function.

<div align="right">*(figure taken from "Python for Data Analysis", p.251)*</div>

We can aggregate in Pandas using the aggregate (or agg, for short) method:
End of explanation
"""


cdystonia_grouped.mean().head()

"""
Explanation: Notice that the treat and sex variables are not included in the aggregation. Since it does not make sense to aggregate non-string variables, these columns are simply ignored by the method.
Some aggregation functions are so common that Pandas has a convenience method for them, such as mean:
End of explanation
"""


cdystonia_grouped.mean().add_suffix('_mean').head()

"""
Explanation: The add_prefix and add_suffix methods can be used to give the columns of the resulting table labels that reflect the transformation:
End of explanation
"""


# Write your answer here

"""
Explanation: Exercise
Use the quantile method to generate the median values of the twstrs variable for each patient.
End of explanation
"""


cdystonia.groupby(['week','site']).mean().head()

"""
Explanation: If we wish, we can easily aggregate according to multiple keys:
End of explanation
"""


normalize = lambda x: (x - x.mean())/x.std()

cdystonia_grouped.transform(normalize).head()

"""
Explanation: Alternately, we can transform the data, using a function of our choice with the transform method:
End of explanation
"""


%timeit cdystonia_grouped['twstrs'].mean().head()

"""
Explanation: It is easy to do column selection within groupby operations, if we are only interested split-apply-combine operations on a subset of columns:
End of explanation
"""


cdystonia_grouped[['twstrs']].mean().head()

"""
Explanation: Or, as a DataFrame:
End of explanation
"""


chunks = dict(list(cdystonia_grouped))

chunks[4]

"""
Explanation: If you simply want to divide your DataFrame into chunks for later use, its easy to convert them into a dict so that they can be easily indexed out as needed:
End of explanation
"""


dict(list(cdystonia.groupby(cdystonia.dtypes, axis=1)))

"""
Explanation: By default, groupby groups by row, but we can specify the axis argument to change this. For example, we can group our columns by dtype this way:
End of explanation
"""


cdystonia2.head(10)

"""
Explanation: Its also possible to group by one or more levels of a hierarchical index. Recall cdystonia2, which we created with a hierarchical index:
End of explanation
"""


cdystonia2.groupby(level='obs', axis=0)['twstrs'].mean()

"""
Explanation: The level argument specifies which level of the index to use for grouping.
End of explanation
"""


def top(df, column, n=5):
    return df.sort_index(by=column, ascending=False)[:n]

"""
Explanation: Apply
We can generalize the split-apply-combine methodology by using apply function. This allows us to invoke any function we wish on a grouped dataset and recombine them into a DataFrame.
The function below takes a DataFrame and a column name, sorts by the column, and takes the n largest values of that column. We can use this with apply to return the largest values from every group in a DataFrame in a single call.
End of explanation
"""


goo  = vessels_merged.groupby('mmsi')

top3segments = vessels_merged.groupby('mmsi').apply(top, column='seg_length', n=3)[['names', 'seg_length']]
top3segments.head(15)

"""
Explanation: To see this in action, consider the vessel transit segments dataset (which we merged with the vessel information to yield segments_merged). Say we wanted to return the 3 longest segments travelled by each ship:
End of explanation
"""


from IPython.core.display import HTML
HTML(filename='../data/titanic.html')

"""
Explanation: Notice that additional arguments for the applied function can be passed via apply after the function name. It assumes that the DataFrame is the first argument.
Exercise
Load the dataset in titanic.xls. It contains data on all the passengers that travelled on the Titanic.
End of explanation
"""


# Write your answer here

"""
Explanation: Women and children first?

Use the groupby method to calculate the proportion of passengers that survived by sex.
Calculate the same proportion, but by class and sex.
Create age categories: children (under 14 years), adolescents (14-20), adult (21-64), and senior(65+), and calculate survival proportions by age category, class and sex.
End of explanation
"""


