Python, Econometrics: Diagnostic measures, studentized deleted residuals, dffits and Cook distances.

# -*- coding: utf-8 -*-
"""
File  diagnostics.py
 
Author: Dr. Ernesto Adorio
           University of the Philippines, Clarkfield
           Angeles City, Pampanga.
Version 0.0.1 december 25 release
License: noncommercial academic use only
"""
from matlib import *
from lm     import *
 
# add the following simple function to lm.py:
def hatmatrix (X):
    """
    Computes the hat matrix.
    """
    VC= matinv(matxtx(X))
    B = matmatt(VC, X)
    return matmul(X, B) 
 
def diaghatmatrix(X):
    """
    returns diagonal elements of hatmatrix as a vector.
    """
    VC  = matinv(matxtx(X))
    hii = [0] * len(X)
    for i in range(len(X)):
       row = X[i][:]
       hii[i] = matvectmatvec(row, VC)
    return hii      
 
 
def VIF(X):
    """
    Returns the vector of VIF for each independent variable.
    X is the model matrix and may include leading constant column
    of all 1s.
 
    revised: dec.25,2009 tested with Neter's example.
    """
    rows = len(X)
    cols = len(X[0])
 
    #Check if first column of X is a column of all ones.
    startcol = 0
    for x in X:
       if x[0]-1.0 < 1.0e-5:
         startcol = 1
         break
 
    # Compute ALL VIFs.
    out = []
    for j in range(startcol, cols):
       Z   = matDelCols(X, [j])
       Xj  = [x[j] for x in X]
       fit = ols(Z, Xj)
       out.append(1.0/(1.0 - fit.R2))
    return out
 
def studdelresids(X,Y, hii):
    """
    Computes the studentized deleted residuals for model matrix X.
 
    Reference formula 11.25 of Neter, Kutner and Wassermann, page 400
    n -number of data
    p - number of columns of X
    hii -leverage values (diagonal of hat matrix) 
    ei - error at ith observation.
    d_i* = e_i [(n-p-1)/{SSE(1-hii) - e_i^2}]**(1/2)
    """
    n     = len(X)
    fit   = ols(X,Y)
    dstar = fit.resids
    if hii is None:
       hii  = diaghatmatrix(X)
    SSE   = fit.RSS
    num = n -len(X[0])-1
    for i in range(n):
       dstar[i] *= sqrt(num/(SSE*(1- hii[i]) - (dstar[i])**2))
    return dstar
 
 
def dffits(X,Y, hii):
    """
    Computes the dffits vector for model matrix X.
    Reference formula 11.25 of Neter, Kutner and Wassermann, page 400
    n -number of data
    p - number of variables
    hii -leverage values (diagonal of hat matrix) 
    ei - error at ith observation.
    d_i* = e_i [(n-p-1)/{SSE(1-hii) - e_i^2}]**(1/2)
    dffits = d_i* * ( hii/ (1 -hii))
    """
    n     = len(X)
    fit   = ols(X,Y)
    dstar = fit.resids
    if hii is None:
       hii  = diaghatmatrix(X)
    SSE   = fit.RSS
    num = n -len(X[0])-1
    for i in range(n):
       dstar[i] *= sqrt(num/(SSE*(1- hii[i]) - (dstar[i])**2))
       dstar[i] *= sqrt(hii[i]/ (1 -hii[i]))
    return dstar
 
def CookDistance(X,Y, hii=None,resids=None,MRSS = None):
    """
    Computes Cook's Distance using formula 11.29a p. 403
    """
    n = len(X)
    p = len(X[0])
    fit = None
    if resids is None:
       fit = ols(X,Y)
       ei  = fit.resids
    else:
       ei = resids
 
    if hii is None:
       hii  = diaghatmatrix(X)
    if MRSS is None:
       if fit is None:
          fit = ols(X,Y)
       MRSS = fit.MRSS
 
 
    Di = [0] * n
    invden = 1.0 / (p * MRSS)
    for i in range(n):
        Di[i] = ei[i]**2 * invden * hii[i]/ (1 - hii[i])**2   
    return Di
 
 
 
def Test():
   print "reading bodyfat data file body.fat"
   # Read data file "bodyfat.dat"
   sdata = open("data/bodyfat.dat").read()
   data = matreadstring(sdata)
   print"splitting matrix:"
   X, Y = matsplitbycolumn(data, -1)
   mataugprint(X, Y)
   """
   Alternately you can use:
   X = matDelCols(data, [-1])
   Y = [d[-1] for d in data]
   print "Input matrix"
   mataugprint(X,Y)
   """
 
   Xmodel = matInsertConstCol(X, 0, 1,inplace=False)
   print "Variance Inflation factors for full model"
   print "matches with p. 410, Table  11.4  of Neter's text."
   print VIF(Xmodel)
 
 
   # Neters's first two variable example:
   #print "Hat matrix for two variable example of Neter."
   Xmodel = matSelectCols(Xmodel, [0,1,2])
   fit = ols(Xmodel, Y)
   resids = fit.resids
   MRSS = fit.MRSS
   hii =   diaghatmatrix(Xmodel)
   dstars = studdelresids(Xmodel, Y, hii)
   DFFITS = dffits(Xmodel, Y, hii)
   D = CookDistance(X,Y, resids= resids, hii = hii, MRSS = MRSS)
   print
   print "Diagnostics for Neter's two variable model example. See pp. "
   print "%2s %8s %8s %8s %8s %8s" % ("i", "ei", "dffit", "dstar", "hii", "Cook")
   for i, (e,dffit, dstar,h,d) in enumerate(zip(resids,DFFITS,dstars,hii,D)):
      print "%2d %8.4f %8.4f %8.4f %8.4f %8.4f" %(i,e,dffit, dstar,h,d)
 
 
if __name__ == "__main__":
   Test()