src/stat-models/regression-sheeple.lisp

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