regression.lsp

   1 ;;; -*- mode: lisp -*-
   2 ;;;
   3 ;;; Copyright (c) 2005--2007, by A.J. Rossini <blindglobe@gmail.com>
   4 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
   5 ;;; Since 1991, ANSI was finally finished.  Modified to match ANSI
   6 ;;; Common Lisp.
   7
   8 ;;;;
   9 ;;;; regression.lsp XLISP-STAT regression model proto and methods
  10 ;;;; XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney
  11 ;;;; Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz
  12 ;;;; You may give out copies of this software; for conditions see the file
  13 ;;;; COPYING included with this distribution.
  14 ;;;;
  15 ;;;;
  16 ;;;; Incorporates modifications suggested by Sandy Weisberg.
  17 ;;;;
  18
  19 (in-package :cl-user)
  20
  21 (defpackage :lisp-stat-regression-linear
  22   (:use :common-lisp
  23         :lisp-stat-object-system
  24         :lisp-stat-basics
  25         :lisp-stat-compound-data
  26         :lisp-stat-sequence
  27         :lisp-stat-matrix)
  28   (:shadowing-import-from :lisp-stat-object-system
  29                           slot-value call-method call-next-method)
  30
  31   (:export regression-model regression-model-proto x y intercept sweep-matrix
  32            basis weights included total-sum-of-squares residual-sum-of-squares
  33            predictor-names response-name case-labels))
  34
  35 (in-package :lisp-stat-regression-linear)
  36
  37 ;;; Regresion Model Prototype
  38
  39 (defproto regression-model-proto
  40           '(x y intercept sweep-matrix basis weights
  41               included
  42               total-sum-of-squares
  43               residual-sum-of-squares
  44               predictor-names
  45               response-name
  46               case-labels)
  47           ()
  48           *object*
  49           "Normal Linear Regression Model")
  50
  51 (defun regression-model (x y &key
  52                            (intercept T)
  53                            (print T)
  54                            weights
  55                            (included (repeat t (length y)))
  56                            predictor-names
  57                            response-name
  58                            case-labels)
  59 "Args: (x y &key (intercept T) (print T) weights
  60           included predictor-names response-name case-labels)
  61 X           - list of independent variables or X matrix
  62 Y           - dependent variable.
  63 INTERCEPT   - T to include (default), NIL for no intercept
  64 PRINT       - if not NIL print summary information
  65 WEIGHTS     - if supplied should be the same length as Y; error variances are
  66               assumed to be inversely proportional to WEIGHTS
  67 PREDICTOR-NAMES, RESPONSE-NAME, CASE-LABELS
  68             - sequences of strings or symbols.
  69 INCLUDED    - if supplied should be the same length as Y, with elements nil
  70               to skip a in computing estimates (but not in residual analysis).
  71 Returns a regression model object. To examine the model further assign the
  72 result to a variable and send it messages.
  73 Example (data are in file absorbtion.lsp in the sample data directory/folder):
  74   (def m (regression-model (list iron aluminum) absorbtion))
  75   (send m :help) (send m :plot-residuals)"
  76   (let ((x (cond
  77             ((matrixp x) x)
  78             ((vectorp x) (list x))
  79             ((and (consp x) (numberp (car x))) (list x))
  80             (t x)))
  81         (m (send regression-model-proto :new)))
  82     (send m :x (if (matrixp x) x (apply #'bind-columns x)))
  83     (send m :y y)
  84     (send m :intercept intercept)
  85     (send m :weights weights)
  86     (send m :included included)
  87     (send m :predictor-names predictor-names)
  88     (send m :response-name response-name)
  89     (send m :case-labels case-labels)
  90     (if print (send m :display))
  91     m))
  92
  93 (defmeth regression-model-proto :isnew ()
  94   (send self :needs-computing t))
  95
  96 (defmeth regression-model-proto :save ()
  97 "Message args: ()
  98 Returns an expression that will reconstruct the regression model."
  99   `(regression-model ',(send self :x)
 100                      ',(send self :y)
 101                      :intercept ',(send self :intercept)
 102                      :weights ',(send self :weights)
 103                      :included ',(send self :included)
 104                      :predictor-names ',(send self :predictor-names)
 105                      :response-name ',(send self :response-name)
 106                      :case-labels ',(send self :case-labels)))
 107
 108 ;;; Computing and Display Methods
 109
 110 (defmeth regression-model-proto :compute ()
 111 "Message args: ()
 112 Recomputes the estimates. For internal use by other messages"
 113   (let* ((included (if-else (send self :included) 1 0))
 114          (x (send self :x))
 115          (y (send self :y))
 116          (intercept (send self :intercept))
 117          (weights (send self :weights))
 118          (w (if weights (* included weights) included))
 119          (m (make-sweep-matrix x y w))
 120          (n (array-dimension x 1))
 121          (p (- (array-dimension m 0) 1))
 122          (tss (aref m p p))
 123          (tol (* 0.001 (reduce #'* (mapcar #'standard-deviation (column-list x)))))
 124          ;; (tol (* 0.001 (apply #'* (mapcar #'standard-deviation (column-list x)))))
 125          (sweep-result
 126           (if intercept
 127               (sweep-operator m (iseq 1 n) tol)
 128               (sweep-operator m (iseq 0 n) (cons 0.0 tol)))))
 129     (setf (slot-value 'sweep-matrix) (first sweep-result))
 130     (setf (slot-value 'total-sum-of-squares) tss)
 131     (setf (slot-value 'residual-sum-of-squares)
 132           (aref (first sweep-result) p p))
 133     (setf (slot-value 'basis)
 134           (let ((b (remove 0 (second sweep-result))))
 135             (if b (- (reduce #'- (reverse b)) 1)
 136                 (error "no columns could be swept"))))))
 137
 138 (defmeth regression-model-proto :needs-computing (&optional set)
 139   (if set (setf (slot-value 'sweep-matrix) nil))
 140   (null (slot-value 'sweep-matrix)))
 141
 142 (defmeth regression-model-proto :display ()
 143 "Message args: ()
 144 Prints the least squares regression summary. Variables not used in the fit
 145 are marked as aliased."
 146   (let ((coefs (coerce (send self :coef-estimates) 'list))
 147         (se-s (send self :coef-standard-errors))
 148         (x (send self :x))
 149         (p-names (send self :predictor-names)))
 150     (if (send self :weights)
 151         (format t "~%Weighted Least Squares Estimates:~2%")
 152         (format t "~%Least Squares Estimates:~2%"))
 153     (when (send self :intercept)
 154           (format t "Constant               ~10f   ~A~%"
 155                   (car coefs) (list (car se-s)))
 156           (setf coefs (cdr coefs))
 157           (setf se-s (cdr se-s)))
 158     (dotimes (i (array-dimension x 1))
 159              (cond
 160                ((member i (send self :basis))
 161                 (format t "~22a ~10f   ~A~%"
 162                         (select p-names i) (car coefs) (list (car se-s)))
 163                 (setf coefs (cdr coefs) se-s (cdr se-s)))
 164                (t (format t "~22a    aliased~%" (select p-names i)))))
 165     (format t "~%")
 166     (format t "R Squared:             ~10f~%" (send self :r-squared))
 167     (format t "Sigma hat:             ~10f~%" (send self :sigma-hat))
 168     (format t "Number of cases:       ~10d~%" (send self :num-cases))
 169     (if (/= (send self :num-cases) (send self :num-included))
 170         (format t "Number of cases used:  ~10d~%" (send self :num-included)))
 171     (format t "Degrees of freedom:    ~10d~%" (send self :df))
 172     (format t "~%")))
 173
 174 ;;; Slot accessors and mutators
 175
 176 (defmeth regression-model-proto :x (&optional new-x)
 177 "Message args: (&optional new-x)
 178 With no argument returns the x matrix as supplied to m. With an argument
 179 NEW-X sets the x matrix to NEW-X and recomputes the estimates."
 180   (when (and new-x (matrixp new-x))
 181         (setf (slot-value 'x) new-x)
 182         (send self :needs-computing t))
 183   (slot-value 'x))
 184
 185 (defmeth regression-model-proto :y (&optional new-y)
 186 "Message args: (&optional new-y)
 187 With no argument returns the y sequence as supplied to m. With an argument
 188 NEW-Y sets the y sequence to NEW-Y and recomputes the estimates."
 189   (when (and new-y (or (matrixp new-y) (sequencep new-y)))
 190         (setf (slot-value 'y) new-y)
 191         (send self :needs-computing t))
 192   (slot-value 'y))
 193
 194 (defmeth regression-model-proto :intercept (&optional (val nil set))
 195 "Message args: (&optional new-intercept)
 196 With no argument returns T if the model includes an intercept term, nil if
 197 not. With an argument NEW-INTERCEPT the model is changed to include or
 198 exclude an intercept, according to the value of NEW-INTERCEPT."
 199   (when set
 200         (setf (slot-value 'intercept) val)
 201         (send self :needs-computing t))
 202   (slot-value 'intercept))
 203
 204 (defmeth regression-model-proto :weights (&optional (new-w nil set))
 205 "Message args: (&optional new-w)
 206 With no argument returns the weight sequence as supplied to m; NIL means
 207 an unweighted model. NEW-W sets the weights sequence to NEW-W and
 208 recomputes the estimates."
 209   (when set
 210         (setf (slot-value 'weights) new-w)
 211         (send self :needs-computing t))
 212   (slot-value 'weights))
 213
 214 (defmeth regression-model-proto :total-sum-of-squares ()
 215 "Message args: ()
 216 Returns the total sum of squares around the mean."
 217   (if (send self :needs-computing) (send self :compute))
 218   (slot-value 'total-sum-of-squares))
 219
 220 (defmeth regression-model-proto :residual-sum-of-squares ()
 221 "Message args: ()
 222 Returns the residual sum of squares for the model."
 223   (if (send self :needs-computing) (send self :compute))
 224   (slot-value 'residual-sum-of-squares))
 225
 226 (defmeth regression-model-proto :basis ()
 227 "Message args: ()
 228 Returns the indices of the variables used in fitting the model."
 229   (if (send self :needs-computing) (send self :compute))
 230   (slot-value 'basis))
 231
 232 (defmeth regression-model-proto :sweep-matrix ()
 233 "Message args: ()
 234 Returns the swept sweep matrix. For internal use"
 235   (if (send self :needs-computing) (send self :compute))
 236   (slot-value 'sweep-matrix))
 237
 238 (defmeth regression-model-proto :included (&optional new-included)
 239 "Message args: (&optional new-included)
 240 With no argument,  NIL means a case is not used in calculating estimates, and non-nil means it is used.  NEW-INCLUDED is a sequence of length of y of nil and t to select cases.  Estimates are recomputed."
 241   (when (and new-included
 242              (= (length new-included) (send self :num-cases)))
 243         (setf (slot-value 'included) (copy-seq new-included))
 244         (send self :needs-computing t))
 245   (if (slot-value 'included)
 246       (slot-value 'included)
 247       (repeat t (send self :num-cases))))
 248
 249 (defmeth regression-model-proto :predictor-names (&optional (names nil set))
 250 "Message args: (&optional (names nil set))
 251 With no argument returns the predictor names. NAMES sets the names."
 252   (if set (setf (slot-value 'predictor-names) (mapcar #'string names)))
 253   (let ((p (array-dimension (send self :x) 1))
 254         (p-names (slot-value 'predictor-names)))
 255     (if (not (and p-names (= (length p-names) p)))
 256         (setf (slot-value 'predictor-names)
 257               (mapcar #'(lambda (a) (format nil "Variable ~a" a))
 258                       (iseq 0 (- p 1))))))
 259   (slot-value 'predictor-names))
 260
 261 (defmeth regression-model-proto :response-name (&optional (name "Y" set))
 262 "Message args: (&optional name)
 263 With no argument returns the response name. NAME sets the name."
 264   (if set (setf (slot-value 'response-name) (if name (string name) "Y")))
 265   (slot-value 'response-name))
 266
 267 (defmeth regression-model-proto :case-labels (&optional (labels nil set))
 268 "Message args: (&optional labels)
 269 With no argument returns the case-labels. LABELS sets the labels."
 270   (if set (setf (slot-value 'case-labels)
 271                 (if labels
 272                     (mapcar #'string labels)
 273                     (mapcar #'(lambda (x) (format nil "~d" x))
 274                             (iseq 0 (- (send self :num-cases) 1))))))
 275   (slot-value 'case-labels))
 276
 277 ;;;
 278 ;;; Other Methods
 279 ;;; None of these methods access any slots directly.
 280 ;;;
 281
 282 (defmeth regression-model-proto :num-cases ()
 283 "Message args: ()
 284 Returns the number of cases in the model."
 285   (length (send self :y)))
 286
 287 (defmeth regression-model-proto :num-included ()
 288 "Message args: ()
 289 Returns the number of cases used in the computations."
 290   (sum (if-else (send self :included) 1 0)))
 291
 292 (defmeth regression-model-proto :num-coefs ()
 293 "Message args: ()
 294 Returns the number of coefficients in the fit model (including the
 295 intercept if the model includes one)."
 296   (if (send self :intercept)
 297       (+ 1 (length (send self :basis)))
 298       (length (send self :basis))))
 299
 300 (defmeth regression-model-proto :df ()
 301 "Message args: ()
 302 Returns the number of degrees of freedom in the model."
 303   (- (send self :num-included) (send self :num-coefs)))
 304
 305 (defmeth regression-model-proto :x-matrix ()
 306 "Message args: ()
 307 Returns the X matrix for the model, including a column of 1's, if
 308 appropriate. Columns of X matrix correspond to entries in basis."
 309   (let ((m (select (send self :x)
 310                    (iseq 0 (- (send self :num-cases) 1))
 311                    (send self :basis))))
 312     (if (send self :intercept)
 313         (bind-columns (repeat 1 (send self :num-cases)) m)
 314         m)))
 315
 316 (defmeth regression-model-proto :leverages ()
 317 "Message args: ()
 318 Returns the diagonal elements of the hat matrix."
 319   (let* ((weights (send self :weights))
 320          (x (send self :x-matrix))
 321          (raw-levs
 322           (matmult (* (matmult x (send self :xtxinv)) x)
 323                    (repeat 1 (send self :num-coefs)))))
 324     (if weights (* weights raw-levs) raw-levs)))
 325
 326 (defmeth regression-model-proto :fit-values ()
 327 "Message args: ()
 328 Returns the fitted values for the model."
 329   (matmult (send self :x-matrix) (send self :coef-estimates)))
 330
 331 (defmeth regression-model-proto :raw-residuals ()
 332 "Message args: ()
 333 Returns the raw residuals for a model."
 334   (- (send self :y) (send self :fit-values)))
 335
 336 (defmeth regression-model-proto :residuals ()
 337 "Message args: ()
 338 Returns the raw residuals for a model without weights. If the model
 339 includes weights the raw residuals times the square roots of the weights
 340 are returned."
 341   (let ((raw-residuals (send self :raw-residuals))
 342         (weights (send self :weights)))
 343     (if weights (* (sqrt weights) raw-residuals) raw-residuals)))
 344
 345 (defmeth regression-model-proto :sum-of-squares ()
 346 "Message args: ()
 347 Returns the error sum of squares for the model."
 348   (send self :residual-sum-of-squares))
 349
 350 (defmeth regression-model-proto :sigma-hat ()
 351 "Message args: ()
 352 Returns the estimated standard deviation of the deviations about the
 353 regression line."
 354   (let ((ss (send self :sum-of-squares))
 355         (df (send self :df)))
 356     (if (/= df 0) (sqrt (/ ss df)))))
 357
 358 ;; for models without an intercept the 'usual' formula for R^2 can give
 359 ;; negative results; hence the max.
 360 (defmeth regression-model-proto :r-squared ()
 361 "Message args: ()
 362 Returns the sample squared multiple correlation coefficient, R squared, for
 363 the regression."
 364   (max (- 1 (/ (send self :sum-of-squares) (send self :total-sum-of-squares)))
 365        0))
 366
 367 (defmeth regression-model-proto :coef-estimates ()
 368 "Message args: ()
 369 Returns the OLS (ordinary least squares) estimates of the regression
 370 coefficients. Entries beyond the intercept correspond to entries in basis."
 371   (let ((n (array-dimension (send self :x) 1))
 372         (indices (if (send self :intercept)
 373                      (cons 0 (+ 1 (send self :basis)))
 374                      (+ 1 (send self :basis))))
 375         (m (send self :sweep-matrix)))
 376     (coerce (compound-data-seq (select m (+ 1 n) indices)) 'list)))
 377
 378 (defmeth regression-model-proto :xtxinv ()
 379 "Message args: ()
 380 Returns ((X^T) X)^(-1) or ((X^T) W X)^(-1)."
 381   (let ((indices (if (send self :intercept)
 382                      (cons 0 (1+ (send self :basis)))
 383                      (1+ (send self :basis)))))
 384     (select (send self :sweep-matrix) indices indices)))
 385
 386 (defmeth regression-model-proto :coef-standard-errors ()
 387 "Message args: ()
 388 Returns estimated standard errors of coefficients. Entries beyond the
 389 intercept correspond to entries in basis."
 390   (let ((s (send self :sigma-hat)))
 391     (if s (* (send self :sigma-hat) (sqrt (diagonal (send self :xtxinv)))))))
 392
 393 (defmeth regression-model-proto :studentized-residuals ()
 394 "Message args:  ()
 395 Computes the internally studentized residuals for included cases and externally studentized residuals for excluded cases."
 396   (let ((res (send self :residuals))
 397         (lev (send self :leverages))
 398         (sig (send self :sigma-hat))
 399         (inc (send self :included)))
 400     (if-else inc
 401              (/ res (* sig (sqrt (pmax .00001 (- 1 lev)))))
 402              (/ res (* sig (sqrt (+ 1 lev)))))))
 403
 404 (defmeth regression-model-proto :externally-studentized-residuals ()
 405 "Message args:  ()
 406 Computes the externally studentized residuals."
 407   (let* ((res (send self :studentized-residuals))
 408          (df (send self :df)))
 409     (if-else (send self :included)
 410              (* res (sqrt (/ (- df 1) (- df (^ res 2)))))
 411              res)))
 412
 413 (defmeth regression-model-proto :cooks-distances ()
 414 "Message args: ()
 415 Computes Cook's distances."
 416   (let ((lev (send self :leverages))
 417         (res (/ (^ (send self :studentized-residuals) 2)
 418                 (send self :num-coefs))))
 419     (if-else (send self :included) (* res (/ lev (- 1 lev) )) (* res lev))))
 420
 421 (defmeth regression-model-proto :plot-residuals (&optional x-values)
 422 "Message args: (&optional x-values)
 423 Opens a window with a plot of the residuals. If X-VALUES are not supplied
 424 the fitted values are used. The plot can be linked to other plots with the
 425 link-views function. Returns a plot object."
 426   (plot-points (if x-values x-values (send self :fit-values))
 427                (send self :residuals)
 428                :title "Residual Plot"
 429                :point-labels (send self :case-labels)))
 430
 431 (defmeth regression-model-proto :plot-bayes-residuals
 432   (&optional x-values)
 433 "Message args: (&optional x-values)
 434 Opens a window with a plot of the standardized residuals and two standard
 435 error bars for the posterior distribution of the actual deviations from the
 436 line. See Chaloner and Brant. If X-VALUES are not supplied  the fitted values
 437 are used. The plot can be linked to other plots with the link-views function.
 438 Returns a plot object."
 439   (let* ((r (/ (send self :residuals) (send self :sigma-hat)))
 440          (d (* 2 (sqrt (send self :leverages))))
 441          (low (- r d))
 442          (high (+ r d))
 443          (x-values (if x-values x-values (send self :fit-values)))
 444          (p (plot-points x-values r
 445                          :title "Bayes Residual Plot"
 446                          :point-labels (send self :case-labels))))
 447 ;; AJR:FIXME
 448 ;; the lambda needs to be something that fits into list
 449 ;;    (map 'list
 450 ;;       #'(lambda (a b c d) (send p :plotline a b c d nil))
 451 ;;       x-values low x-values high)
 452     (send p :adjust-to-data)
 453     p))