home *** CD-ROM | disk | FTP | other *** search
/ Garbo / Garbo.cdr / mac / science / xlspstr1.sit / regression.lsp < prev    next >
Lisp/Scheme  |  1990-08-16  |  18KB  |  443 lines
  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 "regression")
  14.  
  15. (require "statistics")
  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. (defmeth regression-model-proto :compute ()
  102. "Message args: ()
  103. Recomputes the estimates. For internal use by other messages"
  104.   (let* ((included (if-else (send self :included) 1 0))
  105.          (x (send self :x))
  106.          (y (send self :y))
  107.          (intercept (send self :intercept))
  108.          (weights (send self :weights))
  109.          (w (if weights (* included weights) included))
  110.          (m (make-sweep-matrix x y w))
  111.          (n (array-dimension x 1))
  112.          (p (- (array-dimension m 0) 1))
  113.          (tss (aref m p p))
  114.          (tol (* .0001 (mapcar #'standard-deviation (column-list x))))
  115.          (sweep-result
  116.           (if intercept
  117.               (sweep-operator m (iseq 1 n) tol)
  118.               (sweep-operator m (iseq 0 n) (cons 0.0 tol)))))
  119.     (setf (slot-value 'sweep-matrix) (first sweep-result))
  120.     (setf (slot-value 'total-sum-of-squares) tss)
  121.     (setf (slot-value 'residual-sum-of-squares) 
  122.           (aref (first sweep-result) p p))
  123.     (setf (slot-value 'basis)
  124.           (let ((b (remove 0 (second sweep-result))))
  125.             (if b 
  126.                 (- (reverse b) 1)
  127.                 (error "no columns could be swept"))))))
  128.  
  129. (defmeth regression-model-proto :needs-computing (&optional set)
  130.   (if set (setf (slot-value 'sweep-matrix) nil))
  131.   (null (slot-value 'sweep-matrix)))
  132.   
  133. (defmeth regression-model-proto :display ()
  134. "Message args: ()
  135. Prints the least squares regression summary. Variables not used in the fit
  136. are marked as aliased."
  137.   (let ((coefs (coerce (send self :coef-estimates) 'list))
  138.         (se-s (send self :coef-standard-errors))
  139.         (x (send self :x))
  140.         (p-names (send self :predictor-names)))
  141.     (if (send self :weights) 
  142.         (format t "~%Weighted Least Squares Estimates:~2%")
  143.         (format t "~%Least Squares Estimates:~2%"))
  144.     (when (send self :intercept)
  145.           (format t "Constant               ~10g   ~A~%"
  146.                   (car coefs) (list (car se-s)))
  147.           (setf coefs (cdr coefs))
  148.           (setf se-s (cdr se-s)))
  149.     (dotimes (i (array-dimension x 1)) 
  150.              (cond 
  151.                ((member i (send self :basis))
  152.                 (format t "~22a ~10g   ~A~%"
  153.                         (select p-names i) (car coefs) (list (car se-s)))
  154.                 (setf coefs (cdr coefs) se-s (cdr se-s)))
  155.                (t (format t "~22a    aliased~%" (select p-names i)))))
  156.     (format t "~%")
  157.     (format t "R Squared:             ~10g~%" (send self :r-squared))
  158.     (format t "Sigma hat:             ~10g~%" (send self :sigma-hat))
  159.     (format t "Number of cases:       ~10d~%" (send self :num-cases))
  160.     (if (/= (send self :num-cases) (send self :num-included))
  161.         (format t "Number of cases used:  ~10d~%" (send self :num-included)))
  162.     (format t "Degrees of freedom:    ~10d~%" (send self :df))
  163.     (format t "~%")))
  164.  
  165. ;;;
  166. ;;; Slot accessors and mutators
  167. ;;;
  168.  
  169. (defmeth regression-model-proto :x (&optional new-x)
  170. "Message args: (&optional new-x)
  171. With no argument returns the x matrix as supplied to m. With an argument
  172. NEW-X sets the x matrix to NEW-X and recomputes the estimates."
  173.   (when (and new-x (matrixp new-x))
  174.         (setf (slot-value 'x) new-x)
  175.         (send self :needs-computing t))
  176.   (slot-value 'x))
  177.  
  178. (defmeth regression-model-proto :y (&optional new-y)
  179. "Message args: (&optional new-y)
  180. With no argument returns the y sequence as supplied to m. With an argument
  181. NEW-Y sets the y sequence to NEW-Y and recomputes the estimates."
  182.   (when (and new-y (or (matrixp new-y) (sequencep new-y)))
  183.         (setf (slot-value 'y) new-y)
  184.         (send self :needs-computing t))
  185.   (slot-value 'y))
  186.  
  187. (defmeth regression-model-proto :intercept (&optional (val nil set))
  188. "Message args: (&optional new-intercept)
  189. With no argument returns T if the model includes an intercept term, nil if
  190. not. With an argument NEW-INTERCEPT the model is changed to include or
  191. exclude an intercept, according to the value of NEW-INTERCEPT."
  192.   (when set 
  193.         (setf (slot-value 'intercept) val)
  194.         (send self :needs-computing t))
  195.   (slot-value 'intercept))
  196.  
  197. (defmeth regression-model-proto :weights (&optional (new-w nil set))
  198. "Message args: (&optional new-w)
  199. With no argument returns the weight sequence as supplied to m; NIL means
  200. an unweighted model. NEW-W sets the weights sequence to NEW-W and
  201. recomputes the estimates."
  202.   (when set 
  203.         (setf (slot-value 'weights) new-w) 
  204.         (send self :needs-computing t))
  205.   (slot-value 'weights))
  206.  
  207. (defmeth regression-model-proto :total-sum-of-squares ()
  208. "Message args: ()
  209. Returns the total sum of squares around the mean."
  210.   (if (send self :needs-computing) (send self :compute))
  211.   (slot-value 'total-sum-of-squares))
  212.  
  213. (defmeth regression-model-proto :residual-sum-of-squares () 
  214. "Message args: ()
  215. Returns the residual sum of squares for the model."
  216.   (if (send self :needs-computing) (send self :compute))
  217.   (slot-value 'residual-sum-of-squares))
  218.  
  219. (defmeth regression-model-proto :basis ()
  220. "Message args: ()
  221. Returns the indices of the variables used in fitting the model."
  222.   (if (send self :needs-computing) (send self :compute))
  223.   (slot-value 'basis))
  224.  
  225. (defmeth regression-model-proto :sweep-matrix ()
  226. "Message args: ()
  227. Returns the swept sweep matrix. For internal use"
  228.   (if (send self :needs-computing) (send self :compute))
  229.   (slot-value 'sweep-matrix))
  230.  
  231. (defmeth regression-model-proto :included (&optional new-included)
  232. "Message args: (&optional new-included)
  233. 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."
  234.   (when (and new-included 
  235.              (= (length new-included) (send self :num-cases)))
  236.         (setf (slot-value 'included) (copy-seq new-included)) 
  237.         (send self :needs-computing t))
  238.   (if (slot-value 'included)
  239.       (slot-value 'included)
  240.       (repeat t (send self :num-cases))))
  241.  
  242. (defmeth regression-model-proto :predictor-names (&optional (names nil set))
  243. "Message args: (&optional (names nil set))
  244. With no argument returns the predictor names. NAMES sets the names."
  245.   (if set (setf (slot-value 'predictor-names) (mapcar #'string names)))
  246.   (let ((p (array-dimension (send self :x) 1))
  247.         (p-names (slot-value 'predictor-names)))
  248.     (if (not (and p-names (= (length p-names) p)))
  249.         (setf (slot-value 'predictor-names)
  250.               (mapcar #'(lambda (a) (format nil "Variable ~a" a)) 
  251.                       (iseq 0 (- p 1))))))
  252.   (slot-value 'predictor-names))
  253.  
  254. (defmeth regression-model-proto :response-name (&optional (name "Y" set))
  255. "Message args: (&optional name)
  256. With no argument returns the response name. NAME sets the name."
  257.   (if set (setf (slot-value 'response-name) (if name (string name) "Y")))
  258.   (slot-value 'response-name))
  259.  
  260. (defmeth regression-model-proto :case-labels (&optional (labels nil set))
  261. "Message args: (&optional labels)
  262. With no argument returns the case-labels. LABELS sets the labels."
  263.   (if set (setf (slot-value 'case-labels) 
  264.                 (if labels 
  265.                     (mapcar #'string labels)
  266.                     (mapcar #'(lambda (x) (format nil "~d" x)) 
  267.                             (iseq 0 (- (send self :num-cases) 1))))))
  268.   (slot-value 'case-labels))
  269.  
  270. ;;;
  271. ;;; Other Methods
  272. ;;; None of these methods access any slots directly.
  273. ;;;
  274.  
  275. (defmeth regression-model-proto :num-cases ()
  276. "Message args: ()
  277. Returns the number of cases in the model."
  278.   (length (send self :y)))
  279.  
  280. (defmeth regression-model-proto :num-included ()
  281. "Message args: ()
  282. Returns the number of cases used in the computations."
  283.   (sum (if-else (send self :included) 1 0)))
  284.  
  285. (defmeth regression-model-proto :num-coefs ()
  286. "Message args: ()
  287. Returns the number of coefficients in the fit model (including the
  288. intercept if the model includes one)."
  289.   (if (send self :intercept)
  290.       (+ 1 (length (send self :basis)))
  291.       (length (send self :basis))))
  292.  
  293. (defmeth regression-model-proto :df ()
  294. "Message args: ()
  295. Returns the number of degrees of freedom in the model."
  296.   (- (send self :num-included) (send self :num-coefs)))
  297.   
  298. (defmeth regression-model-proto :x-matrix ()
  299. "Message args: ()
  300. Returns the X matrix for the model, including a column of 1's, if
  301. appropriate. Columns of X matrix correspond to entries in basis."
  302.   (let ((m (select (send self :x) 
  303.                    (iseq 0 (- (send self :num-cases) 1)) 
  304.                    (send self :basis))))
  305.     (if (send self :intercept)
  306.         (bind-columns (repeat 1 (send self :num-cases)) m)
  307.         m)))
  308.  
  309. (defmeth regression-model-proto :leverages ()
  310. "Message args: ()
  311. Returns the diagonal elements of the hat matrix."
  312.   (let* ((weights (send self :weights))
  313.          (x (send self :x-matrix))
  314.          (raw-levs 
  315.           (matmult (* (matmult x (send self :xtxinv)) x)
  316.                    (repeat 1 (send self :num-coefs)))))
  317.     (if weights (* weights raw-levs) raw-levs)))
  318.  
  319. (defmeth regression-model-proto :fit-values ()
  320. "Message args: ()
  321. Returns the fitted values for the model."
  322.   (matmult (send self :x-matrix) (send self :coef-estimates)))
  323.  
  324. (defmeth regression-model-proto :raw-residuals () 
  325. "Message args: ()
  326. Returns the raw residuals for a model."
  327.   (- (send self :y) (send self :fit-values)))
  328.  
  329. (defmeth regression-model-proto :residuals () 
  330. "Message args: ()
  331. Returns the raw residuals for a model without weights. If the model
  332. includes weights the raw residuals times the square roots of the weights
  333. are returned."
  334.   (let ((raw-residuals (send self :raw-residuals))
  335.         (weights (send self :weights)))
  336.     (if weights (* (sqrt weights) raw-residuals) raw-residuals)))
  337.  
  338. (defmeth regression-model-proto :sum-of-squares () 
  339. "Message args: ()
  340. Returns the error sum of squares for the model."
  341.   (send self :residual-sum-of-squares))
  342.  
  343. (defmeth regression-model-proto :sigma-hat ()
  344. "Message args: ()
  345. Returns the estimated standard deviation of the deviations about the 
  346. regression line."
  347.   (let ((ss (send self :sum-of-squares))
  348.         (df (send self :df)))
  349.     (if (/= df 0) (sqrt (/ ss df)))))
  350.  
  351. ;; for models without an intercept the 'usual' formula for R^2 can give
  352. ;; negative results; hence the max.
  353. (defmeth regression-model-proto :r-squared ()
  354. "Message args: ()
  355. Returns the sample squared multiple correlation coefficient, R squared, for
  356. the regression."
  357.   (max (- 1 (/ (send self :sum-of-squares) (send self :total-sum-of-squares)))
  358.        0))
  359.  
  360. (defmeth regression-model-proto :coef-estimates ()
  361. "Message args: ()
  362. Returns the OLS (ordinary least squares) estimates of the regression
  363. coefficients. Entries beyond the intercept correspond to entries in basis."
  364.   (let ((n (array-dimension (send self :x) 1))
  365.         (indices (if (send self :intercept) 
  366.                      (cons 0 (+ 1 (send self :basis)))
  367.                      (+ 1 (send self :basis))))
  368.         (m (send self :sweep-matrix)))
  369.     (coerce (compound-data-seq (select m (+ 1 n) indices)) 'list)))
  370.  
  371. (defmeth regression-model-proto :xtxinv () 
  372. "Message args: ()
  373. Returns ((X^T) X)^(-1) or ((X^T) W X)^(-1)."
  374.   (let ((indices (if (send self :intercept) 
  375.                      (cons 0 (1+ (send self :basis))) 
  376.                      (1+ (send self :basis)))))
  377.     (select (send self :sweep-matrix) indices indices)))
  378.  
  379. (defmeth regression-model-proto :coef-standard-errors ()
  380. "Message args: ()
  381. Returns estimated standard errors of coefficients. Entries beyond the
  382. intercept correspond to entries in basis."
  383.   (let ((s (send self :sigma-hat)))
  384.     (if s (* (send self :sigma-hat) (sqrt (diagonal (send self :xtxinv)))))))
  385.  
  386. (defmeth regression-model-proto :studentized-residuals ()
  387. "Message args:  ()
  388. Computes the internally studentized residuals for included cases and externally studentized residuals for excluded cases."
  389.   (let ((res (send self :residuals))
  390.         (lev (send self :leverages))
  391.         (sig (send self :sigma-hat))
  392.         (inc (send self :included)))
  393.     (if-else inc
  394.              (/ res (* sig (sqrt (pmax .00001 (- 1 lev)))))
  395.              (/ res (* sig (sqrt (+ 1 lev)))))))
  396.  
  397. (defmeth regression-model-proto :externally-studentized-residuals ()
  398. "Message args:  ()
  399. Computes the externally studentized residuals."
  400.   (let* ((res (send self :studentized-residuals))
  401.          (df (send self :df)))
  402.     (if-else (send self :included)
  403.              (* res (sqrt (/ (- df 1) (- df (^ res 2)))))
  404.              res)))
  405.  
  406. (defmeth regression-model-proto :cooks-distances ()
  407. "Message args: ()
  408. Computes Cook's distances."
  409.   (let ((lev (send self :leverages))
  410.         (res (/ (^ (send self :studentized-residuals) 2)
  411.                 (send self :num-coefs))))
  412.     (if-else (send self :included) (* res (/ lev (- 1 lev) )) (* res lev))))
  413.  
  414. (defmeth regression-model-proto :plot-residuals (&optional x-values)
  415. "Message args: (&optional x-values)
  416. Opens a window with a plot of the residuals. If X-VALUES are not supplied 
  417. the fitted values are used. The plot can be linked to other plots with the 
  418. link-views function. Returns a plot object."
  419.   (plot-points (if x-values x-values (send self :fit-values))
  420.                (send self :residuals)
  421.                :title "Residual Plot"
  422.                :point-labels (send self :case-labels)))
  423.  
  424. (defmeth regression-model-proto :plot-bayes-residuals 
  425.   (&optional x-values)
  426. "Message args: (&optional x-values)
  427. Opens a window with a plot of the standardized residuals and two standard
  428. error bars for the posterior distribution of the actual deviations from the
  429. line. See Chaloner and Brant. If X-VALUES are not supplied  the fitted values
  430. are used. The plot can be linked to other plots with the link-views function.
  431. Returns a plot object."
  432.   (let* ((r (/ (send self :residuals) (send self :sigma-hat)))
  433.          (d (* 2 (sqrt (send self :leverages))))
  434.          (low (- r d))
  435.          (high (+ r d))
  436.          (x-values (if x-values x-values (send self :fit-values)))
  437.          (p (plot-points x-values r :title "Bayes Residual Plot"
  438.                          :point-labels (send self :case-labels))))
  439.     (map 'list #'(lambda (a b c d) (send p :plotline a b c d nil))
  440.                x-values low x-values high)
  441.     (send p :adjust-to-data)
  442.     p))
  443.