## Car Emissions Study¶

This subset (9x series) of notebooks was created for publication in a series of blogposts for Applied AI. They are part of a wider project primarily serving as a coherent collection of my preferred techniques for data preparation & analysis and Bayesian inference in Python.

A set of notebooks designed to investigate car emissions data from the PoV of the Volkswagen Emissions Scandal which seems to have meaningfully damaged their sales. The motivation is to investigate the data and see if we can see any unsual behaviour in Volkswagen.

Using data from the UK VCA (Vehicle Type Approval) Car Fuel and Emissions Information for August 2015. Dataset available here and also included in the repo since it's small.

# 92_FeatureSelectionAndModelEvaluation¶

#### Demonstrate linear regression, Lasso regularization for feature selection, and posterior predictive checks¶

Notes:

• Python 3.5 project using the latest available PyMC3
• Developed using ContinuumIO Anaconda distribution on a Macbook Pro 3GHz i7, 16GB RAM, OSX 10.10.5.
• If execution becomes unstable or Theano throws weird errors, try clearing the cache $> theano-cache clear and rerunning the notebook. Package Requirements (shown as a conda-env YAML): $> less conda_env_pymc3_examples.yml

name: pymc3_examples
channels:
- defaults
dependencies:
- python=3.5
- jupyter
- ipywidgets
- numpy
- scipy
- matplotlib=1.4.3
- pandas
- scikit-learn
- seaborn
- patsy
- pip

$> conda env create --file conda_env_pymc3_examples.yml$> source activate pymc3_examples

#### Declare full modelspec¶

In [86]:
fml_all = '{} ~ '.format(ft_endog) + ' + '.join(fts_num + fts_cat)
fml_all

Out[86]:
'emissions_nox_mgkm ~ metric_combined + metric_extra_urban + metric_urban_cold + engine_capacity + emissions_co_mgkm + mfr_owner_is_vw + trans + fuel_type + is_tdi'
##### Create design matrices for statsmodels¶
In [87]:
(mx_en, mx_ex) = pt.dmatrices(fml_all, dfs, return_type='dataframe', NA_action='raise')
# custom_describe(mx_ex, 2, )


## Frequentist OLS Regression¶

For later comparison, first let's use statsmodels to run a Freqentist OLS

In [88]:
smfit = sm.OLS(mx_en, mx_ex).fit()
smfit.summary()

Out[88]:
Dep. Variable: R-squared: emissions_nox_mgkm 0.467 OLS 0.464 Least Squares 209.0 Mon, 29 Feb 2016 0.00 10:40:26 -10585. 2641 2.119e+04 2629 2.126e+04 11 nonrobust
coef std err t P>|t| [95.0% Conf. Int.] 46.2939 0.668 69.332 0.000 44.985 47.603 4.2480 0.883 4.809 0.000 2.516 5.980 -1.1535 0.635 -1.816 0.070 -2.399 0.092 -2.5677 1.001 -2.564 0.010 -4.532 -0.604 -15.8664 2.114 -7.505 0.000 -20.012 -11.721 -19.3263 1.009 -19.151 0.000 -21.305 -17.348 5.7168 1.198 4.772 0.000 3.368 8.066 -52.0603 21.036 -2.475 0.013 -93.310 -10.811 19.2034 9.017 2.130 0.033 1.523 36.884 24.2640 12.603 1.925 0.054 -0.449 48.977 8.6972 1.252 6.947 0.000 6.242 11.152 2.3939 0.621 3.853 0.000 1.176 3.612
 Omnibus: Durbin-Watson: 10.327 1.06 0.006 9.358 0.1 0.00929 2.787 129

Observe

• That was easy! statsmodels is great for basic stuff
• The R-squared of 0.467 isn't too bad considering the possible range (-inf,1) (NOTE I explain a little more about r-squared theory towards the end of this notebook)
• The condition number of 129 is far above 20, the recommended theshold at which we should consider the effects of multicollinearity

I won't get into the actual interpretation of the coefficient values yet.

NOTE

• Just in case you missed it, I used patsy above to create 'design matrices' for the data, prior to modelling with statsmodels. This converted the main dataframe to the same 'modelspec' as I will use throughout the Frequentist and Bayesian modelling.
• The categorical features have been binarised (a.k.a one-hot encoded) and the Intercept coefficient is overloaded with the first value from each categorical feature to allow for proper identifiability: i.e:
• if a datapoint had raw feature value trans == manual, that is now indicated by a boolean True in the new column trans[T.manual], and a boolean False in the new column trans[T.semiauto]
• if a datapoint had raw feature value trans == auto, that is now indicated by a boolean False in the new columns trans[T.manual] and trans[T.semiauto]: the Intercept column always has value 1 aka True, meaning that trans == auto is represented by the Intercept.
• The overloading means that the Intercept represents cars with categorical values trans == auto, mfr_owner_is_vw == False, fuel_type == diesel, and is_tdi == False.
• Due to our standardization (mean-centering and dividing by the std.dev.) the Intercept also represents the mean of the numeric features.

## Bayesian OLS Regression¶

Okay, time to use pymc3! Lets create the same OLS model, again using:

• the glm submodule for convenience
• the NUTS sampler for 'better' convergence

### Define and run model¶

In [89]:
with pm.Model() as mdl_ols:

## Use GLM submodule for simplified model specification
## Betas are Uniform (for OLS)
## Likelihood is Normal (with HalfCauchy for error prior)

pm.glm.glm('{} ~ '.format(ft_endog) + ' + '.join(fts_num + fts_cat)
,dfs
,intercept_prior=pm.Uniform.dist(lower=-1e6, upper=1e6)
,regressor_prior=pm.Uniform.dist(lower=-1e6, upper=1e6)
,family=pm.glm.families.Normal())

## find MAP using Powell
start_MAP = pm.find_MAP(fmin=optimize.fmin_powell)

## Sample using NUTS
trc_ols = pm.sample(2000, start=start_MAP, step=pm.NUTS())

# convenience: declare Random Variables (RVs) _not_ created by the PyMC3 backend
rvs = [rv.name for rv in mdl_ols.unobserved_RVs]
_ = [rvs.remove(rv.name) for rv in mdl_ols.free_RVs]

Applied interval-transform to Intercept and added transformed Intercept_interval to model.
Applied interval-transform to mfr_owner_is_vw[T.True] and added transformed mfr_owner_is_vw[T.True]_interval to model.
Applied interval-transform to trans[T.manual] and added transformed trans[T.manual]_interval to model.
Applied interval-transform to trans[T.semiauto] and added transformed trans[T.semiauto]_interval to model.
Applied interval-transform to fuel_type[T.hybrid] and added transformed fuel_type[T.hybrid]_interval to model.
Applied interval-transform to fuel_type[T.petrol] and added transformed fuel_type[T.petrol]_interval to model.
Applied interval-transform to is_tdi[T.True] and added transformed is_tdi[T.True]_interval to model.
Applied interval-transform to metric_combined and added transformed metric_combined_interval to model.
Applied interval-transform to metric_extra_urban and added transformed metric_extra_urban_interval to model.
Applied interval-transform to metric_urban_cold and added transformed metric_urban_cold_interval to model.
Applied interval-transform to engine_capacity and added transformed engine_capacity_interval to model.
Applied interval-transform to emissions_co_mgkm and added transformed emissions_co_mgkm_interval to model.
Applied log-transform to sd and added transformed sd_log to model.
[-----------------100%-----------------] 2000 of 2000 complete in 41.0 sec

Observe

• PyMC's default behaviour when creating the theano-based model is to be quite verbose, and you can see several printout lines informing us of various transforms added to the model
• The NUTS sampler ran 2000 iterations, with a single chain, taking under a minute.

### View feature coefficients¶

In [90]:
pm.df_summary(trc_ols[-1000:], varnames=rvs)

Out[90]:
mean sd mc_error hpd_2.5 hpd_97.5
Intercept 46.244592 0.691005 0.030061 44.765587 47.531293
mfr_owner_is_vw[T.True] 4.259546 0.943761 0.027260 2.308817 6.016226
trans[T.manual] -1.111725 0.631056 0.020960 -2.291213 0.197674
trans[T.semiauto] -2.528842 1.010407 0.027052 -4.494458 -0.545857
fuel_type[T.hybrid] -15.877565 2.238481 0.066838 -20.660791 -11.871254
fuel_type[T.petrol] -19.275976 1.015030 0.045213 -21.193061 -17.204596
is_tdi[T.True] 5.709995 1.282720 0.039323 3.230674 8.177947
metric_combined -53.521139 22.330069 1.836654 -95.289728 -13.675488
metric_extra_urban 19.805878 9.706478 0.791371 1.887635 36.715806
metric_urban_cold 25.112368 13.169643 1.071202 0.929582 49.626974
engine_capacity 8.723662 1.296313 0.048746 6.311740 11.279146
emissions_co_mgkm 2.355926 0.618600 0.018136 1.128695 3.493506
sd 13.345757 0.183700 0.006675 12.985165 13.694030

Observe:

• The above table summarises the final 1000 steps of the traces, giving us the basic statistics of the distributions of the parameter estimates.
• You can see the mean values are very similar to the statsmodels OLS model, which is good to see, and for reference they shown in the following cell:
In [91]:
## recap on the values from statsmodels OLS
pd.DataFrame.from_csv(StringIO(smfit.summary().tables[1].as_csv())).iloc[:,:1]

Out[91]:
coef
Intercept 46.2939
mfr_owner_is_vw[T.True] 4.2480
trans[T.manual] -1.1535
trans[T.semiauto] -2.5677
fuel_type[T.hybrid] -15.8664
fuel_type[T.petrol] -19.3263
is_tdi[T.True] 5.7168
metric_combined -52.0603
metric_extra_urban 19.2034
metric_urban_cold 24.2640
engine_capacity 8.6972
emissions_co_mgkm 2.3939

### View traceplots¶

Rather than just get point-estimate statistics from the traces, let's take a look at the traceplots using PyMC3's built-in pm.traceplot() function.

As I mentioned in the first blogpost:

• each feature coefficient is shown on a single row
• the right-hand-side plot is a simple timeseries of each value on the trace over the 1000 samples
• the left-hand-side plot is a density plot of the traces (mean shown in red) this is a marginal posterior distribution on each coefficient
In [92]:
plot_traces(trc_ols, retain=1000, varnames=rvs)