home *** CD-ROM | disk | FTP | other *** search
/ Garbo / Garbo.cdr / mac / science / xlspstrd.sit / glim / glim.lsp < prev    next >
Text File  |  1991-02-19  |  24KB  |  796 lines
  1. ;;;;
  2. ;;;;
  3. ;;;;                  A Simple GLIM Implementation
  4. ;;;;
  5. ;;;;
  6.  
  7. (provide "glim")
  8.  
  9. ;;;;
  10. ;;;;                         Link Prototypes
  11. ;;;;
  12. ;;;;
  13. ;;;; Links are objects responding to three messages:
  14. ;;;;
  15. ;;;;     :eta     takes a set of mean values and returns the linear
  16. ;;;;              predictor values
  17. ;;;;     :means   takes a set of linear predictor values and returns
  18. ;;;;              the mean values (the inverse of :eta)
  19. ;;;;     :derivs  takes a set of mean values and returns the values of
  20. ;;;;              the derivatives of the linear predictors at the mean
  21. ;;;;              values
  22. ;;;;
  23. ;;;; The arguments should be sequences. The glim-link-proto prototype
  24. ;;;; implements an identity link. Links for binomial errors are defined
  25. ;;;; for means in the unit interval [0, 1], i. e. for n = 1 trials.
  26. ;;;;
  27.  
  28. (defproto glim-link-proto)
  29.  
  30. (defmeth glim-link-proto :eta (mu)
  31. "Method args: (mu)
  32. Returns linear predictor values at MU."
  33.   mu)
  34.  
  35. (defmeth glim-link-proto :means (eta)
  36. "Method args: (eta)
  37. Returns mean values for linear predictor ETA."
  38.   eta)
  39.  
  40. (defmeth glim-link-proto :derivs (mu)
  41. "Method args: (mu)
  42. Returns d(eta)/d(mu) values at MU."
  43.   (repeat 1 (length mu)))
  44.  
  45. (defmeth glim-link-proto :print (&optional (stream t))
  46.   (format stream "#<Glim Link Object: ~s>" (slot-value 'proto-name)))
  47.  
  48. (defmeth glim-link-proto :save ()
  49.   (let ((proto (slot-value 'proto-name)))
  50.     (if (eq self (eval proto)) proto `(send ,proto :new))))
  51.  
  52. ;;;;
  53. ;;;; Identity Link Prototype
  54. ;;;;
  55.  
  56. (defproto identity-link () () glim-link-proto)
  57.  
  58. ;;;;
  59. ;;;; Log Link Prototype
  60. ;;;;
  61.  
  62. (defproto log-link () () glim-link-proto)
  63.  
  64. (defmeth log-link :eta (mu) (log mu))
  65. (defmeth log-link :means (eta) (exp eta))
  66. (defmeth log-link :derivs (mu) (/ mu))
  67.  
  68. ;;;;
  69. ;;;; Inverse-link-prototype
  70. ;;;;
  71.  
  72. (defproto inverse-link () () glim-link-proto)
  73.  
  74. (defmeth inverse-link :eta (mu) (/ mu))
  75. (defmeth inverse-link :means (eta) (/ eta))
  76. (defmeth inverse-link :derivs (mu) (- (/ (^ mu 2))))
  77.  
  78. ;;;;
  79. ;;;; Square Root Link
  80. ;;;;
  81.  
  82. (defproto sqrt-link () () glim-link-proto)
  83.  
  84. (defmeth sqrt-link :eta (mu) (sqrt mu))
  85. (defmeth sqrt-link :means (eta) (^ eta 2))
  86. (defmeth sqrt-link :derivs (mu) (/ 0.5 (sqrt mu)))
  87.  
  88. ;;;;
  89. ;;;; Power Link Prototype
  90. ;;;;
  91.  
  92. (defproto power-link-proto '(power) () glim-link-proto)
  93.  
  94. (defmeth power-link-proto :isnew (power) (setf (slot-value 'power) power))
  95. (defmeth power-link-proto :power () (slot-value 'power))
  96.  
  97. (defmeth power-link-proto :print (&optional (stream t))
  98.   (format stream "#<Glim Link Object: Power Link (~s)>" (send self :power)))
  99.  
  100. (defmeth power-link-proto :save ()
  101.   `(send power-link-proto :new ,(send self :power)))
  102.  
  103. (defmeth power-link-proto :eta (mu) (^ mu (send self :power)))
  104. (defmeth power-link-proto :means (eta) (^ eta (/ (slot-value 'power))))
  105. (defmeth power-link-proto :derivs (mu)
  106.   (let ((p (slot-value 'power)))
  107.     (* p (^ mu (- p 1)))))
  108.  
  109. ;;;;
  110. ;;;; Logit Link Prototype
  111. ;;;;
  112.  
  113. (defproto logit-link () () glim-link-proto)
  114.  
  115. (defmeth logit-link :eta (p) (log (/ p (- 1 p))))
  116. (defmeth logit-link :means (eta)
  117.   (let ((exp-eta (exp eta)))
  118.     (/ exp-eta (+ 1 exp-eta))))
  119. (defmeth logit-link :derivs (p) (+ (/ p) (/ (- 1 p))))
  120.  
  121. ;;;;
  122. ;;;; Probit Link Prototype
  123. ;;;;
  124.  
  125. (defproto probit-link () () glim-link-proto)
  126.  
  127. (defmeth probit-link :eta (p) (normal-quant p))
  128. (defmeth probit-link :means (eta) (normal-cdf eta))
  129. (defmeth probit-link :derivs (p) (/ 1 (normal-dens (normal-quant p))))
  130.  
  131. ;;;;
  132. ;;;; Complimentary Log-Log Link Prototype
  133. ;;;;
  134.  
  135. (defproto cloglog-link () () glim-link-proto)
  136.  
  137. (defmeth cloglog-link :eta (p) (log (- (log (- 1 p)))))
  138. (defmeth cloglog-link :means (eta) (- 1 (exp (- (exp eta)))))
  139. (defmeth cloglog-link :derivs (p)
  140.   (let ((q (- 1 p)))
  141.     (/ -1 (log q) q)))
  142.  
  143. ;;;;
  144. ;;;;
  145. ;;;;                The General GLIM Prototype
  146. ;;;;         (Uses Normal Errors and an Identity Link)
  147. ;;;;
  148. ;;;;
  149.  
  150. (defproto glim-proto 
  151.   '(yvar link offset pweights scale est-scale
  152.          epsilon epsilon-dev count-limit verbose recycle
  153.          eta deviances)
  154.   '() 
  155.   regression-model-proto)
  156.  
  157. ;;;;
  158. ;;;; Slot Accessors
  159. ;;;;
  160.  
  161. (defmeth glim-proto :yvar (&optional (new nil set))
  162. "Message args: (&optional new)
  163. Sets or returns dependent variable."
  164.   (when set
  165.         (setf (slot-value 'yvar) new)
  166.         (send self :needs-computing t))
  167.   (slot-value 'yvar))
  168.  
  169. (defmeth glim-proto :link (&optional (new nil set))
  170. "Message args: (&optional new)
  171. Sets or returns link object."
  172.   (when set
  173.         (setf (slot-value 'link) new)
  174.         (send self :needs-computing t))
  175.   (slot-value 'link))
  176.  
  177. (defmeth glim-proto :offset (&optional (new nil set))
  178. "Message args: (&optional (new nil set))
  179. Sets or returns offset values."
  180.   (when set
  181.         (setf (slot-value 'offset) new)
  182.         (send self :needs-computing t))
  183.   (slot-value 'offset))
  184.  
  185. (defmeth glim-proto :pweights (&optional (new nil set))
  186. "Message args: (&optional (new nil set))
  187. Sets or returns prior weights."
  188.   (when set
  189.         (setf (slot-value 'pweights) new)
  190.         (send self :needs-computing t))
  191.   (slot-value 'pweights))
  192.  
  193. ;; changing the scale does not require recomputing the estimates
  194. (defmeth glim-proto :scale (&optional (new nil set))
  195. "Message args: (&optional (new nil set))
  196. Sets or returns value of scale parameter."
  197.   (if set (setf (slot-value 'scale) new))
  198.   (slot-value 'scale))
  199.  
  200. (defmeth glim-proto :estimate-scale (&optional (val nil set))
  201. "Message args: (&optional (val nil set))
  202. Sets or returns value of ESTIMATE-SCALE option."
  203.   (if set (setf (slot-value 'est-scale) val))
  204.   (slot-value 'est-scale))
  205.  
  206. (defmeth glim-proto :epsilon (&optional new)
  207. "Message args: (&optional new)
  208. Sets or returns tolerance for relative change in coefficients."
  209.   (if new (setf (slot-value 'epsilon) new))
  210.   (slot-value 'epsilon))
  211.  
  212. (defmeth glim-proto :epsilon-dev (&optional new)
  213. "Message args: (&optional new)
  214. Sets or returns tolerance for change in deviance."
  215.   (if new (setf (slot-value 'epsilon-dev) new))
  216.   (slot-value 'epsilon-dev))
  217.  
  218. (defmeth glim-proto :count-limit (&optional new)
  219. "Message args: (&optional new)
  220. Sets or returns maximum number of itrations."
  221.   (if new (setf (slot-value 'count-limit) new))
  222.   (slot-value 'count-limit))
  223.  
  224. (defmeth glim-proto :recycle (&optional (new nil set))
  225. "Message args: (&optional new)
  226. Sets or returns recycle option. If option is not NIL, current values
  227. are used as initial values by :COMPUTE method."
  228.   (when set
  229.         (setf (slot-value 'recycle) new))
  230.   (slot-value 'recycle))
  231.  
  232. (defmeth glim-proto :verbose (&optional (val nil set))
  233. "Message args: (&optional (val nil set))
  234. Sets or returns VERBOSE option. Iteration info is printed if option
  235. is not NIL."
  236.   (if set (setf (slot-value 'verbose) val))
  237.   (slot-value 'verbose))
  238.  
  239. (defmeth glim-proto :eta ()
  240. "Message args: ()
  241. Returns linear predictor values for durrent fit."
  242.   (slot-value 'eta))
  243.  
  244. (defmeth glim-proto :set-eta (&optional val)
  245.   (if val
  246.       (setf (slot-value 'eta) val)
  247.       (setf (slot-value 'eta)
  248.             (+ (send self :offset) (send self :fit-values)))))
  249.  
  250. (defmeth glim-proto :deviances ()
  251. "Message args: ()
  252. Returns deviances for durrent fit."
  253.   (slot-value 'deviances))
  254.  
  255. (defmeth glim-proto :set-deviances ()
  256.   (setf (slot-value 'deviances) 
  257.         (send self :fit-deviances (send self :fit-means))))
  258.  
  259. ;;;;
  260. ;;;; Overrides for Regression Methods
  261. ;;;;
  262.  
  263. ;; A variant of this method should work for any object whose slot values
  264. ;; have valid printed representations.
  265. (defmeth glim-proto :save ()
  266.   (let* ((proto (slot-value 'proto-name))
  267.      (slots (remove 'link (send self :own-slots)))
  268.      (values (mapcar #'slot-value slots)))
  269.     `(let ((object (make-object ,proto))
  270.        (slots ',slots)
  271.        (values ',values))
  272.        (flet ((add-slot (s v) (send object :add-slot s v)))
  273.      (mapcar #'add-slot slots values)
  274.      (add-slot 'link ,(send (send self :link) :save)))
  275.        object)))
  276.  
  277. (defmeth glim-proto :sigma-hat () (sqrt (send self :scale)))
  278.  
  279. ;; This override is only used to modify the documentation string.
  280. (defmeth glim-proto :fit-values ()
  281. "Message args: ()
  282. Returns Xb, the linear predictor values without the offset.
  283. The :fit-means method returns fitted means for the current estimates."
  284.   (call-next-method))
  285.  
  286. ;; this should be merged with the regression-model method
  287. (defmeth glim-proto :x (&optional x)
  288.   (if x
  289.       (let ((x (cond
  290.                 ((matrixp x) x)
  291.                 ((vectorp x) (list x))
  292.                 ((and (consp x) (numberp (car x))) (list x))
  293.                 (t x))))
  294.         (call-next-method (if (matrixp x) x (apply #'bind-columns x)))))
  295.   (call-next-method))
  296.  
  297. (defmeth glim-proto :raw-residuals () 
  298. "Message args: ()
  299. Returns the raw residuals for a model."
  300.   (- (send self :yvar) (send self :fit-means)))
  301.  
  302. ;; This override is needed because regression-model-proto defines its
  303. ;; residuals in terms of :raw-residuals.
  304. (defmeth regression-model-proto :residuals ()
  305. "Message args: ()
  306. Returns the Pearson residuals."
  307.   (let ((raw-residuals (- (send self :y) (send self :fit-values)))
  308.         (weights (send self :weights)))
  309.     (if weights (* (sqrt weights) raw-residuals) raw-residuals)))
  310.  
  311. ;;;;
  312. ;;;; Computing methods
  313. ;;;;
  314.  
  315. (defmeth glim-proto :compute ()
  316.   (let* ((epsilon (send self :epsilon))
  317.          (epsilon-dev (send self :epsilon-dev))
  318.          (maxcount (send self :count-limit))
  319.          (low-lim (* 2 (/ machine-epsilon epsilon)))
  320.          (verbose (send self :verbose)))
  321.     (unless (and (send self :eta) (send self :recycle))
  322.             (send self :initialize-search))
  323.     (send self :compute-step)
  324.     (do ((count 1 (+ count 1))
  325.          (beta 0 (send self :coef-estimates))
  326.          (last-beta -1 beta)
  327.          (dev  0 (send self :deviance))
  328.          (last-dev  -1 dev))
  329.         ((or (> count maxcount) 
  330.              (< (max (abs (/ (- beta last-beta)
  331.                              (pmax (abs last-beta) low-lim))))
  332.                 epsilon)
  333.              (< (abs (- dev last-dev)) epsilon-dev)))
  334.         (if verbose 
  335.             (format t "Iteration ~d: deviance = ~a~%" 
  336.                     count (send self :deviance)))
  337.         (send self :compute-step))))
  338.  
  339. (defmeth glim-proto :compute-step ()
  340. "Args: ()
  341. Executes one iteratively reweighted least squares step."
  342.   (let* ((yvar (send self :yvar))
  343.          (offset (send self :offset))
  344.          (eta (send self :eta))
  345.          (mu (send self :fit-means eta))
  346.          (d-eta (send self :fit-link-derivs mu))
  347.          (z (- (+ eta (* (- yvar mu) d-eta)) offset))
  348.          (v (send self :fit-variances mu))
  349.          (w-inv (* d-eta d-eta v))
  350.          (pw (send self :pweights)))
  351.     (send self :y z)
  352.     (send self :weights (if pw (/ pw w-inv) (/ w-inv)))
  353.     (call-method regression-model-proto :compute)
  354.     (send self :set-eta)
  355.     (send self :set-deviances)
  356.     (if (send self :estimate-scale) 
  357.         (send self :scale (send self :fit-scale)))))
  358.  
  359. (defmeth glim-proto :deviance () 
  360. "Message args: ()
  361. Returns deviance for included cases."
  362.   (sum (if-else (send self :included) (send self :deviances) 0)))
  363.  
  364. (defmeth glim-proto :mean-deviance ()
  365. "Message args: ()
  366. Returns mean deviance for included cases, adjusted for degrees of
  367. freedom."
  368.   (/ (send self :deviance) (send self :df)))
  369.  
  370. (defmeth glim-proto :initialize-search (&optional eta)
  371.   (send self :set-eta
  372.     (if eta eta (send (send self :link) :eta (send self :initial-means))))
  373.   (send self :needs-computing t))
  374.  
  375. (defmeth glim-proto :fit-means (&optional (eta (send self :eta)))
  376. "Message args: (&optional (eta (send self :eta)))
  377. Retruns mean values for current or supplied ETA."
  378.   (send (send self :link) :means eta))
  379.  
  380. (defmeth glim-proto :fit-link-derivs (mu)
  381. "Message args: ()
  382. Returns link derivative values at MU."
  383.   (send (send self :link) :derivs mu))
  384.  
  385. (defmeth glim-proto :display ()
  386. "Message args: ()
  387. Prints the IRWLS regression summary. Variables not used in the fit are
  388. marked as aliased."
  389.   (let ((coefs (coerce (send self :coef-estimates) 'list))
  390.         (se-s (send self :coef-standard-errors))
  391.         (x (send self :x))
  392.         (p-names (send self :predictor-names)))
  393.     (if (send self :weights) 
  394.         (format t "~%Weighted Least Squares Estimates:~2%")
  395.         (format t "~%Least Squares Estimates:~2%"))
  396.     (when (send self :intercept)
  397.           (format t "Constant               ~10g   ~A~%"
  398.                   (car coefs) (list (car se-s)))
  399.           (setf coefs (cdr coefs))
  400.           (setf se-s (cdr se-s)))
  401.     (dotimes (i (array-dimension x 1)) 
  402.              (cond 
  403.                ((member i (send self :basis))
  404.                 (format t "~22a ~10g   ~A~%"
  405.                         (select p-names i) (car coefs) (list (car se-s)))
  406.                 (setf coefs (cdr coefs) se-s (cdr se-s)))
  407.                (t (format t "~22a    aliased~%" (select p-names i)))))
  408.     (format t "~%")
  409.     (if (send self :estimate-scale)
  410.         (format t "Scale Estimate:        ~10g~%" (send self :scale))
  411.         (format t "Scale taken as:        ~10g~%" (send self :scale)))
  412.     (format t "Deviance:              ~10g~%" (send self :deviance))
  413.     (format t "Number of cases:       ~10d~%" (send self :num-cases))
  414.     (if (/= (send self :num-cases) (send self :num-included))
  415.         (format t "Number of cases used:  ~10d~%" (send self :num-included)))
  416.     (format t "Degrees of freedom:    ~10d~%" (send self :df))
  417.     (format t "~%")))
  418.  
  419. ;;;;
  420. ;;;; Error-Dependent Methods (Normal Errors)
  421. ;;;;
  422.  
  423. (defmeth glim-proto :initial-means ()
  424. "Message args: ()
  425. Returns initial means estimate for starting the iteration."
  426.   (send self :yvar))
  427.  
  428. (defmeth glim-proto :fit-variances (mu)
  429. "Message args: (mu)
  430. Returns variance function values at MU."
  431.   (repeat 1 (length mu)))
  432.  
  433. (defmeth glim-proto :fit-deviances (mu)
  434. "Message args: (mu)
  435. Returns deviance values at MU."
  436.   (let ((raw-dev (^ (- (send self :yvar) mu) 2))
  437.                   (pw (send self :pweights)))
  438.        (if pw (* pw raw-dev) raw-dev)))
  439.  
  440. (defmeth glim-proto :fit-scale ()
  441. "Message args: ()
  442. Returns estimate of scale parameter."
  443.   (send self :mean-deviance))
  444.  
  445. ;;;;
  446. ;;;; Initial values for the prototype
  447. ;;;;
  448.  
  449. (send glim-proto :scale 1.0)
  450. (send glim-proto :offset 0.0)
  451. (send glim-proto :link identity-link)
  452. (send glim-proto :estimate-scale t)
  453. (send glim-proto :epsilon .000001)
  454. (send glim-proto :epsilon-dev .001)
  455. (send glim-proto :count-limit 30)
  456. (send glim-proto :verbose t)
  457.  
  458. ;;;;
  459. ;;;; :ISNEW method
  460. ;;;;
  461.  
  462. (defmeth glim-proto :isnew (&key x 
  463.                                  y
  464.                  link
  465.                                  (offset 0)
  466.                                  (intercept t)
  467.                                  included
  468.                                  pweights
  469.                                  (print (and x y))
  470.                                  (verbose t)
  471.                                  predictor-names
  472.                                  response-name
  473.                                  (recycle nil)
  474.                                  case-labels)
  475.   (send self :x x)
  476.   (send self :y y)
  477.   (send self :yvar y)
  478.   (if link (send self :link link))
  479.   (send self :offset offset)
  480.   (send self :intercept intercept)
  481.   (send self :pweights pweights)
  482.   (send self :recycle recycle)
  483.   (send self :verbose verbose)
  484.   (if included (send self :included included))
  485.   (if predictor-names (send self :predictor-names predictor-names))
  486.   (if response-name (send self :response-name response-name))
  487.   (if (or y case-labels) (send self :case-labels case-labels)) ; needs fixing
  488.   (if print (send self :display)))
  489.  
  490. ;;;;
  491. ;;;; Some Additional Residual Methods
  492. ;;;;
  493.  
  494. (defmeth glim-proto :chi-residuals ()
  495. "Message args: ()
  496. Returns the components of Pearson's chi-squared residuals."
  497.   (send self :residuals))
  498.  
  499. (defmeth glim-proto :standardized-chi-residuals ()
  500. "Message args: ()
  501. Returns the components of Standardized Pearson Residuals (Williams, 1987)."
  502.   (send self :studentized-residuals))
  503.  
  504. (defmeth glim-proto :deviance-residuals ()
  505. "Message args: ()
  506. Returns the components of deviance residuals for non binomial models."
  507.   (let* ((dev (sqrt (send self :deviances)))
  508.          (sign (if-else (< (send self :yvar) (send self :fit-means)) -1 1)))
  509.     (* sign dev)))
  510.  
  511. (defmeth glim-proto :standardized-deviance-residuals ()
  512. "Message args: ()
  513. Returns the standardized deviance residuals, (Davison and Tsai, 1989)."
  514.   (let* ((dev (send self :deviance-residuals))
  515.          (inc (send self :included))
  516.          (h (send self :leverages)))
  517.     (if-else inc
  518.              (/ dev (sqrt (* (send self :scale) (- 1 h))))
  519.              (/ dev (sqrt (* (send self :scale) (+ 1 h)))))))
  520.  
  521. (defmeth glim-proto :g2-residuals ()
  522. "Message args: ()
  523. Returns  a weighted combination of the standardized deviance and chi
  524. residuals, (Davison and Tsai, 1989)."
  525.   (let* ((dev (send self :standardized-deviance-residuals))
  526.          (chi (send self :standardized-chi-residuals))
  527.          (inc (send self :included))
  528.          (h (send self :leverages))
  529.          (sign (if-else (< dev 0) -1 1)))
  530.     (* sign (sqrt (+ (* (1- h) (^ dev 2))
  531.                      (*   h    (^ chi 2)))))))
  532.  
  533. ;;;;
  534. ;;;;
  535. ;;;;                 Normal Regression Model Prototype
  536. ;;;;
  537. ;;;;
  538.  
  539. (defproto normalreg-proto () () glim-proto)
  540.  
  541. ;;;;
  542. ;;;; Normal Model Constructor Function
  543. ;;;;
  544.  
  545. (defun normalreg-model (x y &rest args)
  546. "Args: (x y &rest args)
  547. Returns a normal regression model. Accepts :LINK, :OFFSET and :VERBOSE
  548. keywords in addition to the keywords accepted by regression-model."
  549.   (apply #'send normalreg-proto :new :x x :y y args))
  550.  
  551. ;;;;
  552. ;;;;
  553. ;;;;             Poisson Regression Model Prototype
  554. ;;;;
  555. ;;;;
  556.  
  557. (defproto poissonreg-proto () () glim-proto)
  558.  
  559. ;;;;
  560. ;;;; Error-Dependent Methods (Poisson Errors)
  561. ;;;;
  562.  
  563. (defmeth poissonreg-proto :initial-means () (pmax (send self :yvar) 0.5))
  564.  
  565. (defmeth poissonreg-proto :fit-variances (mu) mu)
  566.  
  567. (defmeth poissonreg-proto :fit-deviances (mu)
  568.   (flet ((log+ (x) (log (if-else (< 0 x) x 1)))) ; to prevent log of zero
  569.     (let* ((y (send self :yvar))
  570.            (raw-dev (* 2 (- (* y (log+ (/ y mu))) (- y mu))))
  571.            (pw (send self :pweights)))
  572.       (if pw (* pw raw-dev) raw-dev))))
  573.  
  574. ;;;;
  575. ;;;; Initial values for the prototype
  576. ;;;;
  577.  
  578. (send poissonreg-proto :estimate-scale nil)
  579. (send poissonreg-proto :link log-link)
  580.  
  581. ;;;;
  582. ;;;; Poisson Model Constructor Functions
  583. ;;;;
  584.  
  585. (defun poissonreg-model (x y &rest args)
  586. "Args: (x y &rest args)
  587. Returns a Poisson regression model. Accepts :LINK, :OFFSET and :VERBOSE
  588. keywords in addition to the keywords accepted by regression-model."
  589.   (apply #'send poissonreg-proto :new :x x :y y args))
  590.  
  591. (defun loglinreg-model (x y &rest args)
  592. "Args: (x y &rest args)
  593. Returns a Poisson regression model with a log link. Accepts :OFFSET and
  594. :VERBOSE keywords in addition to the keywords accepted by regression-model."
  595.   (apply #'send poissonreg-proto :new :x x :y y :link log-link args))
  596.  
  597. ;;;;
  598. ;;;;
  599. ;;;;             Binomial Regression Model Prototype
  600. ;;;;
  601. ;;;;
  602.  
  603. (defproto binomialreg-proto '(trials) () glim-proto)
  604.  
  605. ;;;;
  606. ;;;; Slot Accessor
  607. ;;;;
  608.  
  609. (defmeth binomialreg-proto :trials (&optional new)
  610. "Message args: ()
  611. Sets or retruns number of trials for each observation."
  612.   (when new
  613.         (setf (slot-value 'trials) new)
  614.         (send self :needs-computing t))
  615.   (slot-value 'trials))
  616.  
  617. ;;;;
  618. ;;;; Overrides for link-related methods to incorporate trials
  619. ;;;;
  620.  
  621. (defmeth binomialreg-proto :fit-means (&optional (eta (send self :eta)))
  622.   (let ((n (send self :trials))
  623.     (p (call-next-method eta)))
  624.     (* n p)))
  625.  
  626. (defmeth  binomialreg-proto :fit-link-derivs (mu)
  627.   (let* ((n (send self :trials))
  628.      (d (call-next-method (/ mu n))))
  629.     (/ d n)))
  630.  
  631. (defmeth binomialreg-proto :initialize-search (&optional eta)
  632.   (call-next-method 
  633.    (if eta eta (send (send self :link) :eta (send self :initial-probs)))))
  634.  
  635. ;;;;
  636. ;;;; Error-Dependent Methods (Binomial Errors)
  637. ;;;;
  638.  
  639. (defmeth binomialreg-proto :initial-probs ()
  640.   (let* ((n (send self :trials))
  641.      (p (/ (pmax (pmin (send self :yvar) (- n 0.5)) 0.5) n)))
  642.     p))
  643.  
  644. (defmeth binomialreg-proto :initial-means ()
  645.   (* (send self :trials) (send self :initial-probs)))
  646.  
  647. (defmeth binomialreg-proto :fit-variances (mu)
  648.   (let* ((n (send self :trials))
  649.          (p (/ mu n)))
  650.     (* n p (- 1 p))))
  651.  
  652. (defmeth binomialreg-proto :fit-deviances (mu)
  653.   (flet ((log+ (x) (log (if-else (< 0 x) x 1)))) ; to prevent log of zero
  654.     (let* ((n (send self :trials))
  655.            (y (send self :yvar))
  656.            (n-y (- n y))
  657.            (n-mu (- n mu))
  658.            (pw (send self :pweights))
  659.            (raw-dev (* 2 (+ (* y (log+ (/ y mu))) 
  660.                             (* n-y (log+ (/ n-y n-mu)))))))
  661.       (if pw (* pw raw-dev) raw-dev))))
  662.  
  663. ;;;;
  664. ;;;; Other Methods
  665. ;;;;
  666.  
  667. (defmeth binomialreg-proto :fit-probabilities ()
  668. "Message args: ()
  669. Returns the fitted probabilities for the model."
  670.   (/ (send self :fit-means) (send self :trials)))
  671.  
  672. ;;;;
  673. ;;;; :ISNEW method
  674. ;;;;
  675.  
  676. (defmeth binomialreg-proto :isnew (&rest args &key trials)
  677.   (send self :trials trials)
  678.   (apply #'call-next-method args))
  679.  
  680. ;;;;
  681. ;;;; Initial values for the prototype
  682. ;;;;
  683.  
  684. (send binomialreg-proto :estimate-scale nil)
  685. (send binomialreg-proto :link logit-link)
  686.  
  687. ;;;;
  688. ;;;; Binomial Model Constructor Functions
  689. ;;;;
  690.  
  691. (defun binomialreg-model (x y n &rest args)
  692. "Args: (x y n &rest args)
  693. Returns a binomial regression model. Accepts :LINK, :OFFSET and :VERBOSE
  694. keywords in addition to the keywords accepted by regression-model."
  695.   (apply #'send binomialreg-proto :new :x x :y y :trials n args))
  696.  
  697. (defun logitreg-model (x y n &rest args)
  698. "Args: (x y n &rest args)
  699. Returns a logistic regression model (binomial regression model with logit
  700. link). Accepts :OFFSET and :VERBOSE keywords in addition to the keywords
  701. accepted by regression-model."
  702.   (apply #'send binomialreg-proto :new 
  703.      :x x :y y :trials n :link logit-link args))
  704.  
  705. (defun probitreg-model (x y n &rest args)
  706. "Args: (x y n &rest args)
  707. Returns a probit regression model (binomial regression model with probit
  708. link). Accepts :OFFSET and :VERBOSE keywords in addition to the keywords
  709. accepted by regression-model."
  710.   (apply #'send binomialreg-proto :new
  711.      :x x :y y :trials n :link probit-link args))
  712.  
  713. ;;;;
  714. ;;;;
  715. ;;;;               Gamma Regression Model Prototype
  716. ;;;;
  717. ;;;;
  718.  
  719. (defproto gammareg-proto () () glim-proto)
  720.  
  721. ;;;;
  722. ;;;; Error-Dependent Methods
  723. ;;;;
  724.  
  725. (defmeth gammareg-proto :initial-means () (pmax (send self :yvar) 0.5))
  726.  
  727. (defmeth gammareg-proto :fit-variances (mu) (^ mu 2))
  728.  
  729. (defmeth gammareg-proto :fit-deviances (mu)
  730.   (let* ((y (send self :yvar))
  731.      (pw (send self :pweights))
  732.          (raw-dev (* 2 (+ (- (log (/ y mu))) (/ (- y mu) mu)))))
  733.     (if pw (* raw-dev pw) raw-dev)))
  734.  
  735. ;;;;
  736. ;;;; Initial values for the prototype
  737. ;;;;
  738.  
  739. (send gammareg-proto :link inverse-link)
  740.  
  741. ;;;;
  742. ;;;; Gamma Model Constructor Function
  743. ;;;;
  744.  
  745. (defun gammareg-model (x y &rest args)
  746. "Args: (x y &rest args)
  747. Returns a Gamma regression model. Accepts :LINK, :OFFSET and :VERBOSE
  748. keywords in addition to the keywords accepted by regression-model."
  749.   (apply #'send gammareg-proto :new :x x :y y args))
  750.  
  751. ;;;;
  752. ;;;;
  753. ;;;;                Some Simple Design Matrix Tools
  754. ;;;;
  755. ;;;;
  756.  
  757. (defun indicators (x &key (drop-first t) (test #'eql))
  758. "Args: (x &key (drop-first t) (test #'eql))
  759. Returns a list of indicators sequences for the levels of X. TEST is
  760. used to check equality of levels. If DROP-FIRST is true, the indicator
  761. for the first level is dropped."
  762.   (let ((levels (remove-duplicates (coerce x 'list))))
  763.     (mapcar #'(lambda (lev) (if-else (map-elements test lev x) 1 0))
  764.             (if drop-first (rest levels) levels))))
  765.  
  766. (defun cross-terms (x &rest args)
  767. "Args: (x &rest args)
  768. Arguments should be lists. Returns list of cross products, with the first
  769. argument list varying slowest."
  770.   (case (length args)
  771.     (0 (error "too few arguments"))
  772.     (1 (let ((y (first args)))
  773.          (apply #'append 
  774.                 (mapcar #'(lambda (a) (mapcar #'(lambda (b) (* a b)) y)) x))))
  775.     (t (cross-terms x (apply #'cross-terms args)))))
  776.  
  777. (defun level-names (x &key (prefix "") (drop-first t))
  778. "Args: (x &key (prefix "") (drop-first t))
  779. Constructs name strings using unique levels in X and PREFIX."
  780.   (let ((levels (remove-duplicates (coerce x 'list))))
  781.     (mapcar #'(lambda (x) (format nil "~a(~a)" prefix x))
  782.             (if drop-first (rest levels) levels))))
  783.  
  784. (defun cross-names (x &rest args)
  785. "Args: (x &rest args)
  786. Arguments should be lists. Constructs cross products of names, separated
  787. by dots. First index varies slowest."
  788.   (flet ((paste (x y) (format nil "~a.~a" x y)))
  789.     (case (length args)
  790.       (0 (error "too few arguments"))
  791.       (1 (let ((y (first args)))
  792.            (apply #'append
  793.                   (mapcar #'(lambda (a) (mapcar #'(lambda (b) (paste a b)) y))
  794.                           x))))
  795.       (t (cross-names x (apply #'cross-names args))))))
  796.