In [1]:

```
from __future__ import division
import pandas as pd
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import mean_squared_error
from sklearn.feature_selection import SelectKBest
from sklearn.feature_selection import f_regression
from sklearn.feature_selection import chi2
from sklearn.cluster import MeanShift
from scipy.interpolate import LSQUnivariateSpline
%matplotlib inline
```

In [2]:

```
# Data comes from ISLR::Wage - written using write.csv(Wage, "wage.csv", row.names=FALSE)
wage_df = pd.read_csv("../data/Wage.csv")
wage_df.head()
```

Out[2]:

In [3]:

```
# factorize non-numeric variables
wage_df["sex"] = pd.factorize(wage_df["sex"])[0]
wage_df["maritl"] = pd.factorize(wage_df["maritl"])[0]
wage_df["race"] = pd.factorize(wage_df["race"])[0]
wage_df["education"] = pd.factorize(wage_df["education"])[0]
wage_df["region"] = pd.factorize(wage_df["region"])[0]
wage_df["health"] = pd.factorize(wage_df["health"])[0]
wage_df["health_ins"] = pd.factorize(wage_df["health_ins"])[0]
wage_df.head()
```

Out[3]:

We attempt to build a model to predict wage given a degree 4 polynomial over age. In R the polynomial features are created using the poly() function. As pointed out in this post from Dave Moore the method creates orthogonal polynomial features, which are different from simple powers of the original feature. Code for the poly() function shown below is a straight copy from the ortho_poly_fit() function code from this page.

In [4]:

```
def poly(x, degree = 1):
n = degree + 1
x = np.asarray(x).flatten()
if(degree >= len(np.unique(x))):
print("'degree' must be less than number of unique points")
return
xbar = np.mean(x)
x = x - xbar
X = np.fliplr(np.vander(x, n))
q, r = np.linalg.qr(X)
z = np.diag(np.diag(r))
raw = np.dot(q, z)
norm2 = np.sum(raw ** 2, axis=0)
alpha = (np.sum((raw ** 2) * np.reshape(x,(-1,1)), axis=0) / norm2 + xbar)[:degree]
Z = raw / np.sqrt(norm2)
return Z, norm2, alpha
X = poly(wage_df["age"].values, 4)[0]
X[0:5, 1:]
```

Out[4]:

In [5]:

```
y = wage_df["wage"].values
reg = LinearRegression()
reg.fit(X, y)
print "Intercepts:", reg.intercept_
print "Coefficients:", reg.coef_
ax = wage_df.plot(x="age", y="wage", style="o")
ax.set_ylabel("wage")
# The poly() method cannot return features for a single X value, so we have
# to plot the raw predictions.
age = wage_df["age"].values
ypred = reg.predict(X)
polyline = np.poly1d(np.polyfit(age, ypred, 4))
xs = range(int(np.min(age)), int(np.max(age)))
ys = polyline(xs)
ax.plot(xs, ys, 'r', linewidth=2.5)
```

Out[5]:

We can repeat this experiment with manually calculated powers of the age feature as follows. Even though the parameters are different, the resulting linear model appears to be the same.

In [6]:

```
x1 = wage_df["age"].values
x2 = np.power(x1, 2)
x3 = np.power(x1, 3)
x4 = np.power(x1, 4)
X = np.vstack((x1, x2, x3, x4)).T
X[0:5, :]
```

Out[6]:

In [7]:

```
y = wage_df["wage"].values
reg = LinearRegression()
reg.fit(X, y)
print "Intercepts:", reg.intercept_
print "Coefficients:", reg.coef_
ax = wage_df.plot(x="age", y="wage", style="o")
ax.set_ylabel("wage")
xs = range(int(np.min(x1)), int(np.max(x1)))
ys = [reg.predict([x, x**2, x**3, x**4]) for x in xs]
ax.plot(xs, ys, 'r', linewidth=2.5)
```

Out[7]:

R's summary() function applied to a linear model returns p-values for each coefficient, which enables one to guage the significance of each feature. Scikit-Learn does not provide anything similar, but it does provide methods for feature selection which provides p-values and scores. Based on p-values and scores below, it appears that the 4th degree feature is the least important of the four. However, R reports a p-value of 0.05 for the 4th power.

In [8]:

```
selector = SelectKBest(score_func=f_regression, k=10)
selector.fit(X, y)
(selector.scores_, selector.pvalues_)
```

Out[8]:

The p-values observed from a linear model's summary() output only works as an indicator of feature importance if a single feature is being used. For multiple features, R advises the use of ANOVA. It turns out that the SelectKBest feature selector is based on ANOVA. The R example manually built models based on different feature sets then applied ANOVA to the models - there is no direct analog to that here.

Fitting a Logistic Regression model to a binary response variable. The response variable is "high-earners" where the salary > 250K.

In [9]:

```
X = np.vstack((x1, x2, x3)).T
y = wage_df["wage"].map(lambda w: 1 if w > 250 else 0).values
reg = LogisticRegression()
reg.fit(X, y)
print "Intercepts:", reg.intercept_
print "Coefficients:", reg.coef_
ypred = reg.predict(X)
print "MSE:", mean_squared_error(y, ypred)
# Plot the predicted probability of being a high-earner over age range
ages = wage_df["age"].values
xs = range(np.min(ages), np.max(ages))
ys = [reg.predict_proba(np.array([x, x*x, x*x*x]))[0][0] for x in xs]
plt.plot(xs, ys)
plt.xlabel("age")
plt.ylabel("p(wage > 250k)")
```

Out[9]:

Splines is a curve approximated by a series of straight lines. The breaks in the straight lines are called knots. The idea is to fit multiple linear models to these ranges, one each for each straight line segment. Scipy offers the UnivariateSpline but it requires manual knot selection. The code snippet below uses Mean-Shift clustering to automatically derive the knots based on this thread on StackOverflow.

In [10]:

```
ages = wage_df["age"].values
wages = wage_df["wage"].values
X = np.vstack((ages, wages)).T
# cluster points to find the knot location
msc = MeanShift()
msc.fit(X)
knots = msc.cluster_centers_[:, 0]
# fit a spline over the points
spl = LSQUnivariateSpline(ages, wages, knots)
xs = range(np.min(ages), np.max(ages))
ys = spl(xs)
# plot the points and the spline
ax = wage_df.plot(x="age", y="wage", style="o")
ax.set_ylabel("wage")
ax.plot(xs, ys, 'r', linewidth=2.5)
```

Out[10]:

Generalized Additive Models are not yet available under Scikit-Learn. It was one of the projects proposed for GSoC 2014 but as far as I know it did not get approved. StatsModel has an untested implementation. Since it is not available in Python, I decided to skip it for now.