<jupyter_start><jupyter_text>Lambda School Data Science - Quantile RegressionRegressing towards the median - or any quantile - as a way to mitigate outliers and control risk.  LectureLet's look at data that has a bit of a skew to it:http://archive.ics.uci.edu/ml/datasets/Beijing+PM2.5+Data<jupyter_code>import pandas as pd
df = pd.read_csv('http://archive.ics.uci.edu/ml/machine-learning-databases/'
                 '00381/PRSA_data_2010.1.1-2014.12.31.csv')
df.head()
df.describe()
df['pm2.5'].plot.hist();
# How does linear regression handle it?
from sklearn.linear_model import LinearRegression

# Let's drop NAs and limit to numeric values
df = df._get_numeric_data().dropna()
X = df.drop('pm2.5', axis='columns')
y = df['pm2.5']

linear_reg = LinearRegression().fit(X, y)
linear_reg.score(X, y)
# Not bad - but what if we wanted to model the distribution more conservatively?
# Let's try quantile
import statsmodels.formula.api as smf

# Different jargon/API in StatsModel documentation
# "endogenous" response var is dependent (y), it is "inside"
# "exogenous" variables are independent (X), it is "outside"
# Bonus points - talk about "exogenous shocks" and you're a bona fide economist

# ~ style formulas look like what R uses
# y ~ x1 + x2 + ...
# Also, these formulas break with . in variable name, so lets change that
df = df.rename(index=str, columns={'pm2.5': 'pm25'})

# Now let's construct the formula string using all columns
quant_formula = 'pm25 ~ ' + ' + '.join(df.drop('pm25', axis='columns').columns)
print(quant_formula)

quant_mod = smf.quantreg(quant_formula, data=df)
quant_reg = quant_mod.fit(q=.5)
quant_reg.summary()  # "summary" is another very R-thing<jupyter_output>pm25 ~ No + year + month + day + hour + DEWP + TEMP + PRES + Iws + Is + Ir<jupyter_text>That fit to the median (q=0.5), also called "Least Absolute Deviation." The pseudo-R^2 isn't really directly comparable to the R^2 from linear regression, but it clearly isn't dramatically improved. Can we make it better?<jupyter_code>help(quant_mod.fit)
quantiles = (.05, .96, .1)

for quantile in quantiles:
    print(quant_mod.fit(q=quantile).summary())<jupyter_output>QuantReg Regression Results                          
==============================================================================
Dep. Variable:                   pm25   Pseudo R-squared:              0.04130
Model:                       QuantReg   Bandwidth:                       8.908
Method:                 Least Squares   Sparsity:                        120.7
Date:                Wed, 23 Jan 2019   No. Observations:                41757
Time:                        12:35:20   Df Residuals:                    41745
                                        Df Model:                           11
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept   3.072e-05    6.4e-06      4.803      0.000    1.82e-05    4.33e-05
No         -6.994e-05   9.59e-06     -7.292      0.0[...]<jupyter_text>"Strong multicollinearity", eh? In other words - maybe we shouldn't throw every variable in our formula. Let's hand-craft a smaller one, picking the features with the largest magnitude t-statistics for their coefficients. Let's also search for more quantile cutoffs to see what's most effective.<jupyter_code>quant_formula = 'pm25 ~ DEWP + TEMP + Ir + hour + Iws'
quant_mod = smf.quantreg(quant_formula, data=df)
for quantile in range(50, 100):
    quantile /= 100
    quant_reg = quant_mod.fit(q=quantile)
    print((quantile, quant_reg.prsquared))
# Okay, this data seems *extremely* skewed
# Let's trying logging
import numpy as np

df['pm25'] = np.log(1 + df['pm25'])
quant_mod = smf.quantreg(quant_formula, data=df)
quant_reg = quant_mod.fit(q=.25)
quant_reg.summary()  # "summary" is another very R-thing<jupyter_output><empty_output><jupyter_text>Overall - in this case, quantile regression is not *necessarily* superior to linear regression. But it does give us extra flexibility and another thing to tune - what the center of what we're actually fitting in the dependent variable.The basic case of `q=0.5` (the median) minimizes the absolute value of residuals, while OLS minimizes the squared value. By selecting `q=0.25`, we're targeting a lower quantile and are effectively saying that we only want to over-estimate at most 25% of the time - we're being *risk averse*.Depending on the data you're looking at, and the cost of making a false positive versus a false negative, this sort of flexibility can be extremely useful.Live - let's consider another dataset! Specifically, "SkillCraft" (data on competitive StarCraft players): http://archive.ics.uci.edu/ml/datasets/SkillCraft1+Master+Table+Dataset<jupyter_code># TODO Live!
# Hint - we may only care about the *top* quantiles here
# Another hint - there are missing values, but Pandas won't see them right away

skill = pd.read_csv('http://archive.ics.uci.edu/ml/machine-learning-databases/00272/SkillCraft1_Dataset.csv')
import numpy as np
skill.replace('?', np.nan, inplace=True)
skill.isnull().sum()
# replacing missing values for Age, HoursPerWeek, TotalHours
skill['Age'] = skill['Age'].fillna(skill['Age'].median())
skill['HoursPerWeek'] = skill['HoursPerWeek'].fillna(skill['HoursPerWeek'].median())
skill['TotalHours'] = skill['TotalHours'].fillna(skill['TotalHours'].median())
print ('Total nulls', skill.isnull().sum().sum())
# check dtypes
skill.dtypes
# correcting object types
skill['Age'] = pd.to_numeric(skill['Age'])
skill['HoursPerWeek'] = pd.to_numeric(skill['HoursPerWeek'])
skill['TotalHours'] = pd.to_numeric(skill['TotalHours'])
# predicting top 10% of players in APM
target = 'APM'
features = skill.columns.drop(target)

X = skill[features]
y = skill[target]
quant_formula_skill = 'APM ~ ' + ' + '.join(X.columns)
print(quant_formula_skill)

quant_mod = smf.quantreg(quant_formula_skill, data=skill)
quant_reg = quant_mod.fit(q=.9)
quant_reg.summary()  # "summary" is another very R-thing<jupyter_output>APM ~ GameID + LeagueIndex + Age + HoursPerWeek + TotalHours + SelectByHotkeys + AssignToHotkeys + UniqueHotkeys + MinimapAttacks + MinimapRightClicks + NumberOfPACs + GapBetweenPACs + ActionLatency + ActionsInPAC + TotalMapExplored + WorkersMade + UniqueUnitsMade + ComplexUnitsMade + ComplexAbilitiesUsed<jupyter_text>LeagueIndex is a strong indicator of APM.<jupyter_code>quant_reg.predict(X.iloc[0, :])<jupyter_output><empty_output><jupyter_text>Assignment - birth weight dataBirth weight is a situation where, while the data itself is actually fairly normal and symmetric, our main goal is actually *not* to model mean weight (via OLS), but rather to identify mothers at risk of having children below a certain "at-risk" threshold weight.Quantile regression gives us just the tool we need. For the data we are using, see: http://people.reed.edu/~jones/141/BirthWgt.html    bwt: baby's weight in ounces at birth    gestation: duration of pregnancy in days    parity: parity indicator (first born = 1, later birth = 0)    age: mother's age in years    height: mother's height in inches    weight: mother's weight in pounds (during pregnancy)    smoke: indicator for whether mother smokes (1=yes, 0=no)     Use this data and `statsmodels` to fit a quantile regression, predicting `bwt` (birth weight) as a function of the other covariates. First, identify an appropriate `q` (quantile) to target a cutoff of 90 ounces - babies above that birth weight are generally healthy/safe, babies below are at-risk.Then, fit and iterate your model. Be creative! You may want to engineer features. Hint - mother's age likely is not simply linear in its impact, and the other features may interact as well.At the end, create at least *2* tables and *1* visualization to summarize your best model. Then (in writing) answer the following questions:- What characteristics of a mother indicate the highest likelihood of an at-risk (low weight) baby?- What can expectant mothers be told to help mitigate this risk?Note that second question is not exactly a data science question - and that's okay! You're not expected to be a medical expert, but it is a good exercise to do a little bit of digging into a particular domain and offer informal but informed opinions.<jupyter_code># imports
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
bwt = pd.read_csv('http://people.reed.edu/~jones/141/Bwt.dat')
bwt.head()
bwt.describe()<jupyter_output><empty_output><jupyter_text>'q' really doesnt have much to do with our 'cutoff' value, it has much more to do with our desired confidence in predictive output. For this analysis, let's say I want to make sure only 5% of true babies' weights are below my prediction.It is much safer to operate with the assumption that a baby may be underweight (false positive) than to have too many false negatives, resulting in a baby not getting the proper treatment.<jupyter_code># predicting baby weights at 5% quantile
target = 'bwt'
features = bwt.columns.drop(target)

X = bwt[features]
y = bwt[target]
quant_formula_bwt = 'bwt ~ ' + ' + '.join(X.columns)
print(quant_formula_bwt)

quant_mod = smf.quantreg(quant_formula_bwt, data=bwt)
quant_reg = quant_mod.fit(q=.05)
quant_reg.summary()  # "summary" is another very R-thing<jupyter_output>bwt ~ gestation + parity + age + height + weight + smoke<jupyter_text>A couple issues from this first regression1. MutlicollinearityThis likely results from age, parity, height, and weight in some combination. Parity does not seem to have much predictive power, maybe it would be best to drop this. Let's see how it behaves in the next regression and go from there.Height and weight may be best to combine into a relative measure of the two, such as BMI.2. Age is not statistically significantThere is likely some inflection point for age. A 45+ year old woman having a child would intuitively seem more likely to have complications. Next time, we will try including this term.<jupyter_code># BMI term
bwt['bmi'] = bwt['weight'] / (bwt['height']**2)
# age squared term
bwt['age_squared'] = bwt['age'] ** 2
# re-run regression
# predicting baby weights at 5% quantile
target = 'bwt'
features = bwt.columns.drop([target, 'weight', 'height'])
X = bwt[features]
y = bwt[target]

quant_formula_bwt = 'bwt ~ ' + ' + '.join(X.columns)
print(quant_formula_bwt)

quant_mod = smf.quantreg(quant_formula_bwt, data=bwt)
quant_reg = quant_mod.fit(q=.05)
quant_reg.summary()  # "summary" is another very R-thing<jupyter_output>bwt ~ gestation + parity + age + smoke + bmi + age_squared<jupyter_text>Parity still isn't statistically significant. Let's drop that.Age, as well as its squared term, are not statistically significant. This is puzzling.Let's take a step back and and plot some of our features against baby weight.<jupyter_code># pairplot with parity as hue
sns.pairplot(data=bwt, x_vars=bwt.columns.drop(['bwt', 'parity', 'smoke']), y_vars='bwt', hue='parity');
# pairplot with smoker as hue
sns.pairplot(data=bwt, x_vars=bwt.columns.drop(['bwt', 'parity', 'smoke']), y_vars='bwt', hue='smoke');<jupyter_output><empty_output><jupyter_text>From the above, we can gain some intuition for relationships we haven't formalized in the model yet, including:- young smokers seem more prone to have babies that weight less, especially the age_squared term- smokers with early births (<260 weeks) seem prone to have babies that weight less- age doesn't really seem to be strongly related to low-weight babies from these charts, confirming some of our intuition<jupyter_code># smokers with early births
bwt['smoke_earlybirth'] = (bwt['gestation'] < 260) * bwt['smoke']
# early births in general
bwt['earlybirth'] = bwt['gestation'] < 260
# age squared and smoke interaction term
bwt['age_squared_smoke'] = bwt['age_squared'] * bwt['smoke']
bwt.head()<jupyter_output><empty_output><jupyter_text>After adding some more variables, let's see what a full regression against the features looks like. We'll drop unnecessary variables from there.<jupyter_code># re-run regression
# predicting baby weights at 5% quantile
target = 'bwt'
features = bwt.columns.drop(target)
X = bwt[features]
y = bwt[target]

quant_formula_bwt = 'bwt ~ ' + ' + '.join(X.columns)
print(quant_formula_bwt)

quant_mod = smf.quantreg(quant_formula_bwt, data=bwt)
quant_reg = quant_mod.fit(q=.05)
quant_reg.summary()  # "summary" is another very R-thing<jupyter_output>bwt ~ gestation + parity + age + height + weight + smoke + bmi + age_squared + smoke_earlybirth + age_squared_smoke + earlybirth<jupyter_text>An early birth seems to be by far the most powerful predictor. Let's drop features that aren't significant at alpha=.1 and re-run.<jupyter_code># re-run regression
# predicting baby weights at 5% quantile
target = 'bwt'
features = bwt.columns.drop([target, 'parity', 'age', 'height', 'weight', 'bmi', 
                             'age_squared', 'age_squared_smoke', 'smoke_earlybirth'])
X = bwt[features]
y = bwt[target]

quant_formula_bwt = 'bwt ~ ' + ' + '.join(X.columns)
print(quant_formula_bwt)

quant_mod = smf.quantreg(quant_formula_bwt, data=bwt)
quant_reg = quant_mod.fit(q=.05)
quant_reg.summary()  # "summary" is another very R-thing<jupyter_output>bwt ~ gestation + smoke + earlybirth<jupyter_text>Our pseudo R-squared barely moves in this case. I know I've been cheating a bit by working with all of my data, but let's try this model for real by splitting it into a train and test set.To validate, we will see what percent of the true test observations fall above and below the predicted value. If our model is working well, around 5% of results will fall below our predicted value.<jupyter_code># randomly split data into train and test
np.random.seed(1)
bwt['split'] = np.random.choice(['train', 'test'], size=bwt.shape[0], p=[0.67, 0.33])

bwt_train, bwt_test = bwt[bwt['split'] == 'train'], bwt[bwt['split'] == 'test']

X_train, y_train = bwt_train[features], bwt_train[target]
X_test, y_test = bwt_test[features], bwt_test[target]

quant_formula_bwt = 'bwt ~ ' + ' + '.join(X_train.columns)
print(quant_formula_bwt)

quant_mod = smf.quantreg(quant_formula_bwt, data=bwt_train)
quant_reg_final = quant_mod.fit(q=.05)
quant_reg_final.summary()  # "summary" is another very R-thing
# validating results output
results = pd.DataFrame({'y_predicted': quant_reg_final.predict(X_test), 'y_true': y_test})
true_lt_predicted = (results['y_true'] < results['y_predicted']).astype(int)

print ('Actual percent of predictions lying below predicted y value:',
       round((true_lt_predicted.sum() * 1.0)/true_lt_predicted.count(), 3))<jupyter_output>Actual percent of predictions lying below predicted y value: 0.921<jupyter_text>Summary Risk Factors and MitigationThe coefficient estimations of our quartile regression model (q=0.05) is printed below.The strongest predictors of the 5th quartile in birth weight are1. An Early Birth : defined as a birth occurring earlier than 260 days2. Gestation Period : the absolute amount of days a pregnancy has lasted3. Smoking : A Mother who smokes decreases the expected 5th quartile birth rateNote: late births likely come with their own complications, so 'gestation' should not be considered an absolutely beneficial feature of pregnancy.<jupyter_code>quant_reg_final.summary()<jupyter_output><empty_output><jupyter_text>To mitigate risks, mothers should avoid smoking and take steps (to the extent possible) to prevent early births from occurring. The Mayo Clinic lists the following risk factors for premature births:- Having a previous premature birth- Pregnancy with twins, triplets or other multiples- An interval of less than six months between pregnancies- Conceiving through in vitro fertilization- Problems with the uterus, cervix or placenta- Smoking cigarettes or using illicit drugs- Some infections, particularly of the amniotic fluid and lower genital tract- Some chronic conditions, such as high blood pressure and diabetes- Being underweight or overweight before pregnancy- Stressful life events, such as the death of a loved one or domestic violence- Multiple miscarriages or abortions- Physical injury or trauma PerformanceAs shown above, we an compare our predictions on the test data to the actual amount of baby weight points lying below our estimates. In our test data, only 5.6% of points lie below our estimated weight.<jupyter_code># validating results output
results = pd.DataFrame({'y_predicted': quant_reg_final.predict(X_test), 'y_true': y_test})
true_lt_predicted = (results['y_true'] < results['y_predicted']).astype(int)

print ('Actual percent of predictions lying below predicted y value:',
       round((true_lt_predicted.sum() * 1.0)/true_lt_predicted.count(), 3))<jupyter_output>Actual percent of predictions lying below predicted y value: 0.056<jupyter_text>Let's visualize our only final continuous feature, gestation period, and compare our predicted values to the actual results.<jupyter_code># create sample data for gestation, assuming smoking and no smoking
gestation = np.arange(175, 350, 0.5)
early_birth = np.where(gestation < 260, 1, 0)
X_plot_smoke = pd.DataFrame({'gestation' : gestation, 'smoke' : 1, 'earlybirth' : early_birth})
X_plot_nosmoke = pd.DataFrame({'gestation' : gestation, 'smoke' : 0, 'earlybirth' : early_birth})

y_plot_smoke, y_plot_nosmoke = quant_reg_final.predict(X_plot_smoke), quant_reg_final.predict(X_plot_nosmoke)
# plotting line of best fit for smokers
fig_smoke, ax_smoke = plt.subplots(figsize=(13,9))

ax_smoke.plot(X_plot_smoke['gestation'], y_plot_smoke, label='Best Fit Line, Smokers', c='blue')
ax_smoke.scatter(X_test['gestation'][X_test['smoke']==True], y_test[X_test['smoke']==True], 
                 alpha=0.2, 
                 c='blue',
                 label='Test Data: Smokers')

# plotting line of best fit and scatter plot for non-smokers
ax_smoke.plot(X_plot_nosmoke['gestation'], y_plot_nosmoke, label='Best Fit Line, Non-Smokers', c='red')
ax_smoke.scatter(X_test['gestation'][X_test['smoke']==False], y_test[X_test['smoke']==False], 
                 alpha=0.2, 
                 c='red',
                 label='Test Data: Non-Smokers')

# cutoff for early birth binary variable
ax_smoke.axvline(x=260, c='black', linestyle='--', label='Early Birth Cutoff (< 260 Days)')

# labeling
ax_smoke.set_ylabel('Baby Weight (oz)')
ax_smoke.set_xlabel('Pregnancy Lenth (days)')
ax_smoke.legend(loc='upper left')
ax_smoke.set_title('Predicting Baby Weight With Quantile Regression (q=0.05)')

# annotations
ax_smoke.text(310, 150, '95% of True Points Should \nLie Above Best Fit Line', fontsize=12)
ax_smoke.text(300, 60, '5% of True Points Should \nLie Below Best Fit Line', fontsize=12)

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1128887480859602, 'fraction_numerical': 0.016209937166693066, 'mean_word_length': 4.605208641609944, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 57, 'https://': 0, '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': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '589407', 'license': 'MIT', 'n_tokens_mistral': 6149, 'n_tokens_neox': 5664, 'n_words': 2024, 'path': 'module3-quantile-regression/LS_DS_233_Quantile_Regression.ipynb', 'repo_name': 'zangell44/DS-Unit-2-Sprint-3-Advanced-Regression', 'size': '18939'}