lisp-matrix direct name conversions.
[CommonLispStat.git] / src / stat-models / regression.lisp
blob4b0923fe424c72e5f284486f1754805569c3cb03
1 ;;; -*- mode: lisp -*-
2 ;;;
3 ;;; Copyright (c) 2008--, 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.
8 ;;;; Originally from:
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 ;;;; Incorporates modifications suggested by Sandy Weisberg.
17 ;;; This version uses lisp-matrix for underlying numerics.
19 (in-package :lisp-stat-regression-linear)
21 ;;; Regresion Model Prototype
23 ;; The general strategy behind the fitting of models using prototypes
24 ;; is that we need to think about want the actual fits are, and then
25 ;; the fits can be used to recompute as components are changes. One
26 ;; catch here is that we'd like some notion of trace-ability, in
27 ;; particular, there is not necessarily a fixed way to take care of the
28 ;; audit trail. save-and-die might be a means of recording the final
29 ;; approach, but we are challenged by the problem of using advice and
30 ;; other such features to capture stages and steps that are considered
31 ;; along the goals of estimating a model.
33 ;; Note that the above is a stream-of-conscience response to the
34 ;; challenge of reproducibility in the setting of prototype "on-line"
35 ;; computation.
37 (defvar regression-model-proto nil
38 "Prototype for all regression model instances.")
40 (defproto regression-model-proto
41 '(x y intercept sweep-matrix basis weights
42 included
43 total-sum-of-squares
44 residual-sum-of-squares
45 predictor-names
46 response-name
47 case-labels
48 doc)
50 *object*
51 "Normal Linear Regression Model")
54 (defun regression-model
55 (x y &key
56 (intercept T)
57 (print T)
58 (weights nil)
59 (included (repeat t (vector-dimension y)))
60 predictor-names
61 response-name
62 case-labels
63 (doc "Undocumented Regression Model Instance")
64 (debug T))
65 "Args: (x y &key (intercept T) (print T) (weights nil)
66 included predictor-names response-name case-labels)
67 X - list of independent variables or X matrix
68 Y - dependent variable.
69 INTERCEPT - T to include (default), NIL for no intercept
70 PRINT - if not NIL print summary information
71 WEIGHTS - if supplied should be the same length as Y; error
72 variances are
73 assumed to be inversely proportional to WEIGHTS
74 PREDICTOR-NAMES, RESPONSE-NAME, CASE-LABELS
75 - sequences of strings or symbols.
76 INCLUDED - if supplied should be the same length as Y, with
77 elements nil to skip a in computing estimates (but not
78 in residual analysis).
79 Returns a regression model object. To examine the model further assign the
80 result to a variable and send it messages.
81 Example (data are in file absorbtion.lsp in the sample data directory):
82 (def m (regression-model (list iron aluminum) absorbtion))
83 (send m :help) (send m :plot-residuals)"
84 (let ((x (cond
85 ((typep x 'matrix-like) x)
86 #| assume only numerical vectors -- but we need to ensure coercion to float.
87 ((or (typep x 'sequence)
88 (and (consp x)
89 (numberp (car x)))
90 (make-vector (length x) :initial-contents x)))
92 (t (error "not matrix-like.");x
93 ))) ;; actually, might should barf.
94 (y (cond
95 ((typep y 'vector-like) y)
97 ((and (consp x)
98 (numberp (car x))) (make-vector (length y) :initial-contents y))
100 (t (error "not vector-like."); y
101 ))) ;; actually, might should barf.
102 (m (send regression-model-proto :new)))
103 (format t "~%")
104 (send m :doc doc)
105 (send m :x x)
106 (send m :y y)
107 (send m :intercept intercept)
108 (send m :weights weights)
109 (send m :included included)
110 (send m :predictor-names predictor-names)
111 (send m :response-name response-name)
112 (send m :case-labels case-labels)
113 (if debug
114 (progn
115 (format t "~%")
116 (format t "~S~%" (send m :doc))
117 (format t "X: ~S~%" (send m :x))
118 (format t "Y: ~S~%" (send m :y))))
119 (if print (send m :display))
123 ;; regression-model is the old API, but regression as a generic will
124 ;; be the new API. We need to distinguish between APIs which enable
125 ;; the user to do clear activities, and APIs which enable developers
126 ;; to do clear extensions and development, and underlying
127 ;; infrastructure to keep everything straight and enabled.
129 ;; There are conflicting theories for how to structure the
130 ;; specification of mathematical models, along with the statistical
131 ;; inference, along with the data which is instantiating the model.
133 ;; i.e.: mathematical model for the relationships between components,
134 ;; between a component and a summarizing parameter, and between
135 ;; parameters.
137 ;; statistical inference describes the general approach for
138 ;; aggregating into a decision and has impliciations for the scale up
139 ;; from the model on a single instance to the generalization.
141 ;; The data represents the particular substantive context that is
142 ;; driving the model/inference combination, and about which we hope to
143 ;; generate knowledge.
145 ;; numerical analysis selects appropriate algorithms/implementations
146 ;; for combining the above 3.
148 ;; the end result is input on the decision being made (which could be
149 ;; specific (decision analysis/testing), risk-analysis (interval
150 ;; estimation) , most likely/appropriate selection (point estimation)
154 (defclass model ()
155 ((type structure)))
157 (defgeneric regression ;; assumes x/y from lisp-matrix -- start of a set of generics.
158 (model dataset)
159 "Args: (x y &key (intercept T) (print T) (weights nil)
160 included predictor-names response-name case-labels)
161 X - list of independent variables or X matrix
162 Y - dependent variable.
163 INTERCEPT - T to include (default), NIL for no intercept
164 PRINT - if not NIL print summary information
165 WEIGHTS - if supplied should be the same length as Y; error
166 variances are
167 assumed to be inversely proportional to WEIGHTS
168 PREDICTOR-NAMES, RESPONSE-NAME, CASE-LABELS
169 - sequences of strings or symbols.
170 INCLUDED - if supplied should be the same length as Y, with
171 elements nil to skip a in computing estimates (but not
172 in residual analysis).
173 Returns a regression model object. To examine the model further assign the
174 result to a variable and send it messages.
175 Example (data are in file absorbtion.lsp in the sample data directory):
176 (def m (regression-model (list iron aluminum) absorbtion))
177 (send m :help) (send m :plot-residuals)"
178 (let ((m (send regression-model-proto :new)))
179 (format t "~%")
180 (send m :doc doc)
181 (send m :x x)
182 (send m :y y)
183 (send m :intercept intercept)
184 (send m :weights weights)
185 (send m :included included)
186 (send m :predictor-names predictor-names)
187 (send m :response-name response-name)
188 (send m :case-labels case-labels)
189 (if debug
190 (progn
191 (format t "~%")
192 (format t "~S~%" (send m :doc))
193 (format t "X: ~S~%" (send m :x))
194 (format t "Y: ~S~%" (send m :y))))
195 (if print (send m :display))
200 (defmeth regression-model-proto :isnew ()
201 (send self :needs-computing t))
203 (defmeth regression-model-proto :save ()
204 "Message args: ()
205 Returns an expression that will reconstruct the regression model."
206 `(regression-model ',(send self :x)
207 ',(send self :y)
208 :intercept ',(send self :intercept)
209 :weights ',(send self :weights)
210 :included ',(send self :included)
211 :predictor-names ',(send self :predictor-names)
212 :response-name ',(send self :response-name)
213 :case-labels ',(send self :case-labels)))
215 ;;; Computing and Display Methods
217 (defmeth regression-model-proto :compute ()
218 "Message args: ()
219 Recomputes the estimates. For internal use by other messages"
220 (let* ((included (if-else (send self :included) 1 0))
221 (x (send self :x))
222 (y (send self :y))
223 (intercept (send self :intercept)) ;; T/nil
224 (weights (send self :weights)) ;; vector-like
225 (w (if weights (* included weights) included))
226 (m (make-sweep-matrix x y w)) ;;; ERROR HERE of course!
227 (n (matrix-dimension x 1))
228 (p (- (matrix-dimension m 0) 1)) ;; remove intercept from # params -- right?
229 (tss (mref m p p))
230 (tol (* 0.001
231 (reduce #'* (mapcar #'standard-deviation
232 (list-of-columns x)))))
233 (sweep-result
234 (if intercept
235 (sweep-operator m (iseq 1 n) tol)
236 (sweep-operator m (iseq 0 n) (cons 0.0 tol)))))
237 (format t
238 "~%REMOVEME: regr-mdl-prto :compute~%Sweep= ~A~%x= ~A~%y= ~A~%m= ~A~%tss= ~A~%"
239 sweep-result x y m tss)
240 (setf (slot-value 'sweep-matrix) (first sweep-result))
241 (setf (slot-value 'total-sum-of-squares) tss)
242 (setf (slot-value 'residual-sum-of-squares)
243 (mref (first sweep-result) p p))
244 ;; SOMETHING WRONG HERE! FIX-ME
245 (setf (slot-value 'basis)
246 (let ((b (remove 0 (second sweep-result))))
247 (if b (- (reduce #'- (reverse b)) 1)
248 (error "no columns could be swept"))))))
250 (defmeth regression-model-proto :needs-computing (&optional set)
251 "Message args: ( &optional set )
253 If value given, sets the flag for whether (re)computation is needed to
254 update the model fits."
255 (send self :nop)
256 (if set (setf (slot-value 'sweep-matrix) nil))
257 (null (slot-value 'sweep-matrix)))
259 (defmeth regression-model-proto :display ()
260 "Message args: ()
262 Prints the least squares regression summary. Variables not used in the fit
263 are marked as aliased."
264 (let ((coefs (coerce (send self :coef-estimates) 'list))
265 (se-s (send self :coef-standard-errors))
266 (x (send self :x))
267 (p-names (send self :predictor-names)))
268 (if (send self :weights)
269 (format t "~%Weighted Least Squares Estimates:~2%")
270 (format t "~%Least Squares Estimates:~2%"))
271 (when (send self :intercept)
272 (format t "Constant ~10f ~A~%"
273 (car coefs) (list (car se-s)))
274 (setf coefs (cdr coefs))
275 (setf se-s (cdr se-s)))
276 (dotimes (i (array-dimension x 1))
277 (cond
278 ((member i (send self :basis))
279 (format t "~22a ~10f ~A~%"
280 (select p-names i) (car coefs) (list (car se-s)))
281 (setf coefs (cdr coefs) se-s (cdr se-s)))
282 (t (format t "~22a aliased~%" (select p-names i)))))
283 (format t "~%")
284 (format t "R Squared: ~10f~%" (send self :r-squared))
285 (format t "Sigma hat: ~10f~%" (send self :sigma-hat))
286 (format t "Number of cases: ~10d~%" (send self :num-cases))
287 (if (/= (send self :num-cases) (send self :num-included))
288 (format t "Number of cases used: ~10d~%" (send self :num-included)))
289 (format t "Degrees of freedom: ~10d~%" (send self :df))
290 (format t "~%")))
292 ;;; Slot accessors and mutators
294 (defmeth regression-model-proto :doc (&optional new-doc append)
295 "Message args: (&optional new-doc)
297 Returns the DOC-STRING as supplied to m.
298 Additionally, with an argument NEW-DOC, sets the DOC-STRING to
299 NEW-DOC. In this setting, when APPEND is T, don't replace and just
300 append NEW-DOC to DOC."
301 (send self :nop)
302 (when (and new-doc (stringp new-doc))
303 (setf (slot-value 'doc)
304 (if append
305 (concatenate 'string
306 (slot-value 'doc)
307 new-doc)
308 new-doc)))
309 (slot-value 'doc))
312 (defmeth regression-model-proto :x (&optional new-x)
313 "Message args: (&optional new-x)
315 With no argument returns the x matrix-like as supplied to m. With an
316 argument, NEW-X sets the x matrix-like to NEW-X and recomputes the
317 estimates."
318 (when (and new-x (typep new-x 'matrix-like))
319 (setf (slot-value 'x) new-x)
320 (send self :needs-computing t))
321 (slot-value 'x))
323 (defmeth regression-model-proto :y (&optional new-y)
324 "Message args: (&optional new-y)
326 With no argument returns the y vector-like as supplied to m. With an
327 argument, NEW-Y sets the y vector-like to NEW-Y and recomputes the
328 estimates."
329 (when (and new-y
330 (typep new-y 'vector-like))
331 (setf (slot-value 'y) new-y) ;; fixme -- pls set slot value to a vector-like!
332 (send self :needs-computing t))
333 (slot-value 'y))
335 (defmeth regression-model-proto :intercept (&optional (val nil set))
336 "Message args: (&optional new-intercept)
338 With no argument returns T if the model includes an intercept term,
339 nil if not. With an argument NEW-INTERCEPT the model is changed to
340 include or exclude an intercept, according to the value of
341 NEW-INTERCEPT."
342 (when set
343 (setf (slot-value 'intercept) val)
344 (send self :needs-computing t))
345 (slot-value 'intercept))
347 (defmeth regression-model-proto :weights (&optional (new-w nil set))
348 "Message args: (&optional new-w)
350 With no argument returns the weight vector-like as supplied to m; NIL
351 means an unweighted model. NEW-W sets the weights vector-like to NEW-W
352 and recomputes the estimates."
353 (when (and set
354 (typep new-w 'vector-like))
355 (setf (slot-value 'weights) new-w)
356 (send self :needs-computing t))
357 (slot-value 'weights))
359 (defmeth regression-model-proto :total-sum-of-squares ()
360 "Message args: ()
362 Returns the total sum of squares around the mean.
363 This is recomputed if an update is needed."
364 (if (send self :needs-computing)
365 (send self :compute))
366 (slot-value 'total-sum-of-squares))
368 (defmeth regression-model-proto :residual-sum-of-squares ()
369 "Message args: ()
371 Returns the residual sum of squares for the model.
372 This is recomputed if an update is needed."
373 (if (send self :needs-computing)
374 (send self :compute))
375 (slot-value 'residual-sum-of-squares))
377 (defmeth regression-model-proto :basis ()
378 "Message args: ()
380 Returns the indices of the variables used in fitting the model, in a
381 sequence.
382 This is recomputed if an update is needed."
383 (if (send self :needs-computing)
384 (send self :compute))
385 (slot-value 'basis))
388 (defmeth regression-model-proto :sweep-matrix ()
389 "Message args: ()
391 Returns the swept sweep matrix. For internal use"
392 (if (send self :needs-computing)
393 (send self :compute))
394 (slot-value 'sweep-matrix))
396 (defmeth regression-model-proto :included (&optional new-included)
397 "Message args: (&optional new-included)
399 With no argument, NIL means a case is not used in calculating
400 estimates, and non-nil means it is used. NEW-INCLUDED is a sequence
401 of length of y of nil and t to select cases. Estimates are
402 recomputed."
403 (when (and new-included
404 (= (length new-included) (send self :num-cases)))
405 (setf (slot-value 'included) (copy-seq new-included))
406 (send self :needs-computing t))
407 (if (slot-value 'included)
408 (slot-value 'included)
409 (repeat t (send self :num-cases))))
411 (defmeth regression-model-proto :predictor-names (&optional (names nil set))
412 "Message args: (&optional (names nil set))
414 With no argument returns the predictor names. NAMES sets the names."
415 (if set (setf (slot-value 'predictor-names) (mapcar #'string names)))
416 (let ((p (matrix-dimension (send self :x) 1))
417 (p-names (slot-value 'predictor-names)))
418 (if (not (and p-names (= (length p-names) p)))
419 (setf (slot-value 'predictor-names)
420 (mapcar #'(lambda (a) (format nil "Variable ~a" a))
421 (iseq 0 (- p 1))))))
422 (slot-value 'predictor-names))
424 (defmeth regression-model-proto :response-name (&optional (name "Y" set))
425 "Message args: (&optional name)
427 With no argument returns the response name. NAME sets the name."
428 (send self :nop)
429 (if set (setf (slot-value 'response-name) (if name (string name) "Y")))
430 (slot-value 'response-name))
432 (defmeth regression-model-proto :case-labels (&optional (labels nil set))
433 "Message args: (&optional labels)
434 With no argument returns the case-labels. LABELS sets the labels."
435 (if set (setf (slot-value 'case-labels)
436 (if labels
437 (mapcar #'string labels)
438 (mapcar #'(lambda (x) (format nil "~d" x))
439 (iseq 0 (- (send self :num-cases) 1))))))
440 (slot-value 'case-labels))
443 ;;; Other Methods
444 ;;; None of these methods access any slots directly.
447 (defmeth regression-model-proto :num-cases ()
448 "Message args: ()
449 Returns the number of cases in the model."
450 (nelts (send self :y)))
452 (defmeth regression-model-proto :num-included ()
453 "Message args: ()
454 Returns the number of cases used in the computations."
455 (sum (if-else (send self :included) 1 0)))
457 (defmeth regression-model-proto :num-coefs ()
458 "Message args: ()
459 Returns the number of coefficients in the fit model (including the
460 intercept if the model includes one)."
461 (if (send self :intercept)
462 (+ 1 (nelts (send self :basis)))
463 (nelts (send self :basis))))
465 (defmeth regression-model-proto :df ()
466 "Message args: ()
467 Returns the number of degrees of freedom in the model."
468 (- (send self :num-included) (send self :num-coefs)))
470 (defmeth regression-model-proto :x-matrix ()
471 "Message args: ()
472 Returns the X matrix for the model, including a column of 1's, if
473 appropriate. Columns of X matrix correspond to entries in basis."
474 (let ((m (select (send self :x)
475 (iseq 0 (- (send self :num-cases) 1))
476 (send self :basis))))
477 (if (send self :intercept)
478 (bind2 (repeat 1 (send self :num-cases)) m)
479 m)))
481 (defmeth regression-model-proto :leverages ()
482 "Message args: ()
483 Returns the diagonal elements of the hat matrix."
484 (let* ((weights (send self :weights))
485 (x (send self :x-matrix))
486 (raw-levs
487 (m* (* (m* x (send self :xtxinv)) x)
488 (repeat 1 (send self :num-coefs)))))
489 (if weights (* weights raw-levs) raw-levs)))
491 (defmeth regression-model-proto :fit-values ()
492 "Message args: ()
493 Returns the fitted values for the model."
494 (m* (send self :x-matrix) (send self :coef-estimates)))
496 (defmeth regression-model-proto :raw-residuals ()
497 "Message args: ()
498 Returns the raw residuals for a model."
499 (- (send self :y) (send self :fit-values)))
501 (defmeth regression-model-proto :residuals ()
502 "Message args: ()
503 Returns the raw residuals for a model without weights. If the model
504 includes weights the raw residuals times the square roots of the weights
505 are returned."
506 (let ((raw-residuals (send self :raw-residuals))
507 (weights (send self :weights)))
508 (if weights (* (sqrt weights) raw-residuals) raw-residuals)))
510 (defmeth regression-model-proto :sum-of-squares ()
511 "Message args: ()
512 Returns the error sum of squares for the model."
513 (send self :residual-sum-of-squares))
515 (defmeth regression-model-proto :sigma-hat ()
516 "Message args: ()
517 Returns the estimated standard deviation of the deviations about the
518 regression line."
519 (let ((ss (send self :sum-of-squares))
520 (df (send self :df)))
521 (if (/= df 0) (sqrt (/ ss df)))))
523 ;; for models without an intercept the 'usual' formula for R^2 can give
524 ;; negative results; hence the max.
525 (defmeth regression-model-proto :r-squared ()
526 "Message args: ()
527 Returns the sample squared multiple correlation coefficient, R squared, for
528 the regression."
529 (max (- 1 (/ (send self :sum-of-squares) (send self :total-sum-of-squares)))
532 (defmeth regression-model-proto :coef-estimates ()
533 "Message args: ()
535 Returns the OLS (ordinary least squares) estimates of the regression
536 coefficients. Entries beyond the intercept correspond to entries in
537 basis."
538 (let ((n (matrix-dimension (send self :x) 1))
539 (indices (flatten-list
540 (if (send self :intercept)
541 (cons 0 (+ 1 (send self :basis)))
542 (list (+ 1 (send self :basis))))))
543 (m (send self :sweep-matrix)))
544 (format t "~%REMOVEME2: Coef-ests: ~% Sweep Matrix: ~A ~% array dim 1: ~A ~% Swept indices: ~A ~% basis: ~A"
545 m n indices (send self :basis))
546 (coerce (compound-data-seq (select m (+ 1 n) indices)) 'list))) ;; ERROR
548 (defmeth regression-model-proto :xtxinv ()
549 "Message args: ()
550 Returns ((X^T) X)^(-1) or ((X^T) W X)^(-1)."
551 (let ((indices (if (send self :intercept)
552 (cons 0 (1+ (send self :basis)))
553 (1+ (send self :basis)))))
554 (select (send self :sweep-matrix) indices indices)))
556 (defmeth regression-model-proto :coef-standard-errors ()
557 "Message args: ()
558 Returns estimated standard errors of coefficients. Entries beyond the
559 intercept correspond to entries in basis."
560 (let ((s (send self :sigma-hat)))
561 (if s (* (send self :sigma-hat) (sqrt (diagonalf (send self :xtxinv)))))))
563 (defmeth regression-model-proto :studentized-residuals ()
564 "Message args: ()
565 Computes the internally studentized residuals for included cases and externally studentized residuals for excluded cases."
566 (let ((res (send self :residuals))
567 (lev (send self :leverages))
568 (sig (send self :sigma-hat))
569 (inc (send self :included)))
570 (if-else inc
571 (/ res (* sig (sqrt (pmax .00001 (- 1 lev)))))
572 (/ res (* sig (sqrt (+ 1 lev)))))))
574 (defmeth regression-model-proto :externally-studentized-residuals ()
575 "Message args: ()
576 Computes the externally studentized residuals."
577 (let* ((res (send self :studentized-residuals))
578 (df (send self :df)))
579 (if-else (send self :included)
580 (* res (sqrt (/ (- df 1) (- df (^ res 2)))))
581 res)))
583 (defmeth regression-model-proto :cooks-distances ()
584 "Message args: ()
585 Computes Cook's distances."
586 (let ((lev (send self :leverages))
587 (res (/ (^ (send self :studentized-residuals) 2)
588 (send self :num-coefs))))
589 (if-else (send self :included) (* res (/ lev (- 1 lev) )) (* res lev))))
592 (defun plot-points (x y &rest args)
593 "need to fix."
594 (declare (ignore x y args))
595 (error "Graphics not implemented yet."))
597 ;; Can not plot points yet!!
598 (defmeth regression-model-proto :plot-residuals (&optional x-values)
599 "Message args: (&optional x-values)
600 Opens a window with a plot of the residuals. If X-VALUES are not supplied
601 the fitted values are used. The plot can be linked to other plots with the
602 link-views function. Returns a plot object."
603 (plot-points (if x-values x-values (send self :fit-values))
604 (send self :residuals)
605 :title "Residual Plot"
606 :point-labels (send self :case-labels)))
608 (defmeth regression-model-proto :plot-bayes-residuals
609 (&optional x-values)
610 "Message args: (&optional x-values)
612 Opens a window with a plot of the standardized residuals and two
613 standard error bars for the posterior distribution of the actual
614 deviations from the line. See Chaloner and Brant. If X-VALUES are not
615 supplied the fitted values are used. The plot can be linked to other
616 plots with the link-views function. Returns a plot object."
618 (let* ((r (/ (send self :residuals)
619 (send self :sigma-hat)))
620 (d (* 2 (sqrt (send self :leverages))))
621 (low (- r d))
622 (high (+ r d))
623 (x-values (if x-values x-values (send self :fit-values)))
624 (p (plot-points x-values r
625 :title "Bayes Residual Plot"
626 :point-labels (send self :case-labels))))
627 (map 'list #'(lambda (a b c d) (send p :plotline a b c d nil))
628 x-values low x-values high)
629 (send p :adjust-to-data)