src/lsp/glim.lsp

   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 = ~,6g~%"
 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~25t~13,6g~40t(~,6g)~%" (car coefs) (car se-s))
 398           (setf coefs (cdr coefs))
 399           (setf se-s (cdr se-s)))
 400     (dotimes (i (array-dimension x 1))
 401              (cond
 402                ((member i (send self :basis))
 403                 (format t "~a~25t~13,6g~40t(~,6g)~%"
 404                         (select p-names i) (car coefs) (car se-s))
 405                 (setf coefs (cdr coefs) se-s (cdr se-s)))
 406                (t (format t "~a~25taliased~%" (select p-names i)))))
 407     (format t "~%")
 408     (if (send self :estimate-scale)
 409         (format t "Scale Estimate:~25t~13,6g~%" (send self :scale))
 410         (format t "Scale taken as:~25t~13,6g~%" (send self :scale)))
 411     (format t "Deviance:~25t~13,6g~%" (send self :deviance))
 412     (format t "Number of cases:~25t~9d~%" (send self :num-cases))
 413     (if (/= (send self :num-cases) (send self :num-included))
 414         (format t "Number of cases used:~25t~9d~%" (send self :num-included)))
 415     (format t "Degrees of freedom:~25t~9d~%" (send self :df))
 416     (format t "~%")))
 417
 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))))))