<jupyter_start><jupyter_text>Interpreting linear models==========================Linear models are not that easy to interpret when variables arecorrelated.See also the [statistics chapter](http://www.scipy-lectures.org/packages/statistics/index.html) of the[scipy lecture notes](http://www.scipy-lectures.org). Data on wages--------------We first download and load some historical data on wages<jupyter_code>import os
import pandas

# Python 2 vs Python 3:
try:
    from urllib.request import urlretrieve
except ImportError:
    from urllib import urlretrieve

if not os.path.exists('wages.txt'):
    # Download the file if it is not present
    urlretrieve('http://lib.stat.cmu.edu/datasets/CPS_85_Wages',
                'wages.txt')

# Give names to the columns
names = [
    'EDUCATION: Number of years of education',
    'SOUTH: 1=Person lives in South, 0=Person lives elsewhere',
    'SEX: 1=Female, 0=Male',
    'EXPERIENCE: Number of years of work experience',
    'UNION: 1=Union member, 0=Not union member',
    'WAGE: Wage (dollars per hour)',
    'AGE: years',
    'RACE: 1=Other, 2=Hispanic, 3=White',
    'OCCUPATION: 1=Management, 2=Sales, 3=Clerical, 4=Service, 5=Professional, 6=Other',
    'SECTOR: 0=Other, 1=Manufacturing, 2=Construction',
    'MARR: 0=Unmarried,  1=Married',
]

short_names = [n.split(':')[0] for n in names]
data = pandas.read_csv('wages.txt', skiprows=27, skipfooter=6, sep=None,
                       header=None, engine='python')
data.columns = short_names

# Log-transform the wages, as they typically increase with
# multiplicative factors
import numpy as np
data['WAGE'] = np.log10(data['WAGE'])<jupyter_output><empty_output><jupyter_text>The challenge of correlated features--------------------------------------------Plot scatter matrices highlighting the links between differentvariables measured<jupyter_code>import seaborn
# The simplest way to plot a pairplot
seaborn.pairplot(data, vars=['WAGE', 'AGE', 'EDUCATION', 'EXPERIENCE', 'SEX'])<jupyter_output><empty_output><jupyter_text>A fancier pair plot<jupyter_code>from matplotlib import pyplot as plt

g = seaborn.PairGrid(data,
                     vars=['WAGE', 'AGE', 'EDUCATION', 'EXPERIENCE', 'SEX'],
                     diag_sharey=False)
g.map_upper(seaborn.kdeplot)
g.map_lower(plt.scatter, s=2)
g.map_diag(seaborn.kdeplot, lw=3)<jupyter_output><empty_output><jupyter_text>Note that age and experience are highly correlatedA link between a single feature and the target is a *marginal* link.Univariate feature selection selects on marginal links.Linear model compute *conditional* links: removing the effects of otherfeatures on each feature. This is hard when features are correlated. Coefficients of a linear model--------------------------------------------<jupyter_code>from sklearn import linear_model
features = [c for c in data.columns if c != 'WAGE']
X = data[features]
y = data['WAGE']
ridge = linear_model.RidgeCV()
ridge.fit(X, y)
# Visualize the coefs
coefs = ridge.coef_
from matplotlib import pyplot as plt
plt.figure(figsize=(6, 4))
plt.barh(np.arange(coefs.size), coefs)
plt.yticks(np.arange(coefs.size), features)
plt.title("Coefficients")
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>**Scaling coefficients**: coefs cannot easily be compared if X is notstandardized: they should be normalized to the variance of X: thegreater the variance of a feature, the large the impact of thecorresponding coefficent on the output.If the different features have differing, possibly arbitrary, scales,then scaling the coefficients by the feature scale often helpsinterpretation.<jupyter_code>X_std = X.std()
plt.figure(figsize=(6, 4))
plt.barh(np.arange(coefs.size), coefs * X_std)
plt.yticks(np.arange(coefs.size), features)
plt.title("Scaled coefficients")
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Now the age and experience can be better compared: and experience doesappear as more important than age. When features are not too correlated and there is plenty of samples,this is the well-known regime of standard statistics in linear models.Machine learning is not needed, and statsmodels is a great tool (see the[statistics chapter in scipy-lectures](http://www.scipy-lectures.org/packages/statistics/index.html)) The effect of regularization--------------------------------------------Sparse models use l1 regularization to puts some variables tozero. This can often help interpretation<jupyter_code>lasso = linear_model.LassoCV(max_iter=100000, cv=5)
lasso.fit(X, y)

coefs = lasso.coef_
plt.barh(np.arange(coefs.size), coefs * X_std)
plt.yticks(np.arange(coefs.size), features)
plt.title("Sparse coefficients")
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>When two variables are very correlated (such as age and experience), itwill put arbitrarily one or the other to zero depending on their SNR.Here we can see that age probably overshadowed experience. Stability to gauge significance--------------------------------Stability of coefficients when perturbing the data helps giving aninformal evaluation of the significance of the coefficients. Note thatthis is not significance testing in the sense of p-values, as a modelthat returns coefficients always at one indepently of the data willappear as very stable though it clearly does not control for falsedetections.We can do this in a cross-validation loop, using the argument"return_estimator" of `sklearn.model_selection.cross_validate`which has been added in version 0.20 of scikit-learn:<jupyter_code>from sklearn.model_selection import cross_validate<jupyter_output><empty_output><jupyter_text>With the lasso estimator<jupyter_code>cv_lasso = cross_validate(lasso, X, y, return_estimator=True, cv=10)

coefs_ = [estimator.coef_ for estimator in cv_lasso['estimator']]<jupyter_output><empty_output><jupyter_text>Plot the results with seaborn:<jupyter_code>coefs_ = pandas.DataFrame(coefs_, columns=features) * X_std
plt.figure(figsize=(6, 4))
seaborn.boxplot(data=coefs_, orient='h')
seaborn.stripplot(data=coefs_, orient='h', color='k')
plt.axvline(x=0, color='.5')  # Add a vertical line at 0
plt.title('Sparse coefficients')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>With the ridge estimator<jupyter_code>cv_ridge = cross_validate(ridge, X, y, return_estimator=True, cv=10)
coefs_ = [estimator.coef_ for estimator in cv_ridge['estimator']]

coefs_ = pandas.DataFrame(coefs_, columns=features) * X_std
plt.figure(figsize=(6, 4))
seaborn.boxplot(data=coefs_, orient='h')
seaborn.stripplot(data=coefs_, orient='h', color='k')
plt.axvline(x=0, color='.5') # Add a vertical line at 0
plt.title('Non-sparse coefficients (ridge)')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>The story is different. Note also that lasso coefficients are much moreinstable than the ridge. It is often the case that sparse models areunstable.  Which is the truth?Note the difference between the lasso and the ridge estimator: we donot have enough data to perfectly estimate conditional relationships,hence the prior (ie the regularization) makes a difference, and its ishard to tell from the data which is the "truth".One reasonnable model-selection criterion is to believe most the modelthat predicts best. For this, we can inspect the prediction scoresobtained via the cross-validation<jupyter_code>scores = pandas.DataFrame({'lasso': cv_lasso['test_score'],
                           'ridge': cv_ridge['test_score']})
plt.figure(figsize=(6, 2))
seaborn.boxplot(data=scores, orient='h')
seaborn.stripplot(data=scores, orient='h', color='k')
plt.title("Model comparison")<jupyter_output><empty_output><jupyter_text>Note also that the limitations of cross-validation explained previouslystill apply. Ideally, we should use a ShuffleSplit cross-validationobject to sample many times and have a better estimate of theposterior, both for the coefficients and the test scores. Conclusion on factors of wages?As always, concluding is hard. That said, it seems that we shouldprefer the scaled ridge coefficients.Hence, the most important factors of wage are education andexperience, followed by sex: at the same education and experiencefemales earn less than males. Note that this last statement is astatement about the link between wage and sex, conditional on educationand experience. Interpreting random forests=======================Interpreting random forests can be done through feature importance. This index explain how much a feature is used to split a tree. This is also known as the Gini importance.<jupyter_code>from sklearn import ensemble

forest = ensemble.RandomForestRegressor(n_estimators=100)
forest.fit(X, y)

# Visualize the feature importance
importance = forest.feature_importances_
from matplotlib import pyplot as plt
plt.barh(np.arange(importance.size), importance)
plt.yticks(np.arange(importance.size), features)
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Feature importances are a proxy for the mutual information between thefeature and the target, conditionally on the other features. Theconditioning is difficult and not well controlled.Higher-cardinality categorical variables will have largerfeature importances<jupyter_code>from sklearn.metrics import per<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11290670114199526, 'fraction_numerical': 0.007110536522301228, 'mean_word_length': 5.297829036635007, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, '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': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '139430', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 2889, 'n_tokens_neox': 2706, 'n_words': 912, 'path': '07_model_interpretation/01_interpreting_linear_models.ipynb', 'repo_name': 'lmercier36/datacamp-master', 'size': '9282'}