<jupyter_start><jupyter_text>Beachwatch Logistic RegressionsThis analysis merges the Beachwatch and Environmental data packages, then does logistic regressions on bacteria counts versus  rainfal and riverflow.<jupyter_code>import matplotlib.pyplot as plt 
import metapack as mp
import pandas as pd
import numpy as np
import seaborn as sns
sns.set(color_codes=True)
%matplotlib inline<jupyter_output><empty_output><jupyter_text>The beachwatch datasets we will be  using is derived from the original data, and includes an extra set of feature columns for measurement groups and the quantile of each measurement relative to its group. Refer to the [package README file for details about these aditional columns.](https://data.sandiegodata.org/dataset/sandiegodata-org-beachwatch)In the schema below, note the columns that have '_group_' in the name, which are the group means and quantiles, and the  columsn that have comparison names, such as '_gt_' ( greater than )  '_gte_' ( greater than or equal ) or '_lte_' ( less than or equal.<jupyter_code>import metapack as mp
pkg = mp.open_package('http://library.metatab.org/sandiegodata.org-beachwatch-4.zip')
pkg.resource('beachwatch')
beachwatch =    pkg.resource('beachwatch').read_csv(parse_dates=True)
measure_codes = pkg.resource('measure_codes').read_csv(parse_dates=True)
stations =      pkg.resource('stations').read_csv(parse_dates=True)
beachwatch.head()
beachwatch.lresult_gt_lmedian.value_counts(dropna=False)<jupyter_output><empty_output><jupyter_text>Most of the time, when a result is creater than the mean, it is also greater than the median, and when it is less than the mean, it is also less than the median.<jupyter_code># Proportion of gt_mean and gt_median that are different. 
_1 = beachwatch[['lresult_gt_lmedian','lresult_gt_lmedian']].dropna()

len(_1[beachwatch.lresult_gt_lmedian != beachwatch.lresult_gt_lmean])/len(_1)<jupyter_output><empty_output><jupyter_text>Merge in environmental dataNow we can merge in the rain, river and tides dataset.<jupyter_code>env_pkg = mp.open_package('http://library.metatab.org/sandiegodata.org-tides_rain_river-6.zip')
env_pkg.resource('tides_river_rain')
trr = env_pkg.resource('tides_river_rain').read_csv(parse_dates=True).set_index('date')
trr.head()<jupyter_output><empty_output><jupyter_text>We're going to use ``recent_precipt``, a three day rolling average, rather than ``dailyprecip`` because the bacteria samples are sporadic, and we don't want to use only measurements that were taken on the same day as rainfall, especially if a recent rainfall will increase the bacteria counts.<jupyter_code>trr.loc['2007':'2008'].recent_precip.plot()
trr.groupby(by=[trr.index.month, trr.index.year]).sum().recent_precip.plot()
bw_env = beachwatch.set_index('sampledate').join(trr)
bw_env['high_flow'] =  (bw_env.recent_total_discharge > bw_env.recent_total_discharge.median()).astype(int)
bw_env['ldischarge_rate'] = np.log(bw_env.discharge_rate)

bw_env.head()<jupyter_output><empty_output><jupyter_text>Sanity Checking and ExplorationLet's look at a few of the variables to see if they are sensible.<jupyter_code># The values are 1 or 0, so a sum is equivalent to a count. 
sum(bw_env.recent_precip > 0.0), sum(bw_env.recent_precip > 0.0)/len(bw_env)<jupyter_output><empty_output><jupyter_text>The mean discharge rate is a lot larger for the above mean results compared to the below mean results.<jupyter_code>bw_env.groupby('lresult_gt_lmean').mean().discharge_rate.reset_index()<jupyter_output><empty_output><jupyter_text>However, it looks like all of that difference is due to a small number of high-flow events. Both plots have very long tails.<jupyter_code>ax = sns.violinplot(x=bw_env.lresult_gt_lmedian, y=bw_env.discharge_rate)<jupyter_output><empty_output><jupyter_text>Looking at the plot of the log discharge reveals the patterns better.<jupyter_code>ax = sns.violinplot(x=bw_env.lresult_gt_lmedian, y=bw_env.ldischarge_rate)
bw_env.groupby('lresult_gt_lmean').discharge_rate.agg(['mean','std','median'])
bw_env.groupby('lresult_gt_lmean').mean().recent_precip
pd.crosstab(bw_env.lresult_gt_lmedian, bw_env.recent_precip > 0.01 ).plot(kind='bar')<jupyter_output><empty_output><jupyter_text>Logistic RegressionsLogistic regression produces a model for estimating the probablity of a binary dependent varaible. We will use this to determine the probability that the bacteria count result is above the group median given the recent ( 3 day rolling average) precipitation. We'll focus on measure code 24, which is the largest group of Enterroccous measurements using Enterlert. However, before running the regression, we'll set up a linear regression as an example to explain how the logistic regression works. The goal of the logistic regression is to create a curve of probabilities that the bacteria measurement will be greater than the median given the precipitation. So, if we know the preciptation, we can estimate the probability that the bacteria measurement will be high. First we'll motivate the process using the simipler method of linear regression. First, prepare a simple dataframe for the regression. Note that the y variable is binary, either 0, for not above the median, and 1 if the measurement is above the median. The X value, precipitation, is continuous.<jupyter_code>_1 = bw_env[(bw_env.measure_code==24)]

_2 = pd.DataFrame()

_2['y'] = _1.lresult_gt_lmedian 
_2['X'] = _1.recent_precip.fillna(0)

_2.sort_values('X', inplace=True)

lr = _2

plt.scatter(lr['X'],lr['y'])<jupyter_output><empty_output><jupyter_text>Using Linear RegressionTo better explain how logistic regression works, we'll start by trying to solve the problem with linear regression. Scikit learn is just one of the many packages that provides linear regression. Here we will extract the regression coefficients ( they are stored in ``linreg`` ), and then, using ``linreg`` create a new prediction of the ``y`` variable from the ``X`` variable. Note that the ``X`` variable is referenced with an array; the fit and predict functions expect dataframes of 1 or more columns, not Series.<jupyter_code># fit a linear regression model
from sklearn.linear_model import LinearRegression
linreg = LinearRegression()

linreg.fit(lr[['X']], lr['y'])

lr['predict_lin'] = linreg.predict(lr[['X']])<jupyter_output><empty_output><jupyter_text>Now that we have a prediction, we can plot the line and the points. The blue points are the values for precipitation,  divided into those with  rows for measurements above the median, and those with measurements below. The red line is the regression through those points, and the green is a decision curve that transitions from 0 to 1 at the point that the regression line has a y value of 0.5, the point where the expected probability of the measurement being above the median is more likely than not. Note that the red regression line extends above y=1. Since our prediction is supposed to be either 0 or 1, this is nonsense. Logistic regression instead tries to fit a logictic curve, which is 0 at x=0 and 1 when x is large, which is much more sensible. Also note that there are a lot of points for y=0 to the right of the decision line. These are points are errors,  for which the regression would predict that the result would be above the median, when in fact it was below.<jupyter_code>plt.scatter(lr['X'],lr['y'])
plt.plot(lr['X'],lr['predict_lin'], color='red')  

# Here is the decision line -- where the prediction linear line
# crosses .5
plt.plot(lr['X'], np.where(_2['predict_lin'] >= 0.5, 1, 0), color='green')
print(linreg.coef_, linreg.intercept_)
plt.xlabel('preciptation')
plt.ylabel('> median')<jupyter_output>[0.33075766] 0.24598895530104306<jupyter_text>Logistic Regression using SklearnMuch of the sklearn examples come from:http://nbviewer.jupyter.org/github/justmarkham/DAT8/blob/master/notebooks/12_logistic_regression.ipynbHere we tackle a similar problem to the linear regression process, but instead use logistic regression. First we will use scikit's LogisticRegression to plot the class predictions curve, then the predictied prbabilities curve.<jupyter_code># fit a logistic regression model and store the class predictions
from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression(C=1e9)

logreg.fit(lr[['X']], lr['y'])
lr['predict_log'] = logreg.predict(lr[['X']])<jupyter_output><empty_output><jupyter_text>Here are the class predictions, which transitions from 0 to 1, as in the linear regression example.<jupyter_code># plot the class predictions

plt.scatter(lr['X'],lr['y']) # The blue dots

plt.plot(lr['X'],lr['predict_log'], color='red') # The red line1
plt.xlabel('preciptation')
plt.ylabel('> median');<jupyter_output><empty_output><jupyter_text>While analysis can use the 0 or 1 class prediction curve, which always predicts either 0 or 1, it may be more realistic to generate probabilities that the value is greater than the median.  The predicted probabilities are usually an "S curve" -- also known as the logistic curve, from which the regression gets its name --  with a left tail that asymptotically approches 0, and 1 on the right. The curve for this dataset is a bit wacked, since it never predicts 0, regardless of how low the precipitation is. This suggests that there is some other factor that we should include in the model. However, the curve is 1 at higher precipitation values, suggesting that for higher precipitation, the other unknown factor is not so important.<jupyter_code># store the predicted probabilites of class 1
lr['predict_prob'] = logreg.predict_proba(lr[['X']])[:, 1]

# plot the predicted probabilities
plt.scatter(lr['X'],lr['y'])

plt.plot(lr['X'],lr['predict_prob'], color='red')
plt.xlabel('preciptation')
plt.ylabel('> median');
def sklearn_lr(y,X):
    from sklearn.linear_model import LogisticRegression
    logreg = LogisticRegression(C=1e9)

    lr = pd.DataFrame()
   
    lr['X'] = X
    lr['y'] = y
    
    lr.sort_values('X', inplace=True)
    
    logreg.fit(lr[['X']], lr['y'])
    lr['predict_log'] = logreg.predict(lr[['X']])
    
    plt.scatter(lr['X'],lr['y'], color='blue') # The blue dots

    plt.plot(lr['X'],lr['predict_log'], color='green') # The red line1
    plt.xlabel('preciptation')
    plt.ylabel('> median');
    
    # store the predicted probabilites of class 1
    lr['predict_prob'] = logreg.predict_proba(lr[['X']])[:, 1]

    # plot the predicted probabilities
    plt.scatter(lr['X'],lr['y'])

    plt.plot(lr['X'],lr['predict_prob'], color='red')
    plt.xlabel('preciptation')
    plt.ylabel('> median');
    

sklearn_lr(lr.y,lr.X)<jupyter_output><empty_output><jupyter_text>Here are the coefficients, which we can compare to those of the statsmodel code in the next section.<jupyter_code>"X coef={}, intercept={}".format(logreg.coef_[0], logreg.intercept_)<jupyter_output><empty_output><jupyter_text>Logistic Regression Using StatsmodelHere is a function to do the same sort of logistic regression using statsmodel.<jupyter_code>def logit(y,X):
    import statsmodels.api as sm

    X_ = pd.DataFrame()
   
    X_['X'] = X
    X_['const'] = 1

    logit_model=sm.Logit(y,X_)
    result=logit_model.fit()
    print(result.summary())

    # Odds ratio and conf interval
    params = result.params
    conf = result.conf_int()
    conf['OR'] = params
    conf.columns = ['2.5%', '97.5%', 'OR']
    print(np.exp(conf))<jupyter_output><empty_output><jupyter_text>Let's do the same regresion as we did in the previous section with sklear and compare the coefficients. They are the same.<jupyter_code>logit(lr['y'],lr['X'])<jupyter_output>Optimization terminated successfully.
         Current function value: 0.558891
         Iterations 6
                           Logit Regression Results                           
==============================================================================
Dep. Variable:                      y   No. Observations:                40349
Model:                          Logit   Df Residuals:                    40347
Method:                           MLE   Df Model:                            1
Date:                Fri, 17 Aug 2018   Pseudo R-squ.:                 0.01177
Time:                        01:21:00   Log-Likelihood:                -22551.
converged:                       True   LL-Null:                       -22819.
                                        LLR p-value:                8.516e-119
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
-----------------------------[...]<jupyter_text>Logistic Regression Exploration<jupyter_code>sklearn_lr(
    y = bw_env.lresult_gt_lmedian.astype(int), 
    X = bw_env.recent_precip.fillna(0)
)<jupyter_output><empty_output><jupyter_text>I think this means that a 1 unit ( 1 inch ) increase in the recent precipitation increases the odds of above median bacteria counts by a factor of 6.<jupyter_code>sklearn_lr(
    y = bw_env.lresult_gt_lmedian.astype(int), 
    X = bw_env.recent_precip.fillna(0)
)<jupyter_output><empty_output><jupyter_text>Rain SeasonalityMaybe there is a stonger correlation to precipitation in the first rains of the season? The guess is that late in the season, pollutants and nutrients have already been washed out by the early rains, so there should be more pollution earlier in the season.<jupyter_code>trr.groupby(trr.index.month).sum().recent_precip.plot()<jupyter_output><empty_output><jupyter_text>Here we calculate rainfall quartiles. A value of 1, the first quartile, includes all of the months of the year that have total rainfall less than the 25% percentile ( actually 30th, see comment ) while a value of 4 is months with rainfall greater than the 75% percentile.<jupyter_code># Compute the quartile of rainfall for each month

_1 = pd.qcut(trr.groupby(trr.index.month).sum().recent_precip, 
             q=[0, .3, .5, .75, 1.], # .3 not .25 to smooth out summer
            labels=[1,2,3,4])

dmap = {e['date']:e['recent_precip'] for _,e in _1.to_frame().reset_index().iterrows()}

#trr.apply( lambda row: dmap[row.index.month])
bw_env['month'] = bw_env.index.month
bw_env['precip_month_q'] = bw_env.month.apply(lambda v: dmap[v])

bw_env.groupby(bw_env.month).lresult_gt_lmedian.sum().plot()<jupyter_output><empty_output><jupyter_text>By grouping by the month quartile, we can see that days in the 75% percentile months get about 50x the rain of the 25% percentile months.<jupyter_code>bw_env.groupby(bw_env['precip_month_q']).recent_precip.mean()<jupyter_output><empty_output><jupyter_text>However, there are more days of above median results in the 1st quartile months than in any of the others.<jupyter_code>bw_env.groupby(bw_env['precip_month_q']).lresult_gt_lmedian.sum()
_1 = bw_env[(bw_env.measure_code==24)]

fig, ax = plt.subplots(figsize=(12,8))

ax = sns.scatterplot(ax=ax, x="recent_precip", y="lresult", hue="precip_month_q", palette="Set2", data=_1)<jupyter_output><empty_output><jupyter_text>The quantization in the chart above -- the vertical strips of points -- is very suspicious. Lets look at the records for the largest value.  It appears that there were many measreuemnts taken at many different stations, so they all share the same value for ``recent_precip``<jupyter_code>bw_env[bw_env.recent_precip == bw_env.recent_precip.max()][['stationcode', 'measure_code','collectiontime']].head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11115421920465568, 'fraction_numerical': 0.014096346589072098, 'mean_word_length': 4.807735636500188, 'pattern_counts': {'":': 0, '<': 99, '<?xml version=': 0, '>': 110, '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': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '91546', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 5023, 'n_tokens_neox': 4680, 'n_words': 1549, 'path': 'Beachwatch Logistic.ipynb', 'repo_name': 'san-diego-water-quality/example-analysis', 'size': '15465'}