src/lsp/regress.lsp

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