Question

Cook’s distance is used to estimate the influence of a data point when performing least squares regression analysis. It is one of the standard plots for linear regression in R and provides another example of the applicationof leave-one-out resampling.

\[D_i = \frac{\sum_{j=1}^n (\hat Y_j - \hat Y_{j(i)})^2}{p\ \text{MSE}}\]

The calculation of Cook’s distance involves the fitting of \(n\) regression models, so we want to do this as efficiently as possible.

def cook_dist(X, y, model):
    """Vectorized version of Cook's distance."""
    n = len(X)
    fitted = model(y, X).fit()
    yhat = fitted.predict(X)
    p = len(fitted.params)
    mse = np.sum((yhat - y)**2.0)/n
    denom = p*mse
    idx = np.arange(n)
    return np.array([np.sum((yhat - model(y[idx!=i], X[idx!=i]).fit().predict(X))**2.0) for i in range(n)])/denom

import statsmodels.api as sm

# create data set with outliers
nobs = 100
X = np.random.random((nobs, 2))
X = sm.add_constant(X)
beta = [1, .1, .5]
e = np.random.random(nobs)
y = np.dot(X, beta) + e
y[[7, 29, 78]] *= 3

# use Cook's distance to identify outliers
model = sm.OLS
d = cook_dist(X, y, model)
plt.stem(d);