In regression problems, we want to find a function that describes the relationship between observations and outputs. Our observations are recorded as rows in the matrix $X$, our response values as the column vector $Y$, and our estimate linear parameters as the column vector $\beta$. We think of our linear regression model as having a structural part and a random part: $Y= X\beta + \epsilon$ where $\epsilon$ ~ $N(0,1)$. The first two terms are structural and capture the linear relationship. The last is noise or variance in the data.
Linear regression doesn't require that the residuals $|r=Y - pred(X)|$ are normally distributed, but if they are, then the least-squares solution has many nice theoretical propoerties. Let's generate some fake data where the predictor and response are linearly related, and the residuals are normally distributed.
In [2]:
import numpy as np
from sklearn import cross_validation as cv
import matplotlib.pyplot as plt
%matplotlib inline
# use 100 samples
n_pts =100
beta= np.array([3.4])
X = np.expand_dims(np.linspace(0.001, 10, n_pts), axis=1)
def lin_resp_norm_err(beta, X):
return np.dot(X, beta[:, np.newaxis]) + np.expand_dims(np.random.randn(n_pts), axis=1)
Y_lin_norm_err = lin_resp_norm_err(beta, X)
Let's split out data into test and train sets so that we can examine how our linear model performs on held out data.
In [3]:
X_train, X_test, Y_train, Y_test = cv.train_test_split(X, Y_lin_norm_err, test_size=0.33, random_state=42)
We will now train a linear model with an ordinary least squares predictor. We will calculate the mean squre error on the hold out set as well as the $R^{2}$ value for the fit.
In [4]:
from sklearn import linear_model
lin_regr = linear_model.LinearRegression()
lin_regr.fit(X_train, Y_train)
residuals = lin_regr.predict(X_test) - Y_test
mean_error = np.mean(residuals ** 2)
regr_score = lin_regr.score(X_test, Y_test)
plt.scatter(X, Y_lin_norm_err)
plt.plot(X, lin_regr.predict(X), color='red')
plt.text(1, 30, '$r^{2}=%s$' % regr_score)
plt.xlabel('Predictor Variable X')
plt.ylabel('Response Variable Y')
plt.show()
The norms of the residuals follow a normal distribution.
In [5]:
# the histogram of the data with histtype='step'
plt.hist(residuals, bins=25)
l = plt.show()
We commonly think of linear regression containing a structural ($X\beta$) part and a random ($\epsilon$) part. In generalized linear models, we instead focus on the structural component, the link function, and the response distribution. We think of $g(\mu)=X\beta$ where $g$ is the link function, $X\beta$ is the structural component, and $\mu$ is the conditional mean of the distribution at that point in feature space. Note that we can pad the feature vectors of our observations with a '1' to capture the offset in the model, exactly like we do in normal linear regression.
If the link function, $g$, is the identity, the repsonse funciton is linear, and the response function is conditionally normally distributed, we are back again at normal linear regression. Here is a pretty great explanation someone wrote about glms overs on stack overflow.
Count data is fundamentally different from continuous data-- it is always discrete (integer values) and nonnegative. Let's grab some real data to analyze. Here, we are grabbing grouped individual data where there is an indicator for geographic "cell" number, residence (1=Suva, 2=Urban, or 3=Rural), education level of the mother (1=none, 2=primary, 3=upper primary, 4=secondary), The mean number of children ever born, variance in children ever born, number of owmen in the cell.
In [56]:
import urllib2
import pandas as pd
import StringIO
response = urllib2.urlopen('http://data.princeton.edu/wws509/datasets/ceb.dat')
children_ever_born = pd.read_csv(StringIO.StringIO(response.read()), sep='\s*')
Out[56]:
In [62]:
sm.GLM(children_ever_born['mean'], children_ever_born[['dur', 'res', 'educ', 'y']], family=sm.families.Poisson())
If we plot the predictor and response variables we can see what we would expect out of the poisson distribution. As the mean increases, so does the variance in the data.
In [28]:
X_train, X_test, Y_p_train, Y_p_test = cv.train_test_split(mean, Y_poisson, test_size=0.33, random_state=42)
plt.scatter(X_test, Y_p_test)
Out[28]:
In [31]:
glm_poisson = sm.GLM(Y_p_train, X_train, family=sm.families.Poisson())
res = glm_poisson.fit()
In [51]:
Out[51]:
In [50]:
Out[50]:
In [ ]: