<jupyter_start><jupyter_text>Tutorial 2a: Tidy data and split-apply-combine(c) 2017 <NAME>. This work is licensed under a [Creative Commons Attribution License CC-BY 4.0](https://creativecommons.org/licenses/by/4.0/). All code contained herein is licensed under an [MIT license](https://opensource.org/licenses/MIT).<jupyter_code>import numpy as np
import pandas as pd

import bokeh.io
import bokeh.plotting
bokeh.io.output_notebook()<jupyter_output><empty_output><jupyter_text>Note: *The code in this tutorial comes from the tutorial by <NAME>, which can be found [here](http://bebi103.caltech.edu.s3-website-us-east-1.amazonaws.com/2017/tutorials/t2a_tidy_data.html).* IntroductionThe data we are investigating come from <NAME>'s lab. A description of their work on the genetic regulation of sleep can be found on the research page of the lab website. There is a movie of the moving/sleeping larvae similar to the one used to produce the data set we are using in this tutorial. The work based on this data set was published in (Gandhi et al., *Neuron*, 85, 1193–1199, 2015). The genotype dataFirst, we will load the genotype file.<jupyter_code>with open('data/130315_1A_genotypes.txt', 'r') as f:
    for _ in range(30):
        print(next(f), end='')<jupyter_output># Genotype data from the Gandhi, et al. experiment ending March 13, 2013
#
# The experiment was performed in a 96 well plate with zebrafish
# embryos.  Gene sequencing was later used to identify the genotype
# of the fish in each well.  Not all fish could be genotyped.
#
# The mutants being studied have deletions in the gene coding for
# arylalkylamine N-acetyltransferase (aanat), which is a key enzyme
# in the rhythmic production of melatonin.  Melatonin is a hormone
# responsible for regulation of circadian rhythms.  It is often taken
# as a drug to treat sleep disorders.  The goal of this study is to
# investigate the effects of aanat deletion on sleep pattern in
# 5+ day old zebrafish larvae.
#
# Each column lists the wells corresponding to the genotype of
# each fish.  If a number is missing (between 1 and 96), the 
# genotype of that fish is not known.
# 
# These data were kindly provided by <NAME> and <NAME>
# from <NAME>'s lab.  They were part of the paper Gandh[...]<jupyter_text>Each column in this file contains a list of wells in the 96-well plate corresponding to each genotype. Apparently, the columns are tab delimited. We can parse this using the `delimiter` keyword argument of `pd.read_csv()`. When specifying delimiters, tabs are denoted as `\t`. There are two header rows, one starting with `Genotype 1` and the other one starting with `WH 17`. The first one is really redundant, so it can be skipped. Nevertheless, we will read them both in as headers using the `header` kwarg of `pd.read_csv()`.<jupyter_code># Load in the genotype file, call it df_gt for genotype DataFrame
df_gt = pd.read_csv('data/130315_1A_genotypes.txt',
                    delimiter='\t', 
                    comment='#',
                    header=[0 ,1])

# Let's look at it
df_gt<jupyter_output><empty_output><jupyter_text>Now we have two header rows for the columns. If we look at the column names, we see that they are a `MultiIndex` instance.<jupyter_code>df_gt.columns<jupyter_output><empty_output><jupyter_text>This is part of Pandas's cool multi-indexing functionality, which is not comething necessary if we have tidy data. However, the messay data with multi-indexing may result in performance boosts for accessing data. This can be important for very large data sets. But in general, tidy data are easier to conceptualize and syntactically much simpler to access.In this case, we do not need the multi-indexing, and zero-level indexing (`Genotype1`, etc.), so we will leave just level one index. This can be accomplished using the `get_level_values` method of Pandas `MultiIndex` objects.<jupyter_code># Reset the columns to be the second level of indexing
df_gt.columns = df_gt.columns.get_level_values(1)

# Check out the new columns
df_gt.columns

# Let's actually clean up the column names
df_gt.columns = ['wt', 'het', 'mut']<jupyter_output><empty_output><jupyter_text>The above approach to renaming is a shortcut to using `rename()` method that has previously been encountered. It could actually have been done before messing around with the multi-indexing, but why not to learn something new? We need to know how to tidy data sets in this format anyways. Tidying the genotype dataAs they are, the data are not tidy. For tidy data, we would have two columns, `location`, which gives the well number for each larva, and `genotype`, which is `wt` for wild type, `het` for heterozygote, or `mut` for mutant.A useful data tidying tool is the `pd.melt()` function. For this simple data set, it takes the column headings and makes them into a column (with repeated entries) and puts the data accordingly in the correct order. We just need to specify the name of the "variable" column, in this case the genotype, and the name of the "value" column, in this case, the fish ID.<jupyter_code># Tidy the DataFrame
df_gt = pd.melt(df_gt, var_name='genotype', value_name='location')

# Take a look
df_gt<jupyter_output><empty_output><jupyter_text>We still have some unhelpful `NaN` entries, so we drop them using `dropna()` method of `DataFrame`s.<jupyter_code># Drop NaNs
df_gt = df_gt.dropna()

# Take a look
df_gt<jupyter_output><empty_output><jupyter_text>We can see that indices have some skips due to using the NaN's. This is not really of concern, since we will not use them. We can use the `reset_index()` to reset the indexing of the `DataFrame`.<jupyter_code># Reset the indexing of the DataFrame
df_gt = df_gt.reset_index(drop=True)<jupyter_output><empty_output><jupyter_text>We now have a tidy `DataFrame`. However, the fish numbers are floates. This happens becuase NaN is a float, so the column got converted to floats. We can reset the data type of the column to `int`s.<jupyter_code># Set data types to be integers
df_gt.loc[:, 'location'] = df_gt.loc[:, 'location'].astype(int)<jupyter_output><empty_output><jupyter_text>We now have a beautiful, tidy `DataFrame` of genotypes. Apparently, the `pd.melt()` may be the most useful function in tidying messy data sets. Now let's take a look at the behavioral data. The behavioral dataWe will first take a look at the behavioral data set. It has been preprocessed by <NAME> (thanks, Justin!)  to some degree, and we will work with the original raw data later on. Let's look at the contents of the data file.<jupyter_code>with open('data/130315_1A_aanat2.csv', 'r') as f:
    for _ in range(40):
        print(next(f), end='')<jupyter_output># Lightly processed data from VideoTracker from the Gandhi, et al. experiment
# concluding on March 15, 2013.
#
# The experiment was performed in a 96 well plate with zebrafish
# embryos. Gene sequencing was later used to identify the genotype
# of the fish in each well. Not all fish could be genotyped.
#
# The mutants being studied have deletions in the gene coding for
# arylalkylamine N-acetyltransferase (aanat), which is a key enzyme
# in the rhythmic production of melatonin. Melatonin is a hormone
# responsible for regulation of circadian rhythms. It is often taken
# as a drug to treat sleep disorders. The goal of this study is to
# investigate the effects of aanat deletion on sleep pattern in
# 5+ day old zebrafish embryos.
#
# Activity is defined as the number of seconds over the one-minute interval 
# in which a given larva was moving.
#
# The column 'zeit' contains the so-called Zeitgeber time in units of hours.
# The Zeitgeber time is zero when the lights come on on the first [...]<jupyter_text>The provided data set is tidy. Each row corresponds to the measured seconds of activity of a single fish for a single minute. Other information includes the time of the measurement and the day in the life of the fish. The Zeitgeber time seem to somehow give us a "reference time" for the experiment (`Zeit` ist "time" in Deutsch. Alles klar?)Let's load the data set<jupyter_code>df = pd.read_csv('data/130315_1A_aanat2.csv', comment='#')

# Let's take a look
df.head()<jupyter_output><empty_output><jupyter_text>Adding the genotype informationFirst of all, we can try to sort out data by including the genotype information for the fish at each location. The `DataFrame` genotyp has two columns, `location` and `genotype`. The `location` column is shared with our behavior `DataFrame` so we can use `pd.merge()` to add the genotype information as a new column.<jupyter_code>df = pd.merge(df, df_gt)

# Take a look
df.head()<jupyter_output><empty_output><jupyter_text>Light or dark?We can attempt to annotate our `DataFrame` to indicate if it is light or dark. Then we can easily separate the diurnal and nocturnal activity. So let's make a column `light` which will have the entry `True` if the lights are on, and `False` if the lights are off. Referring again to the comments in the header of the data set, the light come on at 9AM and turn off at 11PM. We can use the `time` column to determine this.We should first look at the time column and determine its data type.<jupyter_code>df['time'].dtype<jupyter_output><empty_output><jupyter_text>The result `dtype('O')` means that the data type is "object", which is a generic catch-all for unknown data types. However, we know that the `time` column are actual clock times. We can tell this to Pandas and unleash its data processing power on our data set. We do this using the `pd.to_datetime()` function.<jupyter_code>df['time'] = pd.to_datetime(df['time'])

# What is the data type now?
df['time'].dtype<jupyter_output><empty_output><jupyter_text>The data type is ` * If data are taken on the greater than nanosecond frequency, Pandas's datetime utility will not help. In such case, we would simply use a number (say femtoseconds) as the time variable.* The time zero can not be more than about a thousand years ago. This can come up in geology or fields like that. We can do lots of time-based things with this datetime format, like extraction of the time from the datetime. We use `.dt.time` to do this.<jupyter_code>df['time'].dt.time.head()<jupyter_output><empty_output><jupyter_text>We can also compare times. For example, let's ask if the time is after 9AM.<jupyter_code>(df['time'].dt.time > pd.to_datetime('9:00:00').time()).head()<jupyter_output><empty_output><jupyter_text>We had to convert the string `9:oo:oo` to datetime and then extract the time using the `time()` method. Now we can just make a column of booleans reporting whether the time is greater than 9:00:00 and less than 23:00:00.<jupyter_code>df['light'] = ( (df['time'].dt.time >= pd.to_datetime('9:00:00').time())
               & (df['time'].dt.time < pd.to_datetime('23:00:00').time()))

# Take a look
df.head()<jupyter_output><empty_output><jupyter_text>We make a quick plot of the Zeitgeber time versus light. We only need to do this for a single location (single well).<jupyter_code>p = bokeh.plotting.figure(plot_height=259,
                          plot_width=700,
                          x_axis_label='Zeitgeber time (hours)',
                          y_axis_label='light')

inds = df['location'] == 1
p.circle(df.loc[inds, 'zeit'], df.loc[inds, 'light'])

bokeh.io.show(p)<jupyter_output><empty_output><jupyter_text>Apparently, the light switch on and off as it should. Split-apply-combineLet's now compute some more things. For instance, let's say we want to compute the average per-minute activity of each fish over the course of the experiment. To do this, we must:1. **Split** the data set according to the `location` field, i. e. split it up so we have a separate data set for each fish.2. **Apply** an averaging function to the activity in these split data sets.3. **Combine** the results of these averages on the split data set into a new, summary data set that contains the locations and means for each location.The strategy we want is a **split-apply-combine** strategy. This idea was put forward by <NAME> in [this paper](https://www.jstatsoft.org/article/view/v040i01) and is used quite often:Split data according to some criterion (maybe according to genotype). Apply some function to split-up data. Combine the results into a new data frame. If the data are tidy, it makes a lot of sense. Choose ther column by which we want to split the data. All rows with like entries in the splitting column are then grouped into a new data set. Pandas's split-apply-combine operations are achieved using the `groupby()` method. It can be thought of as a the splitting part. We can then apply functions to the resulting `DataFrameGroupBy` object.<jupyter_code>gb = df.groupby('location')

# Let's look at it
gb<jupyter_output><empty_output><jupyter_text>There is really nothing to see in the `DataFrameGroupBy` object. It is a split-up data set. We will now apply the `np.mean()` function to this `DataFrameGroupBy` object. To apply a funciton, we use the `agg()` method. There also is an `apply()` method that is more generic at the cost of performance. For computing summary statsitics, we are aggregating, hence the use of the `agg()` method, and use of this function gives better performance.)<jupyter_code># Compute the mean of the data grouped by location
df_mean = gb.agg(np.mean)

# Take a look
df_mean.head()<jupyter_output><empty_output><jupyter_text>We got back a `DataFrame`. So the `agg()` method of the `DataFrameGroupBy` object did both the **apply** and **combine** steps. This makes sense - we would never want to apply, and then not combine it into something useful.So in Pandas, split-apply-combine to compute summary statistics is achieved by first doing a `groupby()` operation to get a `DataFrameGroupBy` object and then an `agg()` operation on it. The result is a `DataFrame`.The index of the `DataFrame` is now called `location`, as the original `DataFrame` was grouped by the `'location'`. It then applied the `np.mean()` function to every column of the original `DataFrame`. However, we only wanted it to be applied to `activity`. Let's try the operation again, this time slicing out the location column before doing the `agg()` step.<jupyter_code>df_mean = gb['activity'].agg(np.mean)

# Take a look
df_mean.head()<jupyter_output><empty_output><jupyter_text>We now have a Pandas `Series`, where the index is `location`. Because we like to work with datasets with Boolean indexing, we would prefer a `DataFrame` with two columns, `location` and `activity`. We can get this using the `reset_index()` method.<jupyter_code>df_mean = df_mean.reset_index()

# Take a look
df_mean.head()<jupyter_output><empty_output><jupyter_text>The above `DataFrame` gives us the average seconds of activity per second of each fish.Quick note: Some of the functions we apply are very common, like the mean. For the functions of this type, Pandas often has built-in methods. Let's try using a `DataFrameGroupBy` object's `mean()` method.<jupyter_code>df_mean = gb['activity'].mean().reset_index()

# Take a look
df_mean.head()<jupyter_output><empty_output><jupyter_text>We got exactly the same result! Grouping by multiple columnsWhat we may *really* want is the mean activity for each genotype for each day and for each night. For this, we want to split the data set according to the `genotype`, `day`, and `light` columns. To do this, we simply pass a list of columns as an the argument to the `groupby()` method.<jupyter_code># Group by three columns
gb = df.groupby(['genotype', 'day', 'light'])

# Apply the mean and reset index
df_mean = gb['activity'].mean().reset_index()

# Take a look
df_mean<jupyter_output><empty_output><jupyter_text>Now we can compute other summary statistics, such as the standard deviation. The `agg()` method allows passing multiple functions.<jupyter_code>df_summary = gb['activity'].agg([np.mean, np.std])

# Take a look
df_summary<jupyter_output><empty_output><jupyter_text>Columns are automatically named with the functions we chose to compute. We can also rename these columns and `reset_index()` them to get the results we want.<jupyter_code>df_summary = (gb['activity'].agg([np.mean, np.std])
                            .reset_index()
                            .rename(columns={'mean': 'mean activity',
                                             'std': 'std activity'}))

# Take a look
df_summary<jupyter_output><empty_output><jupyter_text>The above example show that we can rapidly compute and then compare summary statistics, even for the large data sets with split-apply-combine approach ResamplingLet's look at a single plot of the activity of a fish.<jupyter_code>p = bokeh.plotting.figure(plot_height=250,
                          plot_width=700,
                          x_axis_label='Zeitgeber time (hours)',
                          y_axis_label='activity (sec/min)')

inds = df['location'] == 1
p.line(df.loc[inds, 'zeit'], df.loc[inds, 'activity'], line_join='bevel')

bokeh.io.show(p)<jupyter_output><empty_output><jupyter_text>We used thr `line_join=bevel` argument. This ensures that we do not get weird point artifacts in line plots with highly variable signal. This plot might look a little nicer if we colored it differently for night versus day. So we will loop through each day, plot a line with light color when the lights were on and a line with dark color when the lights were of ("days" begin and end at lights on, not a midnight). For coloring, we will use one of Bokeh's built-in color palettes, based on the work of Cynthia Brewer. Since we might need to do this over and over again, let's write a fuction for this.<jupyter_code>def plot_trace(df, p=None, activity='activity', y_axis_label='activity (sec/min)',
               colors=bokeh.palettes.brewer['Paired'][4]):
    """Plot a time trace of fish activity with differently colored
    light and dark periods.
    """
    
    if p is None:
        p = bokeh.plotting.figure(plot_height=250,
                                  plot_width=700,
                                  x_axis_label='Zeutgeber time (hours)',
                                  y_axis_label=y_axis_label)
        
        for day in df['day'].unique():
            # Light color lines for day
            inds = (df['day'] == day) & (df['light'])
            p.line(df.loc[inds, 'zeit'],
                   df.loc[inds, activity],
                   line_join='bevel',
                   color=colors[0])
            
            # Dark lines for night
            inds = (df['day'] == day) & (~df['light'])
            p.line(df.loc[inds, 'zeit'],
                   df.loc[inds, activity],
                   line_join='bevel',
                   color=colors[1])
            
        return p<jupyter_output><empty_output><jupyter_text>Let's use the function.<jupyter_code>p = plot_trace(df.loc[(df['location'] == 1), :])
bokeh.io.show(p)<jupyter_output><empty_output><jupyter_text>The data are quite noisy, as the fish are sometimes quiescent even when they are awake and darting around the well. We can consider a lower sampling frequency for the activity in such case. We can sum the number of seconds of movement (we call activity) over a longer time window than the default one minute - say ten minutes. We could do it by hand, looping over each location and then performing the sum over each ten minute interval. It would be quite lengthy, though. Instead, we will do it using built-in Pandas functions. Using rolling aggregationWe first need to group data set by `location`, as we want to resample the activity for each fish. The `DataFrameGroupBy` object has a `rolling()` method that will create an object over which you may compute rolling averages. Before we build this object, it is convenient to have the `DataFrame` sorted by `location` and then by `zeit_ind`, which is the index of the time point. This is a nice ordering, as each fish has their time trace in order, and the fishes are in order too. We can accomplish this sorting using the `DataFrame`s `sort_values()` along with the `by` kwarg.<jupyter_code>df = df.sort_values(by=['location', 'zeit_ind'])

# Take a look
df.head()<jupyter_output><empty_output><jupyter_text>Now when data is in order, let's create the `DataFrameGroupBy` object and then generate a `RollingGroupby` object. We only need `zeit_ind` and `activity` columns. To make the `RollingGroupby` object, we use the `on` keyworkd argument to indicate that we are going to use the `zeit_ind` column instead of the index to determine the window.<jupyter_code># Create GroupBy object
gb = df.groupby('location')['zeit_ind', 'activity']

# Make a RollingGroupby with window size of 10
rolling = gb.rolling(window=10, on='zeit_ind')

# Look at rolling object
rolling<jupyter_output><empty_output><jupyter_text>Having the `RollingGroupby` object, we can use its `sum()` method to compute a rolling sum.<jupyter_code>df_rolling = rolling.sum()

# Take a look
df_rolling.head(n=20)<jupyter_output><empty_output><jupyter_text>We have two indices, and we want to reset them so that they are columns in the `DataFrame`. We use the `level=0` kwarg to indicate that we only want to location to be reindexed.<jupyter_code>df_rolling = df_rolling.reset_index(level=0)

# Take a look
df_rolling.head(n=20)<jupyter_output><empty_output><jupyter_text>The first nine entries are `NaN`. This is because we cannot compute anything for a window of length 10 when we have less than 10 entries.We can now insert the rolling summed activity into the original `DataFrame` because we have taken care to have it properly sorted.<jupyter_code>df['rolling activity'] = df_rolling['activity']

# Take a look
df.head(n=20)<jupyter_output><empty_output><jupyter_text>Now let's plot the rolling activity against the Zeitgeber time.<jupyter_code>p = plot_trace(df.loc[df['location'] == 1, :],
               activity = 'rolling activity',
               y_axis_label = 'activity (sec/ 10 min)')

bokeh.io.show(p)<jupyter_output><empty_output><jupyter_text>The plot with 10-min windows definitely has less noise, but there is still an entry for each time point. This is because the `rolling()` object allows for overlapping rolling windows. We therefore want to keep only every 10th entry for each fish in the `DataFrame`. We can build a new, resampled `DataFrame` by taking every tenth time point using `zeit_ind[9::10]`. We can check and see if a given `zeit_ind` is in the list of `zeit_ind`s we want to keep using `isin()` method.<jupyter_code># Get all zeit_inds, sorted
zeit_inds = np.sort(df['zeit_ind'].unique())

# Select every tenth, starting at tenth
zeit_inds = zeit_inds[9::10]

# Keep all entries matching the zeit_inds in the DataFrame
df_resampled = df.loc[df['zeit_ind'].isin(zeit_inds), :]

# Drop orignial activity column
del df_resampled['activity']

# Rename the rolling activity column to activity
df_resampled = df_resampled.rename(columns={'rolling activity': 'activity'})

# Take a look
df_resampled.head(n=20)<jupyter_output><empty_output><jupyter_text>This looks good! Now let's make a plot of the trace for the fish in location 1.<jupyter_code>p = plot_trace(df_resampled.loc[df['location']==1, :],
               activity='activity',
               y_axis_label='activity (sec/10 min)')

bokeh.io.show(p)<jupyter_output><empty_output><jupyter_text>Resampling using the resample methodThere are several problems with this approach:1. We might not want to resample past boundaries between light and dark (because then we get time points that contain both light and dark moments).2. If the time points are not evenly sampled, by using a fixed window of 10 time points, we may be sampling time windows of varying lengths.For the second point, we can check the difference between respective time point to see if the data are evenly sampled.<jupyter_code># Compute difference between time points in units of minutes
delta_t = np.diff(df.loc[df['location']==1, 'time']) / 60 / 1e3

# Plot the difference
p = bokeh.plotting.figure(plot_height=200,
                          plot_width=500,
                          x_axis_label='obseration number',
                          y_axis_label='Δt (min)')
p.circle(np.arange(len(delta_t)), delta_t)
bokeh.io.show(p)<jupyter_output><empty_output><jupyter_text>HA! We can see that there are two measurements which are almost 5 minutes apart between time points. This is because the experimenter, in this case, had to check the instrument and momentarily pause acquisition. Here is a big statement:         Always validate the dataset. We need to know about the possible discrepancies in the data **before** performing data analysis. If we are assuming the equal time steps, and they are not equal, we need to know that before doing the analyisis.Now, we can alleviate these problems by using Pandas's built-in `resample()` method, which works on `DataFrameGroupBy` objects. It uses datetime columns to resample data to different time windows. Like we created a `RollingGroupby` object before, we can make an analogous resampling object.<jupyter_code>gb = df.groupby('location')
resampler = gb.resample('10min', on='time')

# Take a look
resampler<jupyter_output><empty_output><jupyter_text>We provided the `resample()` method with the argument `10min` which tells it we want to resample to ten-minute intervals and the `on='time'` kwarg says to use the `time` column for resampling. For instructions on the strings for specifying the resampling frequency, we can look [here](http://pandas.pydata.org/pandas-docs/stable/timeseries.htmloffset-aliases).Now that we know that we are resampling, we can compute the sums. Let's start by summing, which is our goal; we want to know how many seconds of activity the fish had in 10 minutes.<jupyter_code>resampler.sum()['activity'].head()<jupyter_output><empty_output><jupyter_text>Apparently, the resampler uses exact 10-minute intervals; that is, it gives the result at 18:30, 18:40, etc. This alleviates the problem of uneven sampling, in that it only aggregates over real time. However, there are some possible problems with this.1. Forcing of the time intervals on the whole-ten-minute intervals causes some clipping, as the minutes measured start and end on fractional minutes. *It is important to understand how this rounding works, and what possible effects it might have on the subsequent analyses.*2. Some time intervals have fewer observations than others. If we are looking at the mean activity per minute on those intervals, things are OK (though some are undersampled compared to others), but if we sum, some intevals will show lower activity when the fish are just as active. To obtain seconds of activity per 10 minutes, we should adjust the sums to scale with the number of measurements contained in the interval.3. Because the lights on and lights off times are 9:00:00 and 23:00:00, we have dealt with the problem of resampling across lights on/of events. If this were not the case, this strategy would not deal with that.To fix the problem (2), we can take the mean and then scale it according to how many observations there were using `count()` method of the resampler.<jupyter_code># Resample with summing
activity_resampled = resampler.sum()['activity']

# Rescale sum with count
activity_resampled *= 10 / resampler.count()['activity']

# Reset the index
activity_resampled = activity_resampled.reset_index()

# Take a look
activity_resampled.head()<jupyter_output><empty_output><jupyter_text>Activity is now approximately resampled to the ten minutes intervals properly. "Approximately", because there is still issue (1) above, that we are focing whole-minute time stamps, resulting in some errors around the time of lights swithching on and off. However, not much can be done about this, as the experiment was not started and ended on whole-minute intervals, whereas lights on and off are exactly on a minute. While we have the `location`, `time` and `activity` times properly resampled, we need to get the other columns as well. We can use the `first()` method of the resampler to get those entries (since we do not want to sum over them). The `first()` method five the first entry in the resampled interval.<jupyter_code># Get a new DatFrame, resampled without summing.
df_resampled = resampler.first()

# Resampling was done on activity and time columns; delete them
del df_resampled['time']
del df_resampled['location']

# Reset index so resampled location and time indices become columns
df_resampled = df_resampled.reset_index()

# Add in the properly resampled activity
df_resampled['activity'] = activity_resampled['activity']

# Take a look 
df_resampled.head()<jupyter_output><empty_output><jupyter_text>However, `rolling activity` column is a remnant from the earlier analysis, so we should delete it. Also, the `zeit` column no longer represents the `time` column, since that has been resampled. When we applied the `first()` method of the resampler, it did not know that `zeit` is connected to the time. So, we need to re-compute the `zeit` columns.<jupyter_code># Set the Zeitgeber time
zeitgeber_0 = pd.to_datetime('2013-03-16 9:00:00')
df_resampled['zeit'] = (df_resampled['time'] - zeitgeber_0).dt.total_seconds() / 3600

# Delete rolling activity
del df_resampled['rolling activity']

# Take a look
df_resampled.head()<jupyter_output><empty_output><jupyter_text>We can now make a plot of these resampled data.<jupyter_code>p = plot_trace(df_resampled.loc[df_resampled['location']==1, :],
               activity='activity',
               y_axis_label='activity (sec/10 min)')
bokeh.io.show(p)<jupyter_output><empty_output><jupyter_text>This seems to be similar to what we have done before, but we took better care of missing time points. Saving work and pipelineWhen wranging data as we have done, it is important to have code and documentation that will perform *exactly* the same steps as we did to organize the data. This is true of the entire data analysis process, from notes on acquisition to validation to wrangling to inference and reporting. Jupyter notebooks can help, but it is often useful to have functions to do the analysis pipeline we would typically do. For example, we could write a function to generate the resampled `DataFrame` we just generated from the original data.<jupyter_code>def load_and_resample(activity_file, genotype_file, zeitgeber_0, resample_r=None):
    """
    Load and resample activity data.
    Assumes genotype file has columns corresponding
    to wild type, heterozygote, and mutant.
    """
    
    # Load in the genotype file, call it df_gt for genotype DataFrame
    df_gt = pd.read_csv(genotype_file,
                        delimiter='\t',
                        comment='#',
                        header=[0, 1])
    
    # Rename columns
    df_gt.columns = ['wt', 'het', 'mut']
    
    # Melt to tidy
    df_gt = pd.melt(df_gt, var_name='genotype', value_name='location').dropna()
    
    # Reset index
    df_gt = df_gt.reset_index(drop=True)
    
    # Integer location names
    df_gt.loc[:, 'location'] = df_gt.loc[:, 'location'].astype(int)
    
    # Read in activity data
    df = pd.read_csv(activity_file, comment='#')
    
    # Merge with genotype data
    df = pd.merge(df, df_gt)
    
    # Convert time to datetime
    df['time'] = pd.to_datetime(df['time'])
    
    # Column for light or dark
    df['light'] = (  (df['time'].dt.time >=pd.to_datetime('9:00:00').time())
                   & (df['time'].dt.time < pd.to_datetime('23:00:00').time()))
    
    if resample_r is None:
        df_resampled = df
    else:
        # Group by location
        gb = df.groupby('location')
        
        # Make resampler
        resampler = gb.resample(resample_r, on='time')
        
        # Resample with summing
        activity_resampled = resampler.sum()['activity']
        
        # Rescale sum with count
        activity_resampled *= 10 / resampler.count()['activity']
        
        # Reset the index
        activity_resampled = activity_resampled.reset_index()
        
        # Get a new DataFrame, resampled without summing.
        df_resampled = resampler.first()
        
        # Resampling happened on activity and time columns; delete them
        del df_resampled['time']
        del df_resampled['location']
        
        # Reset index so resampled location and time indices become columns
        df_resampled = df_resampled.reset_index()
        
        # Add in the properly resampled activity
        df_resampled['activity'] = activity_resampled['activity']
        
    # Set the Zeitgeber time
    zeitgeber_0 = pd.to_datetime(zeitgeber_0)
    df_resampled['zeit'] = (df_resampled['time'] - zeitgeber_0).dt.total_seconds() / 3600
    
    return df_resampled<jupyter_output><empty_output><jupyter_text>This function makes strong assumptions about the structure of the data files, which need to be verified first in order to use the function. Also ,it should have a much more descriptive doc string explaining these things.Now that we have a function, we can use it to load in the original data sets and then generate the resampled tidy data frame.<jupyter_code>df = load_and_resample('data/130315_1A_aanat2.csv', 
                       'data/130315_1A_genotypes.txt', 
                       '2013-03-16 9:00:00',
                       resample_r='10min')
# Take a look
df.head()<jupyter_output><empty_output><jupyter_text>To save the processed `DataFrame` to a CSV file, we can use `to_csv()` method of a `DataFrame`.<jupyter_code>df.to_csv('data/130315_1A_aanat2_resampled.csv', index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08425246008859291, 'fraction_numerical': 0.00984035437167405, 'mean_word_length': 4.160785517029764, 'pattern_counts': {'":': 0, '<': 197, '<?xml version=': 0, '>': 198, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '47', 'dir': 'jupyter-structured-clean-dedup', 'id': '326775', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 10061, 'n_tokens_neox': 9529, 'n_words': 4206, 'path': 't2a_tidy_data.ipynb', 'repo_name': 'MiroGasparek/DataAnalysis_intro', 'size': '33682'}