;;------------------------------------------------------------------;; ;; ;; ;; Extra-Poisson Variation in Log-linear Models ;; ;; ;; ;; Luca Scrucca, Jan 2000 ;; ;;------------------------------------------------------------------;; (arc-setting (QUOTE *poisson-extra-variance*) (QUOTE T) "If T use Breslow method to estimate extra Poisson variance (see Breslow N.E. (1984), Extra-Poisson Variation in Log-linear Models, Applied Statistics, pp. 38-44)." (QUOTE T)) (defmeth poissonreg-proto :fit-extra-variation (&key (maxiter 30)) "Uses Breslow quasi likelihood model to estimate extra Poisson variance. The model assumes: E(y|x) = mu = exp{b'x} Var(y|x) = mu(1+phi*mu) If phi=0 there is no overdispersion. This algorithm provides an estimate of the extra-variation parameter phi, and the weights to use in subsequent modelling stages. (see Breslow N.E. (1984), Extra-Poisson Variation in Log-linear Models, Applied Statistics, pp. 38-44)" (unless (equal (send (send self :link) :slot-value 'proto-name) 'log-link) (message-dialog "Extra variation fit available only for poisson regression model with log link function") (return-from :fit-extra-variation)) (let* ((data (send self :data)) (pw (send self :pweights)) (inc (which (send self :included))) (y (select (send self :yvar) inc)) (mu (select (send self :fit-means) inc)) (h (diagonal (matmult (send self :x-matrix) (send self :xtxinv) (transpose (send self :x-matrix))))) (p (send self :num-coefs)) (n (send self :num-included)) ;; initial estimate of sigma^2 (sigma2 (div (- (send self :pearson-x2) (- n p)) (sum (* mu (- 1 (* mu h)))))) ) (send self :verbose nil) (format t "~&Extra-Poisson Variation in Log-linear Models~%Fitting .") ;; loop until X2 approx equal to 1 (do* ((i 1 (+ i 1)) (mu mu (select (send self :fit-means) inc)) (sigma2 sigma2 (div (sum (/ (^ (- y mu) 2) (* mu (+ mu (^ sigma2 -1))))) (- n p))) (w (^ (+ 1 (* sigma2 mu)) -1) (^ (+ 1 (* sigma2 mu)) -1)) ) ((< (- (abs (/ (send self :pearson-x2) (- n p))) 1) (send self :epsilon)) (format t "~&Converged after ~a iteration(s) ~%" i) (send self :add-slot 'extra-variation sigma2) (send data :poisson-wts w) (send self :args (append (list :weights "Extra-poisson-var") (send self :args))) (defmeth self :display (&rest args) (apply #'call-next-method args) (format t "Extra variation: ~10,3f~%" (send self :extra-variation))) (send self :display)) (when (> i maxiter) (format t "~&Algoritm not converged after ~a iterations!~%" maxiter) (send self :pweights pw) ;; restore initial pweights (return-from :fit-extra-variation)) (format t " .") (force-output) (send self :pweights w) (send self :compute) ) ) ) (defmeth poissonreg-proto :extra-variation (&optional (new nil set)) "Set or retrieves the extra-variation parameter estimate for poisson log-linear model (see :estimate-extra-variation method for references)." (when set (setf (slot-value 'extra-variation) new)) (slot-value 'extra-variation)) (defmeth datalist-proto :poisson-wts (&rest what) (apply #'send self :variate what)) (defmeth dataset-proto :poisson-wts (data) "Method args: (&rest args) Creates a new datalist of type poisson-wts." (send self :add-datalist (send datalist-proto :new :poisson-wts :data data :name "Extra-poisson-var")) (send self :find-datalist "Extra-poisson-var")) ;; Modify GLM menu for adding the extra variation fit (defmeth glim-proto :menu-items () (mapcar #'(lambda (a) (make-menu-item a self)) (append *arc-glm-regression-menu-items1* (if (and *binomial-extra-variance* (eql (send self :slot-value 'proto-name) 'binomialreg-proto)) (list 'extra-bin-var-menu-item)) (if (and *poisson-extra-variance* (eql (send self :slot-value 'proto-name) 'poissonreg-proto)) (list 'extra-poisson-var-menu-item)) *arc-glm-regression-menu-items2* *arc-linear-regression-menu-items2*))) (rc-menu-item 'extra-poisson-var-menu-item "Fit extra Poisson variance" :fit-extra-variation) ;; fun1.lsp: correction of div function for handling the case that a1 a2 are single values ;; (defun div (a1 &optional a2) (let* ((num (if a2 a1 (repeat 1 (length a1)))) (denom (if a2 a2 a1)) (result nil)) (cond ((eq *computer-type* 'windows) (unless (listp a1) (setf num (list num)) (setf denom (list denom))) (setf result (mapcar #'(lambda (n d) (if (is-missing d) not-a-number (/ n d))) num denom)) (if (listp a1) result (car result))) (t (/ num denom)))))