;; ;; Put this file in the Extras directory in the directory in ;; which you start Arc. ;; Source: http://www.stat.umn.edu/arc ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; Arc lisp code for fitting normal dispersion models ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; $Revision 1.3$ ;; ;; ;;; disp-identity-link is a duplicate of identity-link, persumably to allow ;;; for modifications that have not yet been included. ;;; (defproto disp-identity-link '() () identity-link) ;; ;; Normal dispersion models are fit using generalized inear models. ;; ;; First, check to see if normal dispersion models are already in the ;; Graph&Fit menu. If not, add an item to the menu: ;; (when (null (position 'norm-disp-menu-item *arc-fit-menu-items* :test #'eq)) (arc-glm-method :menu-item 'fit-normdisp-menu-item :menu-title "Fit normal-dispersion" :fitname :normdisp :mean-functions '("Identity" "Inverse" "Exponential") :link-functions (list disp-identity-link inverse-link log-link) :proto 'normal-dispersion-model)) ;; ;; The function arc-glm-method adds the new method to the other glim methods. ;; it assigns the menu item the name 'fit-normdisp-menu-item, and the title ;; "Fit normal-dispersion". The fitname is the dataset-proto method that is ;; called when the menu item is selected. The list of mean functions are ;; the mean functions that will appear on the list, and the link-functions ;; are the actual inverse mean functions that exist within the program. ;; Finally, the proto is prototype of the model that is fit. ;;; ;;; Dataset-proto methods ;;; ;; When an item in the Graph&Fit method is selected, a dataset-proto method ;; is executed. The dataset method is then responsible for doing the ;; appropriate action. (defmeth dataset-proto :normdisp () (send self :fit-norm-disp :normdisp)) (defmeth dataset-proto :fit-norm-disp (&optional (type :normdisp)) "Method args (&optional (type :normdisp)) This method is responsible for assembling the information necessary to get the dispersion-dialog, and if the dialog is correctly completed, calling the make-disperion method. The method is apparently set up for Binomial, Gamma or Poisson dispersion as well, but these don't seem to be tested." (send self :ones) (let* ((name (send self :get-model-name (case type (:binomial "BD") (:gamma "GD") (:poisson "PD") (:normdisp "ND") (t "M")))) (old-mod (send self :last-model :type type)) (names (send self :types (append '(:variate :factor :factor1 :interaction) (if (eql type :binomial) (list :ones) nil) (if (eql type :poisson) (list :quasi-indep :quasi-symmetry))))) (n (if old-mod (reverse (set-difference names (combine (send old-mod :predictors) (send old-mod :response-name)))) names)) (om (if old-mod (send old-mod :args) nil)) (ans (apply #'dispersion-dialog names :name name :candidates n :type type om))) (if ans (apply #'send self :make-dispersion ans)))) (defmeth dataset-proto :make-dispersion (&rest args &key predictors response variance-model intercept weights name offset trials type mean-function) "Method args: (&rest args) This method takes the specifications from the dispersion-dialog, and sets up the regression by send the :new method to the normal-dispersion-model proto." (let* ((x (apply #'append (mapcar #'(lambda (n) (send self :data n)) predictors))) (predictor-names (combine (mapcar #'(lambda (a) (send self :labels a)) predictors))) (z (apply #'append (mapcar #'(lambda (n) (send self :data n)) variance-model))) (z-names (combine (mapcar #'(lambda (a) (send self :labels a)) variance-model))) (y (first (send self :data response))) (response-name response) (pweights (first (send self :data weights))) (weights-name weights) (dataset self) (not-missing (apply #'send self :make-not-missing-indicator (append (if weights (list weights)) (list response) predictors))) (offset-name offset) (offset (first (send self :data offset))) (trials-name trials) (trials (first (send self :data trials))) (link (select (link-functions type) (position mean-function (mean-functions type) :test #'equal))) (r (send normal-dispersion-model :new :x x :y y :trials trials :predictor-names predictor-names :response-name response-name ;:pweights pweights :z z :z-names z-names :data self :link disp-identity-link :intercept intercept :not-missing not-missing :offset (if offset offset 0) :print nil))) (set (intern (string-upcase name)) r) (send r :type type) (send r :add-slot 'args args) (send r :delete-slot 'included) (defmeth r :included (&rest args) (let* ((inc (apply #'send dataset :included args)) (obs (send self :not-missing))) (if obs (mapcar #'(lambda (a b) (and a b)) inc obs) inc))) (send r :menu name) (send self :models r) (when *show-display-fit* (send r :display)) r)) ;;; ;;; This is the normal-dispersion model prototype. ;;; (defproto normal-dispersion-model '(z z-names var-model) () normalreg-proto) (defmeth normal-dispersion-model :isnew (&rest args &key z z-names) "Method args: (&rest args &key z z-names) Differs from the normalreg-proto :isnew by keeping the predictors and their names for the variance terms in z and z-names." (send self :z z) (send self :z-names z-names) (apply #'call-next-method args)) (defmeth normal-dispersion-model :z (&optional new) "Method args: (&optional new) Accessor function." (when new (setf (slot-value 'z) new)) (slot-value 'z)) (defmeth normal-dispersion-model :z-names (&optional new) "Method args: (&optional new) Accessor function." (when new (setf (slot-value 'z-names) new)) (slot-value 'z-names)) (defmeth normal-dispersion-model :compute () "Method args: () The compute method, based on a paper by Aitkin (198?)." (send self :scale 1) ; <-- set the scale = 1 (send self :estimate-scale nil) ; <-- turn off estimation of scale (send self :verbose nil) ; <-- turn off iteration printing (call-next-method) ; <-- estimate mean via OLS. (let* ((z (send self :z)) ; <-- z = variance predictors (itmax (send self :count-limit)) (n (send self :num-included)) (eps (send self :epsilon)) (g (gammareg-model z (^ (send self :raw-residuals) 2) :print nil :link log-link :verbose nil :response-name (send self :response-name) :included (send self :included) :predictor-names (send self :z-names)))) (send g :scale 2) (send g :estimate-scale nil) (send self :var-model g) ; <-- save g in the 'var-model slot (send g :compute) ; <-- estimate variance linear combination (send self :add-slot 'initial-deviance (* n (+ (log (* 2 pi)) 1 (log (/ (send self :residual-sum-of-squares) n))))) (do ((iter 1 (+ iter 1)) (fnew (send self :disp-dev) (send self :disp-dev)) (fold positive-infinity fnew)) ((or (= iter itmax) (> eps (abs (- fold fnew))))) ; <-- stopping crit. (format t "Iteration ~d: deviance = ~,6g~%" iter fnew) (send self :pweights (/ (exp (send g :eta)))) ; <-- working weights (call-next-method) ; <-- reestimate mean (send g :yvar (^ (send self :raw-residuals) 2)) (send g :compute) ))) (defmeth normal-dispersion-model :disp-dev () "Method args: () Compute the deviance for dispersion models." (let* ( (sel (which (send self :included))) (s (* (/ (send self :residual-sum-of-squares) (send self :num-included) ) (exp (select (send (send self :var-model) :eta) sel))))) (sum (+ (log (* 2 pi)) (log s) (/ (^ (select (send self :raw-residuals) sel) 2) s))))) (defmeth normal-dispersion-model :var-model (&optional new) (when new (setf (slot-value 'var-model) new)) (slot-value 'var-model)) (defmeth normal-dispersion-model :display () "Method args: () Display the results." (format t "Normal Dispersion Model: Model for the mean~%") (call-method regression-model-proto :display :anova nil) (format t "~%Normal Dispersion Model: Model for the variance~%") (send (send self :var-model) :display) (format t "-2*L(max), const var: ~23t~13,6g~%" (slot-value 'initial-deviance)) (format t "-2*L(max), model: ~23t~13,6g~%" (send self :disp-dev))) (defmeth normal-dispersion-model :variance-function () (* (send self :scale) (exp (send (send self :var-model) :eta)))) ;;; ;;; dialog box stuff ;;; (defun dispersion-dialog (&rest args) (let* ((ans (apply #'dd1 args))) (if ans (apply #'dd2 ans)))) (defun dd1 (names &key (type :normal) (name "reg") (offset "") (weights "") (trials "") (variance-model "") (response "") mean-function (candidates names) (predictors (repeat "" (length names))) (intercept t)) "Function args: (names &key reg-name type) Opens a dialog box to define a regression." (let* ((version (send text-item-proto :new (format nil "Arc ~a" *arc-version*))) (all (send allocate-items-proto :new)) (candidates-l (send subordinate-list-proto :new (pad-list candidates (length names)) 'primary all)) (candidates-n (send text-item-proto :new "Candidates")) (predictors-l (send subordinate-list-proto :new (pad-list predictors (length names)) 'secondary all)) (predictors-n (send text-item-proto :new "Pred. for mean")) (varmod (send choice-item-proto :new (list "Use mean predictors for variance" "Choose variance predictors in next dialog"))) (response-l (send subordinate-list-proto :new (list response) 'oneitem all)) (response-n (send text-item-proto :new "Response... ")) (offset-l (send subordinate-list-proto :new (if offset (list offset) '("")) 'oneitem all)) (offset-n (send text-item-proto :new "Offset... ")) (trials-l (when (eql type :binomial) (list (send subordinate-list-proto :new (if trials (list trials) '("")) 'oneitem all)))) (trials-n (when (eql type :binomial) (list (send text-item-proto :new "Trials... ")))) (intercept-n (send toggle-item-proto :new "Fit Intercept" :value intercept)) ; (mod-type (format nil "~a-dispersion" (string-capitalize type))) (mod-type (format nil "normal-dispersion")) (name-l (send edit-text-item-proto :new name :text-length 14 )) (name-n (send text-item-proto :new (format nil "Name for ~%~a" mod-type))) (mfs (mean-functions type)) (mfs-l (send choice-item-proto :new mfs :value (position mean-function mfs :test #'equal))) (mfs-n (send text-item-proto :new "Mean function")) (get-answer #'(lambda () (let* ((ans (send all :value)) (raw (send all :all-values))) (cond ((or (null (select ans 1)) (null (select ans 2)) (when trials-n (null (select ans 5)))) (message-dialog "Predictor, response or trials missing") (dispersion-dialog names :mod-type mod-type :name (send name-l :text) :candidates (select raw 0) :predictors (select raw 1) :varmod (select (list 'same 'own 'yhat) (send varmod :value)) :response (first (select raw 2)) :offset (first (select raw 3)) :trials (when trials-n (first (select raw 4))) :intercept (send intercept-n :value))) (t (list :predictors (select ans 1) :candidates (select ans 0) :type type :response (first (select ans 2)) :varmod (select (list 'same 'own 'yhat) (send varmod :value)) :offset (first (select ans 3)) :trials (when trials-n (first (select ans 4))) :mean-function (select mfs (send mfs-l :value)) :intercept (send intercept-n :value) :name (send name-l :text))))))) (cancel (send modal-button-proto :new "Cancel")) (done (send modal-button-proto :new "Done" :action #'(lambda () (funcall get-answer)))) (dialog (send modal-dialog-proto :new (list (list version name-n name-l) (list (append (list candidates-n candidates-l response-n offset-n) trials-n) (append (list predictors-n predictors-l response-l offset-l) trials-l) (list intercept-n mfs-n mfs-l done cancel)) (list varmod))))) (send dialog :default-button done) (send dialog :modal-dialog) )) (defun dd2 (&rest args &key response candidates predictors varmod variance-model) (cond ((eq varmod 'same) (append (list :variance-model predictors) args)) ((eq varmod 'yhat) (append (list :variance-model 'yhat) args)) (t (let* ((n (length (combine candidates predictors))) (ans (twolist-dialog candidates :left-title "Candidates" :left-initial (pad-list candidates n) :right-title "Variance model" :right-initial (pad-list predictors n) :title (format nil "Select predictors for variance model. Response is ~a" response)))) (cond ((null ans) nil) ((null (second ans)) (message-dialog "Specify at least one variable or Cancel") (apply #'dd2 args)) (t (append (list :variance-model (second ans)) args)))))))