;
; Rate model for the ship damage data of McCullagh & Nelder (1989)
; "Generalized Linear Models", page 205
;

(require "glim")
(require "qvcalc")

(def aggregate-months-service 
  (list
   127 63 1095 1095 1512 3353 0 2244 
   44882 17176 28609 20370 7064 13099 0 7117
   1179 552 781 676 783 1948 0 274
   251 105 288 192 349 1208 0 2051
   45 0 789 437 1157 2161 0 542))

(def count-is-informative (> aggregate-months-service 0))

(def damage-incident-count (select
 (list
   0 0 3 4 6 18 0 11 
   39 29 58 53 12 44 0 18 
   1 1 0 1 6 2 0 1 
   0 0 0 0 2 11 0 4 
   0 0 7 7 5 12 0 1)
 (which count-is-informative)))

(def ship-type (select (repeat '(A B C D E) (repeat 8 5)) 
                  (which count-is-informative)))
(def year-built (select 
    (repeat (repeat '(60-64 65-69 70-74 75-79) (repeat 2 4)) 5) 
    (which count-is-informative)))
(def operating-period (select (repeat '(60-74 75-79) 20) 
                              (which count-is-informative)))
(def log-exposure (log (select aggregate-months-service 
                               (which count-is-informative))))

(defmeth poissonreg-proto :fit-scale () 
   (/ (sum (^ (send self :chi-residuals) 2)) (send self :df)))
; use mean X^2 rather than mean deviance to estimate overdispersion

(send poissonreg-proto :estimate-scale t) 
; for a "quasi-Poisson" model as reported by McCullagh & Nelder,
; with Poisson-based covariance matrix scaled up to allow for
; the apparent overdispersion

(send poissonreg-proto :epsilon 0.00000001)      ; use strict
(send poissonreg-proto :epsilon-dev 0.00000001)  ; convergence criteria

(def ship-model 
	(poissonreg-model 
	    (append (indicators ship-type) 
	            (indicators year-built) 
	            (indicators operating-period))
		damage-incident-count 
		:offset log-exposure
		:predictor-names 
		     (append (level-names ship-type :prefix 'type)
		             (level-names year-built :prefix 'built)
		             (level-names operating-period :prefix 'operated))
		:response-name 'damage-incident-count))

; Now define quasi-variance objects for the factors SHIP-TYPE and YEAR-BUILT

; Note that there is no point in quasi-variances for the third factor, 
; OPERATING-PERIOD, since it only has two levels and so a single standard 
; error suffices for the one and only contrast that exists.

(def ship-type-qvs (get-qvs ship-model "type"  
      :model-name "main effects model for ship damage rates"
      :factor-name "ship type"))
      
; Here we could alternatively have used (list 1 2 3 4) in place of the 
; prefix "type", since the four relevant coefficients appear in the 2nd thru 5th
; positions in the list for this model (and indexing in Lisp starts at 0). 

(send ship-type-qvs :summary)

(def ship-built-qvs (get-qvs ship-model "built" 
      :model-name "main effects model for ship damage rates"
      :factor-name "year built"))
      
; Here we could alternatively have used (list 5 6 7) in place of the 
; prefix "built". 

(send ship-built-qvs :summary)