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