;;
;;    GEEs for XLISP-Stat. Modelled on glim.lsp
;;
;;                Version 1.3.
;;
;; Reference: Zeger & Liang, Biometrics 1986
;;            Liang & Zeger, Biometrika 1986
;;            Preisser & Qaqish, Biometrika 1996 for the diagnostics
;;
;;      Copyright (c) thomas lumley 1995-1996
;;
;; Use and copying permitted under the terms of the 
;; GNU public license: see file COPYING
;;
;;
;; For documentation see the file gee.{tex,ps,dvi}
;;
;;

(require "glim") ;; the link functions from glim.lsp are needed

(require "modelformula") ;; model formula functions are needed

(provide "gee")

;;
;; additional methods :valideta :validmu for link function.
;;
;;  :valideta eta  tests whether the vector eta lies entirely in the domain
;;                 of the link function
;;  :validmu  mu   tests whether the vector mu lies entirely in the range 
;;                 of the link function

(defmeth logit-link :valideta (eta) t)

(defmeth probit-link :valideta (eta) t)

(defmeth log-link :valideta (eta) t)

(defmeth cloglog-link :valideta (eta) t)

(defun power-link (k) (send power-link-proto :new k))

(defmeth power-link-proto :valideta (eta) (> (min eta) 0))

(defmeth identity-link :valideta (eta) t)

(defmeth inverse-link :valideta (eta) (> (min eta) 0))

(defmeth sqrt-link :valideta (eta) (> (min eta) 0))
;
(defmeth logit-link :validmu (mu) (and (< (max mu) 1) (> (min mu) 0)))

(defmeth probit-link :validmu (mu) (and (< (max mu) 1) (> (min mu) 0)))

(defmeth cloglog-link :validmu (mu) (and (< (max mu) 1) (> (min mu) 0)))

(defmeth log-link :validmu (mu) (> (min mu) 0))

(defmeth power-link-proto :validmu (mu) (> (min mu) 0))

(defmeth identity-link :validmu (mu) t)

(defmeth inverse-link :validmu (mu) (> (min mu) 0))

(defmeth sqrt-link :validmu (mu) t)

;;
;; Miscellaneous useful functions
;;

(defun unique (x)
"Arguments: x
Returns the unique elements of a sequence x in no very obvious order"
 (remove-duplicates x :test #'equalp))

(defun argselect (xs x) 
"Arguments: xs x
returns indices where list xs = element x"
(which (= x xs))) 

(defun argrank (x) 
"Arguments: x
returns a sequence such that (select x (argrank x)) is in ascending order"
(elt  (column-list (apply #'bind-rows (sort (row-list (bind-columns (copy-seq x) (iseq (length x)))) #'(lambda (r s) (< (elt r 0) (elt  s 0)))))) 1))


(defun zero-matrix (k)
"Arguments: k
Returns a kxk zero-matrix"
 (* 0 (identity-matrix k)))

(defun one-matrix (k)
"Arguments: k
Returns a kxk matrix of 1s"
 (outer-product (repeat 1 k) (repeat 1 k)))

(defun logg (x) (if (compound-data-p x) (map-elements #'logg x)
		  (if (> x 0) (log x) 0)))

(defun band-matrix (width mat) 
"Arguments: width mat
Zeroes elements of mat to make it band-diagonal. 
width is the half-width (width=0 is diagonal)"
  (let* ((mdim (length (row-list mat))))
	 (map-elements #'(lambda (x y) (cond (x y) (t 0)))
		       (<= (abs (outer-product (iseq 1 mdim) (iseq 1 mdim) #'-)) width)
		       mat)
	 ))

(defun cat (ll) (apply #'append (mapcar #'(lambda (x) (coerce x 'list)) ll)))

(defun band-average (mat) 
"Arguments: mat
Averages elements of mat parallel to the main diagonal for symmetric matrix"
(let* ((cols (column-list mat))
       (collist (cat cols))
       (p (length cols))
       (steps (* (iseq 0 (- p 1)) (+ p 1)))
       (newlist (* 0 collist))
       (junk (dotimes (i p newlist)
		(let* (
		       (thesesteps (+ i  (select steps (iseq 0 (- (- p 1) i))))) 
		       )
		  (setf (select newlist thesesteps) (repeat (mean (select collist thesesteps)) (- p i))))
		))
       (newmat (make-array (list p p) :displaced-to (coerce newlist 'vector)))
       )
   (+ newmat (transpose newmat) (- (diagonal (diagonal newmat))))
))


(defun matrix-trace (matrix)
"Arguments: matrix
Returns its trace (of course)"
 (sum (diagonal matrix)))


(defun sign (x) (cond ((= x 0) 0)
		      ((< x 0) -1)
		      ((> x 0) +1)
		      (t (error "Can't happen in (sign)"))
		      ))


;;
;; some diagnostic plot functions
;; and extra methods for scatterplot-proto
;;

(defmeth scatterplot-proto :group-list (&optional (new nil set))
"Adds or returns group list"
(when set
      (setf (slot-value 'group-list) new))
(slot-value 'group-list)
)

(defmeth scatterplot-proto :do-group-hilite (x y m1 m2)
"Highlights all points in the same group as the selected one
Only sensible when called by the appropriate mouse mode"
(send self :points-selected nil)
(let* ((cr (send self :click-range))
       (groups (send self :group-list))
       (p (first (send self :points-in-rect
		       (- x (round (/ (first cr) 2)))
		       (- y (round (/ (second cr) 2)))
		       (first cr)
		       (second cr)))))
  (if (and p groups)
      (let* ((pgroup (elt groups p))
	     (thisgroup (argselect groups pgroup))
	     )
	(send self :points-selected thisgroup)
	))
))


(defmeth scatterplot-proto :add-group-hilite (grouplist)
"add a Group Highlight mode"
(send self :add-mouse-mode 'group-hilite
      :title "Highlight Groups"
      :cursor 'finger
      :click :do-group-hilite
)
(send self :add-slot 'group-list)
(send self :group-list grouplist)
)


(defun scatter-smooth (x y &rest graphargs &key (symmetric nil) (span .6)) 
"Arguments: x y &rest graphargs &key (symmetric nil) (span .6)
Draws a scatterplot with a lowess smoother. If symmetric is t uses the
two-iteration robust version (not good for skewed residuals).
Any extra arguments are passed to plot-points"
  (let* (
	(steps (if symmetric 2 1))
	(smooth (lowess x y :f span :steps steps))
	(plot (apply #'plot-points (append (list x y) graphargs)))
	)
    (send plot :add-lines (first smooth)  (second smooth) )
    plot
  )
  
)

(defun index-plot (y &rest graphargs)
"Arguments: y &rest graphargs
Plots y against 1:(length y). Other arguments are passed to plot-points"
  (apply #'plot-points (append (list (iseq 1 (length y)) y) graphargs)))


;;;;
;;;;                Variance Function Prototypes
;;;;
;;;;
;;;; Variance functions are objects responding to five messages:
;;;;
;;;;    :varweight (mu)  Return the variance function evaluated at
;;;;          the sequence of means mu
;;;;
;;;;    :pearson-resid (y mu) Return (y-mu)/sqrt(V(mu))
;;;;
;;;;    :quasideviance (y mu) the quasideviance (ignoring correlation)
;;;;
;;;;    :validmu (mu) tests if the vector mu lies entirely in the permitted
;;;;        range
;;;;
;;;;    :prettyprint displays the variance function nicely
;;;;
;;;;  the prototype is a Normal variance function
;;;; 

(defproto gee-varfun-proto)

(defmeth gee-varfun-proto :varweight (mu) (repeat 1 (length mu)))

(defmeth gee-varfun-proto :pearson-resid (y mu) (- y mu))

(defmeth gee-varfun-proto :quasi-deviance (y mu) (sum (^ (- y mu) 2)))

(defmeth gee-varfun-proto :validmu (mu) t)

(defmeth gee-varfun-proto :prettyprint (&optional (stream t))
  (format stream "~%Variance function:~25tPrototype"))

(defmeth gee-varfun-proto :print (&optional (stream t))
  (format stream "#<GEE Variance Object: ~s>" (slot-value 'proto-name)))


;;;;
;;;; Normal error structure
;;;;

(defproto normal-error () () gee-varfun-proto)

(defproto gaussian-error () () gee-varfun-proto)  ;;;synonyms are nice

(defmeth normal-error  :prettyprint (&optional (stream t))
  (format stream "~%Variance function:~25tNormal error"))
(defmeth gaussian-error :prettyprint (&optional (stream t))
  (format stream "~%Variance function:~25tGaussian error"))

;;;;
;;;; Binomial error structure
;;;;

(defproto binomial-error () () gee-varfun-proto)

(defmeth binomial-error :varweight (mu)  (* mu (- 1 mu)))
(defmeth binomial-error :pearson-resid (y mu) (/ (- y mu) (sqrt (send self :varweight mu))))
(defmeth binomial-error :quasi-deviance (y mu) (* 2 (sum (+ (* (logg (/ y mu)) y) (* (logg (/ (- 1 y) (- 1 mu))) (- 1 y)) ))))

(defmeth binomial-error :validmu (mu) (and (< (max mu) 1) (> (min mu) 0)))

(defmeth binomial-error :prettyprint (&optional (stream t))
  (format stream "~%Variance function:~25tBinomial error"))


;;;;
;;;; Poisson error structure
;;;;
;;;;
(defproto poisson-error () () gee-varfun-proto)

(defmeth poisson-error :varweight (mu) (copy-list mu))
(defmeth poisson-error :pearson-resid (y mu) (/ (- y mu) (sqrt mu)))
(defmeth poisson-error :quasi-deviance (y mu) (* 2 (sum (- (* y (logg (/ y mu))) (- y mu)))))
(defmeth poisson-error :validmu (mu) (> (min mu) 0))

(defmeth poisson-error :prettyprint (&optional (stream t))
  (format stream "~%Variance function:~25tPoisson error"))


;;;;
;;;; Gamma error structure
;;;;
;;;;
(defproto gamma-error () () gee-varfun-proto)

(defmeth gamma-error :varweight (mu) (^ mu 2))
(defmeth gamma-error :pearson-resid (y mu) (/ (- y mu)  mu))
(defmeth gamma-error :quasi-deviance (y mu) (* 2 (sum (- (/ (- y mu) mu) (log (/ y mu))))))
(defmeth gamma-error :validmu (mu) (> (min mu) 0))

(defmeth gamma-error :prettyprint (&optional (stream t))
  (format stream "~%Variance function:~25tGamma error"))






;;;
;;; Power variance functions
;;;
;;;

(defun power-error (k) (send power-error-proto :new :k k))

(defproto power-error-proto '(k) '() gee-varfun-proto)

(defmeth power-error-proto :varweight (mu) (^ mu (send self :k)))

(defmeth power-error-proto :pearson-resid (y mu)  (/ (- y mu) (sqrt (send self :varweight mu))))

(defmeth power-error-proto :validmu (mu) (> (min mu) 0))

(defmeth power-error-proto :isnew (&key k) (send self :k k))

(defmeth power-error-proto :k (&optional (new nil set))
"Message args: (&optional new)
Sets or returns variance power parameter."
  (when set
	(setf (slot-value 'k) new)
	)
  (slot-value 'k))

(defmeth power-error-proto :prettyprint (&optional (stream t))
  (format stream "~%Variance function:~25tmu^~6g" (send self :k)))

(defmeth power-error-proto :quasi-deviance (y mu)  
 (let* ((k (send self :k))
	 (denom (* (- 1 k) (- 2 k)))
	 )
   (cond ((= k 1) (send poisson-error :quasi-deviance y mu))
	 ((= k 2) (send gamma-error :quasi-deviance y mu))
	 (t  (* -2 (/ (sum  (^ y (- 2 k)) (* (^ mu (- 1 k)) (- (* (- 1 k) y) (* (- 2 k) mu)))))))
	 )
   )
)
;;;;
;;;; default link-variance pairs (an association list)
;;;;

(def default-link `((,binomial-error . logit-link)
		    (,normal-error . identity-link)
		    (,gaussian-error . identity-link)
		    (,poisson-error . log-link)
		    (,gamma-error . inverse-link)
		    ))

(def default-error  `((,logit-link . binomial-error)
		      (,probit-link . binomial-error)
		      (,cloglog-link . binomial-error)
		      (,identity-link . normal-error)
		      (,log-link . poisson-error)
		      (,inverse-link . gamma-error)
		      ))





;;;;
;;;;               Correlation Structure Prototypes
;;;;
;;;;
;;;; Correlation structures are objects responding to three messages:
;;;;
;;;;    :rinv (k alpha times )  returns the inverse of the
;;;;           correlation matrix for a group with k observations, 
;;;;           correlation parameters alphas and times/types times.
;;;;           Some correlation structures would require 
;;;;          more info, such as a GEE1 type structure needing betas and predictors.
;;;;    :estimate (r group  times allow-missing) uses Pearson residuals r,
;;;;          group identifiers  group and time/type numbers t to
;;;;          estimate alpha and the scale parameter phi and returns
;;;;          (alpha phi morealpha) where morealpha is anything else you
;;;;          might want to  know about the correlation structure.
;;;;          eg for the exchangeable structure it gives the estimated
;;;;          correlation in each group. 
;;;;          In fact you can put anything in this list provided phi is the
;;;;          second element as all the messages which need alpha pass back the
;;;;          whole thing. As a general rule the first component should be the
;;;;          thing that is printed out by :prettyprint and should contain
;;;;          enough information to fit the model.
;;;;          If allow-missing is nil the structure may stop with an error 
;;;;          message if it detects missing observations and it cares. 
;;;;          With allow-missing t
;;;;          the correlations will be estimated assuming missing completely at
;;;;          random.
;;;;          Exchangeable and independence structures don't care about missing
;;;;          observations. The distinction is that the natural complete-data
;;;;          algorithm for estimating the correlation still functions when
;;;;          are missing.
;;;;
;;;;    :prettyprint (alpha stream) Print out a nice description of the
;;;;          working model and its parameter estimates into stream.  
;;;; 
;;;;    The prototype is the independence structure
;;;;


(defproto gee-corr-proto)

(defmeth gee-corr-proto :rinv (k alpha times ) (identity-matrix k))

(defmeth gee-corr-proto :estimate (r g times allow-missing) (list nil (/ (length r)(sum (^ r 2)) )))

(defmeth gee-corr-proto :prettyprint (&optional  (alpha nil) (stream t) ) 
  (format stream "~% Working Model Prototype~%"))

;;
;;independence
;;

(defproto independence-corr () () gee-corr-proto)
(defmeth independence-corr :prettyprint (&optional  (alpha nil) (stream t)) 
  (format stream "~%Independence Working Model~%"))

;;;
;;; user-specified fixed correlation matrix
;;;
;;; note: this estimates the scale parameter as 1 to save effort. This has no
;;; effect on the parameters or their standard errors but may be disconcerting.
;;; to the user.

(defun fixed-corr (rmat)  (send fixed-corr-proto :new :rmat rmat))

(defproto fixed-corr-proto '(rmat) '() gee-corr-proto)

(defmeth fixed-corr-proto :rmat (&optional (new nil set))
  (when set 
	(setf (slot-value 'rmat) new))
  (slot-value 'rmat)
)

(defmeth fixed-corr-proto :rinv (k alpha times)
    (let* (
	   (decoder (first (third alpha)))
	   (indices (mapcar #'(lambda (x) (cdr (assoc x decoder))) times))
	   )
       (if (= (length indices) (length decoder)) (select (second (third alpha)) indices indices)  (inverse (select (first alpha) indices indices)))))

(defmeth fixed-corr-proto :estimate  (r g times allow-missing)
 (let* ( 
	(alltimes (if (null times) 
		      (error "You must specify times for fixed correlation structure") 
		    (unique times)))
	(nn (mapcar #'(lambda (i) (count i times)) alltimes))
	(ntimes (cond (allow-missing (length alltimes))
		      ( (apply #'= nn) (length alltimes))
		      (t (error "Missing (unbalanced) data found with :allow-missing nil"))
		      ))
	(ranks (rank alltimes))
	(decoder (mapcar #'(lambda (x) (cons (elt x 0) (elt x 1))) (row-list (bind-columns alltimes ranks))))
	)
 (list (send self :rmat) (/ (length r) (sum (^ r 2))) (list decoder (inverse (send self :rmat))))
))


(defmeth fixed-corr-proto  :prettyprint (&optional  (alpha nil) (stream t) ) 
  (format stream "~%Known Correlation  Working Model: ~% correlation matrix =~%" )
  (print-matrix (first alpha) stream :float-digits *gee-corr-digits*)
)



;;;
;;;     exchangeable correlation structure
;;;


(defproto exchangeable-corr () () gee-corr-proto)

(defmeth exchangeable-corr :rinv (k alpha times) 
  (let* ( (alphaa (first alpha))
	 (aa (+ 1 (* (- k 1) alphaa)))
	 (R1 (- (identity-matrix k) (* (/ alphaa aa) (one-matrix k))))
	 )
    (/ R1 (- 1 alphaa)))
)

(defun exchcorr-one-group (ri) 
"Correlation and variance for exchangeable structure
from one group of Pearson residuals"
 (let* ( (rri (outer-product ri ri))
	 (rii (sum (diagonal rri)))
	 (rij (- (sum rri) rii))
	 (ni (length ri))
	 (onlyone (if (= ni 1) 1 0))
	 )
    (list  (if (= ni 1) 0 (/ rij (* ni (- ni 1)))) (/ rii ni) onlyone) 
 ))



(defmeth exchangeable-corr :estimate (r g times allow-missing)
  (let* (  (grouplist (unique g))
	   (n (length grouplist))
	  (corrs nil)
	  (corrlist (dolist (gi (coerce grouplist 'list) corrs)
		  (setf corrs (append  corrs (list (exchcorr-one-group (select r (argselect g gi)))) ))
		  ))
	  (corrtotal (apply #'+ corrlist))
	  (phi (/ n (second corrtotal)))
	  (alpha (/ (* phi (first corrtotal)) (- n (third corrtotal)))) 
	  (groupcorrs (mapcar #'(lambda (x) (if (equalp (first x) 0) alpha (/ (first x) (second x) ))) corrlist)) 
	  )
    (list alpha phi groupcorrs ))
)

(defmeth exchangeable-corr :prettyprint (&optional  (alpha nil) (stream t) ) 
  (format stream "~%Exchangeable Working Model:  correlation = ~7,4f~%" (first alpha)))

;;
;; clustered-corr as an alias for exchangeable-corr
;;

(defproto clustered-corr () () exchangeable-corr)
(defmeth clustered-corr  :prettyprint (&optional  (alpha nil) (stream t) ) 
  (format stream "~%Clustered (exchangeable) Working Model:  correlation = ~7,4f~%" (first alpha)))


;;;
;;;  saturated correlation structure
;;;

; note. To avoid redundant computation we store the inverse of the correlation 
; matrix in (third alpha) and the matrix itself in (first alpha). The 
; uninverted matrix is only used for printing out. We also store an 
; association list with all the times and their ranks
; in (third alpha) so that we can find the right submatrix
; when there are missing observations

(defproto saturated-corr () () gee-corr-proto)

(defmeth saturated-corr  :prettyprint (&optional  (alpha nil) (stream t) ) 
  (format stream "~%Saturated Working Model: ~% correlation matrix =~%" )
  (print-matrix (first alpha) stream :float-digits *gee-corr-digits*)
)

; This allows missing data, but only if the flag :allow-missing is set
; Small amounts of missing data are ok. The data must be missing completely at
; random for the method to be valid.

(defmeth saturated-corr :estimate  (r g times allow-missing)
 (let* ( 
	(alltimes (if (null times) 
		      (error "You must specify times for saturated correlation structure") 
		    (unique times)))
	(nn (mapcar #'(lambda (i) (count i times)) alltimes))
	(ntimes (cond (allow-missing (length alltimes))
		      ( (apply #'= nn) (length alltimes))
		      (t (error "Missing (unbalanced) data found with :allow-missing nil"))
		      ))
	(ranks (rank alltimes))
	(decoder (mapcar #'(lambda (x) (cons (elt x 0) (elt x 1))) (row-list (bind-columns alltimes ranks))))
	(grouplist (coerce (unique g) 'list))
	(ngroups (length grouplist))
	(ssp-matrix (zero-matrix ntimes))
	(count-matrix (zero-matrix ntimes))
	(junk (dolist (gg grouplist)
	   (let* ( (thisgroup (argselect g gg))
		   (indices  (mapcar #'(lambda (x) (cdr (assoc x decoder))) (select times thisgroup)))
		   (ri (select r thisgroup))
		   (ssp-inc (outer-product ri ri))
		   (ssp-subs (+ ssp-inc (select ssp-matrix indices indices)))
		   (count-inc (one-matrix (length thisgroup)))
		   (count-subs (+ count-inc (select count-matrix indices indices)))
		   )
	     (setf (select ssp-matrix indices indices) ssp-subs)
	     (setf (select count-matrix indices indices) count-subs)
	     )
	   ))
	(cov-matrix (/ ssp-matrix  count-matrix))
	(vars (diagonal cov-matrix))
	(phi (/ (mean vars)))
	(corrs (* cov-matrix phi))
	(corr-matrix (+ corrs (diagonal (- (repeat 1 ntimes) (diagonal corrs)))))
	)
   (list corr-matrix phi (list decoder (inverse corr-matrix)))
  ))


(defmeth saturated-corr :rinv (k alpha times)
    (let* (
	   (decoder (first (third alpha)))
	   (indices (mapcar #'(lambda (x) (cdr (assoc x decoder))) times))
	   )
       (if (= (length indices) (length decoder)) (select (second (third alpha)) indices indices) (inverse (select (first alpha) indices indices)))))



;;
;;  saturated correlation structure - ML version
;;
;;  this differs from the above in the handling of the diagonal matrix
;;
;;  for Normal data it is ML

; note. To avoid redundant computation we store the inverse of the correlation 
; matrix in (third alpha) and the matrix itself in (first alpha). The 
; uninverted matrix is only used for printing out. We also store an 
; association list with all the times and their ranks
; in (third alpha) so that we can find the right submatrix
; when there are missing observations

(defproto saturated-ml-corr () () gee-corr-proto)

(defmeth saturated-ml-corr  :prettyprint (&optional  (alpha nil) (stream t) ) 
  (format stream "~%Saturated Working Model (Normal-ML): ~% correlation matrix =~%" )
  (print-matrix (first alpha) stream :float-digits *gee-corr-digits*)
)

; This allows missing data, but only if the flag :allow-missing is set
; Small amounts of missing data are ok. The data must be missing completely at
; random for the method to be valid.

(defmeth saturated-ml-corr :estimate  (r g times allow-missing)
 (let* ( 
	(alltimes (if (null times) 
		      (error "You must specify times for saturated correlation structure") 
		    (unique times)))
	(nn (mapcar #'(lambda (i) (count i times)) alltimes))
	(ntimes (cond (allow-missing (length alltimes))
		      ( (apply #'= nn) (length alltimes))
		      (t (error "Missing (unbalanced) data found with :allow-missing nil"))
		      ))
	(ranks (rank alltimes))
	(decoder (mapcar #'(lambda (x) (cons (elt x 0) (elt x 1))) (row-list (bind-columns alltimes ranks))))
	(grouplist (coerce (unique g) 'list))
	(ngroups (length grouplist))
	(ssp-matrix (zero-matrix ntimes))
	(count-matrix (zero-matrix ntimes))
	(junk (dolist (gg grouplist)
	   (let* ( (thisgroup (argselect g gg))
		   (indices  (mapcar #'(lambda (x) (cdr (assoc x decoder))) (select times thisgroup)))
		   (ri (select r thisgroup))
		   (ssp-inc (outer-product ri ri))
		   (ssp-subs (+ ssp-inc (select ssp-matrix indices indices)))
		   (count-inc (one-matrix (length thisgroup)))
		   (count-subs (+ count-inc (select count-matrix indices indices)))
		   )
	     (setf (select ssp-matrix indices indices) ssp-subs)
	     (setf (select count-matrix indices indices) count-subs)
	     )
	   ))
	(cov-matrix (/ ssp-matrix  count-matrix))
	(vars (diagonal cov-matrix))
	(phi (/ (mean vars)))
	(corr-matrix (* cov-matrix phi))
	)

   (list corr-matrix phi (list decoder (inverse corr-matrix)))
  ))

(defmeth saturated-ml-corr :rinv (k alpha times)
    (let* (
	   (decoder (first (third alpha)))
	   (indices (mapcar #'(lambda (x) (cdr (assoc x decoder))) times))
	   )
       (if (= (length indices) (length decoder)) (select (second (third alpha)) indices indices) (inverse (select (first alpha) indices indices)))))


;;;
;;;  nonstationary m-dependence correlation structure
;;;
;;;  (first alpha) contains the matrix of estimated correlations
;;;  (third alpha) contains a decoder for the times and the inverse correlation matrix
;;;

(defun m-dependence (m &key (stationary nil)) (if stationary (send stat-m-dependence-proto :new :m m) (send m-dependence-proto :new :m m)))

(defproto m-dependence-proto '(m) '() gee-corr-proto)

(defmeth m-dependence-proto :m (&optional (new nil set))
"Message args: (&optional new)
Sets or returns m-dependence parameter."
  (when set
	(setf (slot-value 'm) new)
	)
  (slot-value 'm))

(defmeth m-dependence-proto :prettyprint (&optional  (alpha nil) (stream t) ) 
  (format stream "~%Nonstationary ~2d - dependence Working Model: ~% correlation matrix =~%" (send self :m) )
  (print-matrix (first alpha) stream :float-digits *gee-corr-digits*))

(defmeth m-dependence-proto :estimate (r g times allow-missing)
 (let* ( 
	(alltimes (if (null times) 
		      (error "You must specify times for nonstationary m-dependence  structure") 
		    (unique times)))
	(nn (mapcar #'(lambda (i) (count i times)) alltimes))
	(ntimes (cond (allow-missing (length alltimes))
		      ( (apply #'= nn) (length alltimes))
		      (t (error "Missing (unbalanced) data found with :allow-missing nil"))
		      ))
	(ranks (rank alltimes))
	(decoder (mapcar #'(lambda (x) (cons (elt x 0) (elt x 1))) (row-list (bind-columns alltimes ranks))))
	(grouplist (coerce (unique g) 'list))
	(ngroups (length grouplist))
	(ssp-matrix (zero-matrix ntimes))
	(count-matrix (zero-matrix ntimes))
	(junk (dolist (gg grouplist)
	   (let* ( (thisgroup (argselect g gg))
		   (indices  (mapcar #'(lambda (x) (cdr (assoc x decoder))) (select times thisgroup)))
		   (ri (select r thisgroup))
		   (ssp-inc (outer-product ri ri))
		   (ssp-subs (+ ssp-inc (select ssp-matrix indices indices)))
		   (count-inc (one-matrix (length thisgroup)))
		   (count-subs (+ count-inc (select count-matrix indices indices)))
		   )
	     (setf (select ssp-matrix indices indices) ssp-subs)
	     (setf (select count-matrix indices indices) count-subs)
	     )
	   ))
	(cov-matrix (/ ssp-matrix  count-matrix))
	(vars (diagonal cov-matrix))
	(phi (/ (mean vars)))
	(corrs (* cov-matrix phi))
	(corr-matrix (band-matrix (send self :m) (+ corrs (diagonal (- (repeat 1 ntimes) (diagonal corrs))))))
	)
   (list corr-matrix phi (list decoder (inverse corr-matrix)))
  )) 
	 

(defmeth m-dependence-proto :rinv (k alpha times)
    (let* (
	   (decoder (first (third alpha)))
	   (indices (mapcar #'(lambda (x) (cdr (assoc x decoder))) times))
	   )
       (if (= (length indices) (length decoder)) (select (second (third alpha)) indices indices) (inverse (select (first alpha) indices indices)))))

;;;
;;;  Stationary m-dependence correlation structure
;;;
;;;  (first alpha) contains the matrix of estimated correlations
;;;  (third alpha) contains a decoder for the times and the inverse correlation matrix
;;;

(defun stat-m-dependence (m) (send stat-m-dependence-proto :new :m m))

(defproto stat-m-dependence-proto '(m) '() gee-corr-proto)

(defmeth stat-m-dependence-proto :m (&optional (new nil set))
"Message args: (&optional new)
Sets or returns m-dependence parameter."
  (when set
	(setf (slot-value 'm) new)
	)
  (slot-value 'm))

(defmeth stat-m-dependence-proto :prettyprint (&optional  (alpha nil) (stream t) ) 
  (format stream "~%Stationary ~2d - dependence Working Model: ~% correlation matrix =~%" (send self :m) )
  (print-matrix (first alpha) stream :float-digits *gee-corr-digits*))


(defmeth stat-m-dependence-proto :estimate (r g times allow-missing)
 (let* ( 
	(alltimes (if (null times) 
		      (error "You must specify times for stationary m-dependence  structure") 
		    (unique times)))
	(nn (mapcar #'(lambda (i) (count i times)) alltimes))
	(ntimes (cond (allow-missing (length alltimes))
		      ( (apply #'= nn) (length alltimes))
		      (t (error "Missing (unbalanced) data found with :allow-missing nil"))
		      ))
	(ranks (rank alltimes))
	(decoder (mapcar #'(lambda (x) (cons (elt x 0) (elt x 1))) (row-list (bind-columns alltimes ranks))))
	(grouplist (coerce (unique g) 'list))
	(ngroups (length grouplist))
	(ssp-matrix (zero-matrix ntimes))
	(count-matrix (zero-matrix ntimes))
	(junk (dolist (gg grouplist)
	   (let* ( (thisgroup (argselect g gg))
		   (indices  (mapcar #'(lambda (x) (cdr (assoc x decoder))) (select times thisgroup)))
		   (ri (select r thisgroup))
		   (ssp-inc (outer-product ri ri))
		   (ssp-subs (+ ssp-inc (select ssp-matrix indices indices)))
		   (count-inc (one-matrix (length thisgroup)))
		   (count-subs (+ count-inc (select count-matrix indices indices)))
		   )
	     (setf (select ssp-matrix indices indices) ssp-subs)
	     (setf (select count-matrix indices indices) count-subs)
	     )
	   ))
	(cov-matrix (band-matrix (send self :m) (/ ssp-matrix  count-matrix)))
	(vars (diagonal cov-matrix))
	(phi (/ (mean vars)))
	(corr-matrix (* (band-average  cov-matrix) phi))
       	)
   (list corr-matrix phi (list decoder (inverse corr-matrix)))
  )) 
	 

(defmeth stat-m-dependence-proto :rinv (k alpha times)
    (let* (
	   (decoder (first (third alpha)))
	   (indices (mapcar #'(lambda (x) (cdr (assoc x decoder))) times))
	   )
       (if (= (length indices) (length decoder)) (select (second (third alpha)) indices indices) (inverse (select (first alpha) indices indices)))))

;;;
;;; AR-1 correlation structure
;;;
;;; (first-order autoregressive)


(defproto AR-1-corr () () gee-corr-proto)


(defmeth AR-1-corr :prettyprint (&optional  (alpha nil) (stream t) ) 
  (format stream "~%Autoregressive (AR-1) Working Model: correlation ~7,4f~%" (first alpha))
)

(defmeth AR-1-corr :estimate (r g times allow-missing)
  (let* (
	 (stat1 (send (stat-m-dependence 1) :estimate r g times allow-missing))
	 (rho (select (first stat1) 1 0))
	 (phi (second stat1))
	 (ntimes (length (unique times)))
	 (decoder (first (third stat1)))
	 (corr-matrix (^ rho (abs (outer-product (iseq 1 ntimes) (iseq 1 ntimes) #'-))))
	 )
    (list rho phi (list decoder (inverse corr-matrix) corr-matrix))
    )
)


(defmeth AR-1-corr :rinv (k alpha times)
    (let* (
	   (decoder (first (third alpha)))
	   (indices (mapcar #'(lambda (x) (cdr (assoc x decoder))) times))
	   )
       (if (= (length indices) (length decoder)) (select (second (third alpha)) indices indices) (inverse (select (third (third alpha)) indices indices)))))




;;;;
;;;;
;;;;                The General GEE Prototype
;;;;
;;;;

(defproto gee-proto 
  '(y x  times  group link varfun corr-structure offset pweights
      alpha beta scale allow-missing
	 epsilon epsilon-dev count-limit verbose recycle
	 eta quasideviance covariance-matrix naive-covariance-matrix needs-computing
	 intercept response-name predictor-names group-names case-labels formula scores
	 working-residuals cluster-influences weight-blocks
	 group-plots observation-plots
	 finite-sample-correction)
  '() 
  *object*
  "GEE -- Correlated data marginal GLM")

;; defined as an *object* because the fitting method is different from the 
;; linear model -- iteratively reweighted least squares with non-diagonal 
;; weight matrix is too slow and complicated. 

;;;;
;;;; Slot Accessors
;;;;

(defmeth gee-proto :y (&optional (new nil set))
"Message args: (&optional new)
Sets or returns dependent variable."
  (when set
	(setf (slot-value 'y) new)
	(send self :needs-computing t))
  (slot-value 'y))

(defmeth gee-proto :alpha (&optional (new nil set))
"Message args: (&optional new)
Sets or returns correlation parameters."
  (when set
	(setf (slot-value 'alpha) new)
	(send self :needs-computing t))
 (slot-value 'alpha))

(defmeth gee-proto :needs-computing (&optional (new nil set))
  (when set (setf (slot-value 'needs-computing) new))
  (slot-value 'needs-computing))

(defmeth gee-proto :intercept (&optional (val nil set))
"Message args: (&optional new-intercept)
With no argument returns T if the model includes an intercept term, nil if
not. With an argument NEW-INTERCEPT the model is changed to include or
exclude an intercept, according to the value of NEW-INTERCEPT."
  (when set 
	(setf (slot-value 'intercept) val)
	(send self :needs-computing t))
  (slot-value 'intercept))

(defmeth gee-proto :times (&optional (new nil set))
"Message args: (&optional new)
Sets or returns time column"
  (when set
	(setf (slot-value 'times) new)
	(send self :needs-computing t))
  (slot-value 'times))

(defmeth gee-proto :scale (&optional (new nil set))
"Message args: (&optional new)
Sets or returns scale parameter."
  (when set
	(setf (slot-value 'scale) new)
	)
  (slot-value 'scale))

(defmeth gee-proto :allow-missing (&optional (new nil set))
"Message args: (&optional new)
Flag to control missing data in correlation estimates.
When nil the program will not estimate correlations if
it detects missing data."
  (when set
	(setf (slot-value 'allow-missing) new)
	)
  (slot-value 'allow-missing))

(defmeth gee-proto :xnames (&optional (new nil set))
"Message args: (&optional new)
Sets or returns predictor names ."
  (when set
	(setf (slot-value 'xnames) new))
  (slot-value 'xnames))

(defmeth gee-proto :ynames (&optional (new nil set))
"Message args: (&optional new)
Sets or returns dependent variable name."
  (when set
	(setf (slot-value 'ynames) new))
  (slot-value 'ynames))

(defmeth gee-proto :link (&optional (new nil set))
"Message args: (&optional new)
Sets or returns link object."
  (when set
	(setf (slot-value 'link) new)
	(send self :needs-computing t))
  (slot-value 'link))

(defmeth gee-proto :varfun (&optional (new nil set))
"Message args: (&optional new)
Sets or returns variance function object."
  (when set
	(setf (slot-value 'varfun) new)
	(send self :needs-computing t))
  (slot-value 'varfun))

(defmeth gee-proto :corr-structure (&optional (new nil set))
"Message args: (&optional new)
Sets or returns correlation structure object."
  (when set
	(setf (slot-value 'corr-structure) new)
	(send self :needs-computing t))
  (slot-value 'corr-structure))

(defmeth gee-proto :offset (&optional (new nil set))
"Message args: (&optional (new nil set))
Sets or returns offset values."
  (when set
	(setf (slot-value 'offset) new)
	(send self :needs-computing t))
  (slot-value 'offset))

(defmeth gee-proto :pweights (&optional (new nil set))
"Message args: (&optional (new nil set))
Sets or returns prior weights."
  (when set
	(setf (slot-value 'pweights) new)
	(send self :needs-computing t))
  (slot-value 'pweights))

(defmeth gee-proto :group (&optional (new nil set))
"Message args: (&optional (new nil set))
Sets or returns group indicator."
  (when set
	(setf (slot-value 'group) new)
	(send self :needs-computing t))
  (slot-value 'group))

(defmeth gee-proto :scale (&optional (new nil set))
"Message args: (&optional (new nil set))
Sets or returns value of scale parameter."
  (if set 
      (setf (slot-value 'scale)  new))
  (slot-value 'scale))

(defmeth gee-proto :epsilon (&optional new)
"Message args: (&optional new)
Sets or returns tolerance for relative change in coefficients."
  (if new (setf (slot-value 'epsilon) new))
  (slot-value 'epsilon))

(defmeth gee-proto :epsilon-dev (&optional new)
"Message args: (&optional new)
Sets or returns tolerance for change in quasideviance."
  (if new (setf (slot-value 'epsilon-dev) new))
  (slot-value 'epsilon-dev))

(defmeth gee-proto :count-limit (&optional new)
"Message args: (&optional new)
Sets or returns maximum number of itrations."
  (if new (setf (slot-value 'count-limit) new))
  (slot-value 'count-limit))

(defmeth gee-proto :recycle (&optional (new nil set))
"Message args: (&optional new)
Sets or returns recycle option. If option is not NIL, current values
are used as initial values by :COMPUTE method."
  (when set
	(setf (slot-value 'recycle) new))
  (slot-value 'recycle))

(defmeth gee-proto :verbose (&optional (val nil set))
"Message args: (&optional (val nil set))
Sets or returns VERBOSE option. Iteration info is printed if option
is not NIL."
  (if set (setf (slot-value 'verbose) val))
  (slot-value 'verbose))

(defmeth gee-proto :formula (&optional (new nil set))
 "Message args: (&optional new)
Sets or returns model formula. "
  (when set
	(setf (slot-value 'formula) new))
  (slot-value 'formula))

(defmeth gee-proto :finite-sample-correction (&optional (new nil set))
 "Message args: (&optional new)
Sets or returns flag for finite-sample correction to variance. "
  (when set
	(setf (slot-value 'finite-sample-correction) new))
  (slot-value 'finite-sample-correction))

(defmeth gee-proto :eta ()
"Message args: ()
Returns linear predictor values for current fit."
    (slot-value 'eta))

(defmeth gee-proto :set-eta (&optional val)
  (if val
      (setf (slot-value 'eta) val)
      (setf (slot-value 'eta)
	    (+ (send self :offset) (matmult (send self :x-matrix) (send self :beta)))))
  (slot-value 'eta))

(defmeth gee-proto :beta ()
"Message args: ()
Returns coefficient values for current fit."
  (slot-value 'beta))

(defmeth gee-proto :coef-estimates () "Message args: ()
Returns coefficient values for current fit."
(slot-value 'beta))

(defmeth gee-proto :set-beta (val)
  "set coefficient estimates: only done from within compute-step-beta
   and :initialize-search"
  (setf (slot-value 'beta) val)
  (slot-value 'beta)
)

(defmeth gee-proto :covariance-matrix ()
 (if  (send self :needs-computing) (send self :compute))
  (slot-value 'covariance-matrix))

(defmeth gee-proto :set-covariance-matrix (val)
    (setf  (slot-value 'covariance-matrix) val))

(defmeth gee-proto :naive-covariance-matrix ()
 (if  (send self :needs-computing) (send self :compute))
  (slot-value 'naive-covariance-matrix))

(defmeth gee-proto :set-naive-covariance-matrix (val)
    (setf  (slot-value 'naive-covariance-matrix) val))

(defmeth gee-proto :quasideviance ()
"Message args: ()
Returns deviances for durrent fit."
  (slot-value 'quasideviance))

(defmeth gee-proto :set-deviances (val)
  (setf (slot-value 'quasideviance) val ))

(defmeth gee-proto :num-cases ()
  (length (first (column-list (send self :x)))))

(defmeth gee-proto :x (&optional x)
"Arguments: &optional x
Sets or returns predictor variables (without intercept). These may be a matrix, a single sequence or a list of sequences. See also :x-matrix"
  (if x
      (let ((x (cond
		((matrixp x) x)
		((vectorp x) (list x))
		((and (consp x) (numberp (car x))) (list x))
		(t x))))
	(setf (slot-value 'x) (if (matrixp x) x (apply #'bind-columns x))))
  (slot-value 'x) )
)

(defmeth gee-proto :conf-interval (&key transform (coverage 0.95) (names t))
"Arguments: &key transform (coverage 0.95) (names t)
Produces a matrix of Wald confidence intervals for beta with predictor
names if requested  and optionally transforms them to a different scale"
(let* ((beta (send self :beta))
       (se (sqrt (diagonal (send self :covariance-matrix))))
       (zalpha (normal-quant (+ coverage (/ (- 1 coverage) 2))))
       (baseci (apply #'bind-columns (list beta (- beta (* zalpha se)) (+ beta (* zalpha se)))))
       (xform (if (null transform) baseci (funcall transform baseci)))
       (name-list (if (send self :intercept) (cons "Intercept" (send self :predictor-names)) (send self :predictor-names)))
       )
  (if names (bind-columns name-list xform) xform)
))


(defmeth gee-proto :x-matrix ()
"Returns the design matrix (including intercept if present). See also :x"
  (cond ( (send self :intercept) (bind-columns (repeat 1 (length (send self :y))) (send self :x)))
	(t (send self :x))
  ))

(defmeth gee-proto :coef-standard-errors ()
(sqrt (diagonal (send self :covariance-matrix))))

(defmeth gee-proto :naive-coef-standard-errors ()
 (sqrt (diagonal (send self :naive-covariance-matrix))))

(defmeth gee-proto :predictor-names (&optional (names nil set))
"Message args: (&optional (names nil set))
With no argument returns the predictor names. NAMES sets the names."
  (if set (setf (slot-value 'predictor-names) (mapcar #'string names)))
  (let ((p (array-dimension (send self :x) 1))
	(p-names (slot-value 'predictor-names)))
    (if (not (and p-names (= (length p-names) p)))
	(setf (slot-value 'predictor-names)
	      (mapcar #'(lambda (a) (format nil "Variable ~a" a)) 
		      (iseq 0 (- p 1))))))
  (slot-value 'predictor-names))

(defmeth gee-proto :response-name (&optional (name "Y" set))
"Message args: (&optional name)
With no argument returns the response name. NAME sets the name."
  (if set (setf (slot-value 'response-name) (if name (string name) "Y")))
  (slot-value 'response-name))

(defmeth gee-proto :case-labels (&optional (labels nil set))
"Message args: (&optional labels)
With no argument returns the case-labels. LABELS sets the labels."
  (if set (setf (slot-value 'case-labels) 
		(if labels 
		    (mapcar #'string labels)
		    (mapcar #'(lambda (x) (format nil "~d" x)) 
			    (iseq 0 (- (send self :num-cases) 1))))))
  (slot-value 'case-labels))

;;;
;;;                   diagnostic methods
;;;

;;
;; quantities used in diagnostic calculations
;;
(defmeth gee-proto :scores (&optional (new nil set))
  (when set
	(setf (slot-value 'scores) new))
  (slot-value 'scores))

(defmeth gee-proto :working-residuals (&optional (new nil set))
  (when set
	(setf (slot-value 'working-residuals) new))
  (slot-value 'working-residuals))

(defmeth gee-proto :cluster-influences (&optional (new nil set))
  (when set
	(setf (slot-value 'cluster-influences) new))
  (slot-value 'cluster-influences))

(defmeth gee-proto :weight-blocks (&optional (new nil set))
  (when set
	(setf (slot-value 'weight-blocks) new))
  (slot-value 'weight-blocks))

;;
;; plot-linking methods
;; 

(defmeth gee-proto :bind-plot-to-gee (plot unit)
  (send plot :add-slot 'linked-gee self)
  (send plot :add-slot 'link-unit unit)
  (defmeth plot :links ()
    (let* ((link-list (send (slot-value 'linked-gee) :linked-plots (slot-value 'link-unit))))
      (if (member self link-list) link-list)
      ))
  (defmeth plot :linked (&optional (link nil set))
    (when set
	  (if link 
	      (send (slot-value 'linked-gee) :add-link (slot-value 'link-unit) self)
	    (send (slot-value 'linked-gee) :drop-link (slot-value 'link-unit) self)
	    )
	  ;(call-next-method link)
	  )
    (call-next-method)
    )
  (send plot :linked t)
)


(defmeth gee-proto :linked-plots (unit &optional (new nil set))
"Arguments: unit  &optional (new nil set)
Sets or returns the list of linked plots
unit is +group+ or +observation+"
(when set
      (if (= unit +group+) 
	  (setf (slot-value 'group-plots) new)
	(setf (slot-value 'observation-plots) new)
      ))
(if (= unit +group+) (slot-value 'group-plots) (slot-value 'observation-plots))
)


(defmeth gee-proto :drop-link (unit plot)
(if (= unit +group+)
    (setf (slot-value 'group-plots) (remove plot (slot-value 'group-plots)))
  (setf (slot-value 'observation-plots) (remove plot (slot-value 'observation-plots)))
  ))


(defmeth gee-proto :add-link (unit plot)
(if (= unit +group+)
    (setf (slot-value 'group-plots) (cons plot (slot-value 'group-plots)))
  (setf (slot-value 'observation-plots) (cons plot (slot-value 'observation-plots)))
  ))


;;;
;;; residuals
;;;

(defmeth gee-proto :raw-residuals () 
"Message args: ()
Returns the raw residuals for a model."
  (- (send self :y) (send self :fit-means)))

(defmeth gee-proto :pearson-residuals  ()
 "Returns Pearson residuals"
(send (send self :varfun) :pearson-resid (send self :y) (send self :fit-means))
)

(defmeth gee-proto :standardised-residuals ()
"Pearson residuals standardised by working correlation" 
(let* ((pres (send self :pearson-residuals))
       (times (send self :times))
       (groups (send self :group))
       (grouplist (unique groups))
       (timelist (unique times))
       (corr (send self :corr-structure))
       (Rinv-old (identity-matrix (length timelist)))
       (sqrtRinv (identity-matrix (length timelist)))
       (alpha (send self :alpha))
       (junk (dolist (g (coerce grouplist 'list) pres)
		     (let* (
			    (thisgroup (argselect groups g))
			    (ni (length thisgroup))
			    (thesetimes (if (null timelist) nil (select times thisgroup)))
			    (Rinv (send corr :rinv ni alpha thesetimes))
			    )
		       (if (null (equal Rinv Rinv-old) ) (setf sqrtRinv (first (chol-decomp Rinv))))
		       (setf (select pres thisgroup) (matmult sqrtRinv (select pres thisgroup)))
		       )))
       )
  pres)
)

(defmeth gee-proto :plot-standardised-residuals (&rest graphargs &key variable (smooth t) (show-labels t))
"Arguments: &rest graphargs &key variable (smooth t) (show-labels t)
Plot decorrelated Pearson residuals against eta or variable
with group hiliting mode available"
(cond ((compound-data-p variable) (mapcar #'(lambda (x) (send self :plot-standardised-residuals :variable x :smooth smooth :show-labels show-labels)) variable))
      (t
       (let* ((stdres (send self :standardised-residuals))
	      (xaxis (if (and variable (> variable 0)) (select (column-list (send self :x-matrix)) variable) (send self :eta)))
	      (xname (if (and variable (> variable 0)) (select (send self :predictor-names) ( - variable 1)) "Linear Predictor"))
	      (plot (if smooth (apply #'scatter-smooth (append (list xaxis stdres :variable-labels (list xname "Standardised Resid")) graphargs)) (apply #'plot-points (append (list xaxis stdres :variable-labels (list xname "Standardised Resid")) graphargs))))
	      )
	 (send plot :add-group-hilite (send self :group))
	 (send plot :mouse-mode 'group-hilite)
	 (if show-labels (send plot :showing-labels t))
	 (send self :bind-plot-to-gee plot +observation+)
	 plot
	 )
       )
      )
)

(defmeth gee-proto :plot-pearson-residuals (&rest graphargs &key variable (smooth t) (show-labels t))
"&rest graphargs &key variable (smooth t) (show-labels t)
Plot Pearson residuals against eta or variable
with group hiliting mode available"
(cond ((compound-data-p variable) (mapcar #'(lambda (x) (send self :plot-pearson-residuals :variable x :smooth smooth :show-labels show-labels)) variable))
      (t 
       (let* ((pres (send self :pearson-residuals))
	      (xaxis (if (and variable (> variable 0)) (select (column-list (send self :x-matrix)) variable) (send self :eta)))
	      (xname (if (and variable (> variable 0)) (select (send self :predictor-names) ( - variable 1)) "Linear Predictor"))
	      (plot (if smooth (apply #'scatter-smooth (append (list xaxis pres :variable-labels (list xname "Pearson Resid")) graphargs)) (apply #'plot-points (append (list xaxis pres :variable-labels (list xname "Pearson Resid")) graphargs))))
	      )
	 (send plot :add-group-hilite (send self :group))
	 (send plot :mouse-mode 'group-hilite)
	 (if show-labels (send plot :showing-labels t))
	 (send self :bind-plot-to-gee plot +observation+)
	 plot
	 )
       )
      )
)

(defmeth gee-proto :deviance-residuals ()
"Independence model deviance residuals"
(let* ((y (send self :y))
       (mu (send self :fit-means))
       (varfun (send self :varfun))
       )
  (mapcar #'(lambda (a b) (* (sign (- a b)) (sqrt (send varfun :quasi-deviance a b)))) y mu)
))


(defmeth gee-proto :plot-deviance-residuals (&rest graphargs &key variable (smooth t) (show-labels t))
"&rest graphargs &key variable (smooth t) (show-labels t)
Plot deviance residuals against eta or variable
with group hiliting mode available"
(cond ((compound-data-p variable) (mapcar #'(lambda (x) (send self :plot-deviance-residuals :variable x :smooth smooth :show-labels show-labels)) variable))
      (t 
       (let* ((res (send self :deviance-residuals))
	      (xaxis (if (and variable (> variable 0)) (select (column-list (send self :x-matrix)) variable) (send self :eta)))
	      (xname (if (and variable (> variable 0)) (select (send self :predictor-names) ( - variable 1)) "Linear Predictor"))
	      (plot (if smooth (apply #'scatter-smooth (append (list xaxis res :variable-labels (list xname "Deviance Resid")) graphargs)) (apply #'plot-points (append (list xaxis res :variable-labels (list xname "Deviance Resid")) graphargs))))
	      )
	 (send plot :add-group-hilite (send self :group))
	 (send plot :mouse-mode 'group-hilite)
	 (if show-labels (send plot :showing-labels t))
	  (send self :bind-plot-to-gee plot +observation+)
	plot
	 )
       )
      )
)


(defmeth gee-proto :group-names ()
(mapcar #'(lambda (xx) (write-to-string xx :escape nil)) (unique (send self :group)))
)


;;
;;  deletion diagnostics: 
;;  Preisser & Qaqish (Biometrika 1996)
;;

(defconstant +not-group+ 0)
(defconstant +group+ 1)
(defconstant +observation+ 0)
(defconstant +cluster+ 1)

(defmeth gee-proto :leverage (&key (unit +observation+))
"Leverage values for group or observation, from Preisser & Qaqish."
(let* ((Hi (send self :cluster-influences)) 
      )
  (cond ((null Hi)                      ;;was it  computed earlier ?
	 (when (send self :verbose) (format t "~&computing diagnostic information~%"))
	 (send self :compute-step-beta t) 
	 (send self :leverage :unit unit)
	 )
	((= unit +group+) (mapcar #'(lambda (x) (list (first x) (matrix-trace (third x)))) Hi)
	 )
	((= unit +observation+)
	 (let* (
		(rankij (argrank (apply #'append (mapcar #'second Hi))))
		(hij (apply #'append  (mapcar #'(lambda (x) (diagonal (third x))) Hi)))
		)
	   (select hij rankij)
	   )
	 )
	(t (error "Can't happen in :leverage"))
	)
  )
)
       
(defmeth gee-proto :plot-leverage (&rest graphargs &key (unit +observation+) (show-labels t))
"Plot group or individual leverages with link to other plots, and group 
hiliting for individual leverages plot"
(cond ((= unit +group+)
       (let* ((infl (send self :leverage :unit +group+))
	      (infseq (mapcar #'second infl))
	      (rankseq (argrank (mapcar #'first infl)))
	      (rankg (argrank (unique (send self :group))))
	      (lbls (select (send self :group-names) rankg))
	      (plot (apply #'plot-points (append (list (iseq 1 (length infseq)) (select infseq rankseq) :variable-labels (list "Index" "Leverage") :point-labels lbls) graphargs)))
	      )
	 (if show-labels (send plot :showing-labels t))
	 (send self :bind-plot-to-gee plot +group+)
	 plot
	 )
       )
      ((= unit +observation+)
       (let* ((res (send self :leverage :unit +observation+))
	      (plot (apply #'index-plot (append (list  res :variable-labels (list "Index" "Leverage") :show-labels show-labels) graphargs)))
	     )
	 (send plot :add-group-hilite (send self :group))
	 (send plot :mouse-mode 'group-hilite)
	 (if show-labels (send plot :showing-labels t))
	 (send self :bind-plot-to-gee plot +observation+)
	 plot
	 )
       )
      (t 
       (error "invalid value of unit in :plot-leverage")
       )
      )
)


(defmeth gee-proto :dbeta (&key (unit +observation+))
"Deletion diagnostics DBETAC, DBETAO of Preisser & Qaqish"
(let* ((H (send self :cluster-influences)) ;;was it  computed earlier ?
      )
  (cond ((null H)	
	 (when (send self :verbose) (format t "~&computing diagnostic information~%"))
	 (send self :compute-step-beta t) 
	 (send self :dbeta :unit unit)
	 )
	((= unit +group+)
	 (let* (
		(W (send self :weight-blocks))
		(groups (send self :group))
		(grouplist (unique groups))
		(E (send self :working-residuals))
		(XWXinv  (send self :naive-covariance-matrix))
		(p (first (array-dimensions XWXinv)))
		(p-1 (- p 1))
		(dbeta nil)
		)
	   (dolist (g grouplist dbeta)
		   (let* ((thisgroup (argselect groups g))
			  (Xi (select (send self :x-matrix) thisgroup (iseq p)))
			  (Ei (second (assoc g E :test #'equal)))
			  (Wi (third (assoc g W :test #'equal)))
			  (Hi (third (assoc g H :test #'equal)))
			  (dbetai (- (matmult XWXinv (transpose Xi) Wi (- (identity-matrix (length thisgroup)) Hi) Ei)))
			  )
		     (setf dbeta (cons (list g dbetai) dbeta))
		     )
		   )
	   dbeta
	   )
	 )
	((= unit +observation+)
	  (let* (
		(W (send self :weight-blocks))
		(groups (send self :group))
		(grouplist (unique groups))
		(E (send self :working-residuals))
		(XWXinv (send self :naive-covariance-matrix))
		(p (first (array-dimensions XWXinv)))
		(ndbeta nil)
		(junk (dolist (g grouplist ndbeta)
			      (let* ((thisgroup (argselect groups g))
				     (Xi (select (send self :x-matrix) thisgroup (iseq p)))
				     (Ei (second (assoc g E :test #'equal)))
				     (Wi (third (assoc g W :test #'equal)))
				     (Hi (third  (assoc g H :test #'equal)))
				     (Vi (inverse Wi))
				     (ni (length thisgroup))
				     )
				(dolist (j (iseq ni))
					(let* (
					       (notj (remove j (iseq ni)))
					       (Xij (select Xi j (iseq p)))
					       (Xi-j (select Xi notj (iseq p)))
					       (Eij (select Ei j))
					       (Ei-j (select Ei notj))
					       (Wij (select Wi j j))
					       (Vi-j (select Vi notj notj))
					       (Vji-j (select Vi j notj))
					       (VVinv (matmult Vji-j (inverse Vi-j)))
					       (tXij (- Xij (matmult VVinv Xi-j)))
					       (tEij (- Eij (matmult VVinv Ei-j)))
					       (tQij (select (matmult tXij XWXinv (transpose tXij)) 0  0))
					       (db (- (coerce (first (row-list (apply #'* Xij (/ tEij (- (/ Wij) tQij))))) 'list)))
					       )
					  (setf ndbeta (cons (append (select thisgroup (list j)) db) ndbeta))
					  )
					)
				)
			      ))
		(coldbeta (column-list (apply #'bind-rows ndbeta)))
		(dbeta (row-list  (matmult  (apply #'bind-columns (cdr coldbeta)) XWXinv  )))
		(rankij (argrank (mapcar #'first ndbeta)))
		)
	    (select dbeta rankij)
	    )
	  )
	(t (error "Can't happen in :dbeta"))
	)
  )
)
	   
(defmeth gee-proto :plot-dbeta (&rest graphargs  &key (unit +observation+) variable (show-labels t))
"Plot the change in coefficients for specified variables on deletion of each observation or group. Default is all variables"
 (cond  ((= unit +group+)
	 (let* ((dbeta (send self :dbeta :unit +group+))
		(rankdb (argrank (mapcar #'first dbeta)))
		(db (column-list (apply #'bind-rows (select (mapcar #'second dbeta) rankdb))))
		(rankg (argrank (unique (send self :group))))
		(gnames (select (send self :group-names) rankg))
		(varnames (if (send self :intercept) (cons "intercept" (send self :predictor-names)) (send self :predictor-names)))
		(indices (if (null variable) (iseq  (length db))  variable))
		(plots (mapcar #'(lambda (i) (apply #'index-plot (append (list (elt db i) :point-labels gnames :title (elt varnames i) :variable-labels (list "Index" "Cluster DBETA")) graphargs)))   indices))
		)
	   (if show-labels (mapcar #'(lambda (x) (send x :showing-labels t)) plots))
	   (mapcar #'(lambda (x) (send self :bind-plot-to-gee x +group+)) plots)
	   plots
	   ))
	
	((= unit +observation+)
	 (let* ((db (column-list (apply #'bind-rows (send self :dbeta :unit +observation+))))
		(indices (if (null variable) (iseq  (length db))  variable))
		(varnames (if (send self :intercept) (cons "intercept" (send self :predictor-names)) (send self :predictor-names)))
		(plots (mapcar #'(lambda (i) (apply #'index-plot (append (list (elt db i)  :title (elt varnames i) :variable-labels (list "Index" "Observation DBETA")) graphargs)))   indices))
		)
	   (mapcar #'(lambda (x)
		       (send x :add-group-hilite (send self :group)) 
		       (send x :mouse-mode 'group-hilite)
		       )
		   plots)
	   (if show-labels (mapcar #'(lambda (x) (send x :showing-labels t)) plots))
	   (mapcar #'(lambda (x) (send self :bind-plot-to-gee x +observation+)) plots)
	plots
	   ))
	(t 
	 (error "Invalid value of unit in :plot-dbeta.")
	 )
	)
)


(defmeth gee-proto :cooks-distance (&key (unit +observation+))
"Generalisation of Cook's distance to GEEs"
(let* ((dbeta (send self :dbeta :unit unit))
       (scale (send self :scale))
       (p (length (send self :beta)))
       (XWX (inverse (send self :naive-covariance-matrix)))
       )
  (cond ((= unit +observation+)
	 (/ (mapcar #'(lambda (dbi) (inner-product dbi (matmult  XWX dbi))) dbeta) (* p scale))
	 )
	((= unit +group+)
	  (mapcar #'(lambda (dbi) (list (first dbi) (/ (inner-product (second dbi)(matmult  XWX (second dbi))) (* p scale)))) dbeta)
	 )
	(t (error "~&Invalid value of unit in :cook-distance~%"))
	)
  )
)


(defmeth gee-proto :plot-cooks-distance (&rest graphargs &key (unit +observation+) (show-labels t))
"Plot group or individual Cook's distance with link to other plots, and group 
hiliting for individual leverages plot"
(cond ((= unit +group+)
       (let* ((infl (send self :cooks-distance :unit +group+))
	      (infseq (mapcar #'second infl))
	      (rankseq (argrank (mapcar #'first infl)))
	      (rankg (argrank (unique (send self :group))))
	      (lbls (select (send self :group-names) rankg))
	      (plot (apply #'plot-points (append (list (iseq 1 (length infseq)) (select infseq rankseq) :variable-labels (list "Index" "Cook's distance") :point-labels lbls) graphargs)))
	      )
	 (if show-labels (send plot :showing-labels t))
	 (send self :bind-plot-to-gee plot +group+)
	 plot
	 )
       )
      ((= unit +observation+)
       (let* ((res (send self :cooks-distance :unit +observation+))
	      (plot (apply #'index-plot (append (list  res :variable-labels (list "Index" "Cook's distance") :show-labels show-labels) graphargs)))
	     )
	 (send plot :add-group-hilite (send self :group))
	 (send plot :mouse-mode 'group-hilite)
	 (if show-labels (send plot :showing-labels t))
	 (send self :bind-plot-to-gee plot +observation+) 
	 plot
	 )
       )
      (t 
       (error "invalid value of unit in :plot-cooks-distance")
       )
      )
)

;;
;; initialisation methods
;;

(defmeth gee-proto :domain-error (eta)
"Checks whether we have invalid values of eta"
(let* ((link (send self :link))
       (varfun (send self :varfun))
       )
  (null (if (send link :valideta eta) (send varfun :validmu (send link :means eta)) nil))))


(defmeth gee-proto :initialize-search ()
  (let* (
	 (p (second (array-dimensions (send self :x-matrix))))
	 )
  (send self :alpha nil)
  (cond ( (and (send self :beta) (null (send self :eta))) ) 
	( (null (send self :domain-error 0)) (send self :set-beta (repeat 0 p)))
	( (send self :intercept) (send self :set-beta (append (list 0.5) (repeat 0 (- p 1)))))
	(t (let* ((y (send self :y))
		  (x (column-list (send self :x-matrix)))
		  (ybar (mean y))
		  (link (send self :link))
		  (etabar (if (send link :validmu ybar) (send link :eta ybar) (error "I need help. Please specify starting values using :init-beta~%")))
		  (beta0 (send self :set-beta (/ etabar (* p (mapcar #'mean x)))))
		  (eta0 (send self :set-eta))
		  )
	     (if (send self :domain-error eta0) 
		 (error "I need help. Please specify starting values using :init-beta~%") )
	     ))
	)
  )
  (send self :set-eta)
)

(defmeth gee-proto :fit-means (&optional (eta (send self :eta)))
"Message args: (&optional (eta (send self :eta)))
Retruns mean values for current or supplied ETA."
  (send (send self :link) :means eta))

;;;
;;; :COMPUTE method -- does all the work
;;;

(defmeth gee-proto :compute ()
"Refits the model"
(let* ((epsilon (send self :epsilon))
       (epsilon-dev (send self :epsilon-dev))
       (maxcount (send self :count-limit))
       (low-lim (* 2 (/ machine-epsilon epsilon)))
       (verbose (send self :verbose)))
  (unless (and (send self :eta) (send self :recycle))
	  (send self :initialize-search))
  (do ((count 1 (+ count 1))
       (beta 0 (send self :beta))
       (last-beta -1 beta)
       (dev  0 (send self :quasideviance))
       (last-dev  -1 dev))
	((or (> count maxcount) 
	     (< (max (abs (/ (- beta last-beta)
			     (pmax (abs last-beta) low-lim))))
		epsilon)
	     (< (abs (- dev last-dev)) epsilon-dev)) (if (> count maxcount) (format t "~%~%Failed to converge in ~d iterations~%Treat results with caution ~%" maxcount)))
	(send self :compute-step-alpha)
	(send self :compute-step-beta)
	(if verbose 
	    (format t "Iteration ~d: quasideviance = ~,6g~%" 
		    count (send self :quasideviance)))
	)
  )
(send self :needs-computing nil)
)


(defmeth gee-proto :compute-step-beta (&optional (save-scores *gee-diagnostics-while-fitting*))
"Does one step of quasi-scoring: called by :compute
  updates beta, quasideviance  and covariance matrix "
  (let* ((y (send self :y))
	 (eta (send self :eta))
	 (mu (send self :fit-means eta))
	 (pw (send self :pweights))
	 (grouplist (unique (send self :group)))
	 (timelist (unique (send self :times)))
	 (ntimes (length timelist))
	 (alpha (send self :alpha))
	 (phi (/ (send self :scale)))
	 (beta (send self :beta)) 
	 (p (length beta))
	 (x (row-list (send self :x-matrix)))
	 (cov (zero-matrix p))
         (scores nil)
	 (working-residuals nil)
	 (Hi nil)
	 (Wi-gp nil)
	 (DVS (repeat 0 p))
	 (DVD (zero-matrix p))
	 (junk (dolist (g (coerce grouplist 'list))
	    (let* (
		(thisgroup (argselect (send self :group) g))
		(thesetimes (if (null timelist) nil (select (send self :times) thisgroup)))
		(ni (length thisgroup))
		(Rinv  (send (send self :corr-structure) :rinv ni alpha thesetimes))
		(mui  (select mu thisgroup))
		(Si  (- (select y thisgroup) mui))
		(sqinvAi  (diagonal (/ (sqrt (* (select pw thisgroup)  (send (send self :varfun) :varweight mui))))))
		(derivi  (/ (send (send self :link) :derivs mui)))
		(Vinvi  (* phi (matmult sqinvAi Rinv sqinvAi)))
		(Xi (apply #'bind-rows (select x thisgroup)))
		(Di  (matmult (diagonal derivi) Xi))
		(DVi (matmult (transpose Di)  Vinvi))
		(DVSi (matmult DVi Si))
		(DVDi (matmult DVi Di))
		 )
		(setf cov  (+ cov (outer-product DVSi DVSi)))
		(setf DVS  (+ DVS DVSi))
		(setf DVD  (+ DVD DVDi))
		(when save-scores 
		      ; compute all sorts of useful diagnostic
		      ; components only as needed.
		      (setf scores (cons (cons g (list DVSi)) scores))
		      (setf working-residuals (cons (cons g (list (* (/ derivi) Si))) working-residuals))
		      (setf Hi  (cons (list  g thisgroup  (matmult Xi (send self :naive-covariance-matrix) (transpose Di) Vinvi (diagonal derivi))) Hi))
		      (setf Wi-gp (cons (list g thisgroup  (matmult (diagonal derivi) Vinvi (diagonal derivi))) Wi-gp))
		      )
		)))
	 (Linv-DVD (inverse (first (chol-decomp DVD))))
	 (DVDinv (matmult (transpose Linv-DVD) Linv-DVD))
	 ;; new value of beta -- needs to be tested for domain errors
	 (newbeta (+ beta (matmult DVDinv DVS )))
	 (neweta  (+ (send self :offset) (matmult (send self :x-matrix) newbeta)))	 
	 )
    ;; try new value of beta, halve step size if it is out of range.
    ;;
    (send self :set-beta 
	  (do* ( 
		(beta1 newbeta (/ (+ beta beta1) 2))
		(eta1 neweta  (+ (send self :offset) (matmult (send self :x-matrix) beta1)))
		)
	       ( (null (send self :domain-error eta1)) beta1)
	        (format t "Step size reduced ~%")
		(if (send self :verbose) (format t "beta: ~a~%" beta1))
	       )
	  )
    ;; update everything else
    (send self :set-naive-covariance-matrix DVDinv)
    (send self :set-covariance-matrix (* (if (send self :finite-sample-correction) (/ (- 1 (/ (length grouplist)))) 1 )  (matmult  DVDinv cov DVDinv)))
    (send self :set-eta)
    (send self :scores scores)
    (send self :cluster-influences Hi)
    (send self :working-residuals working-residuals)
    (send self :weight-blocks Wi-gp)
    (send self :set-deviances (send (send self :varfun) :quasi-deviance y mu))
  
 )
)



(defmeth gee-proto :compute-step-alpha ()
"Updates the estimates of the correlation parameters
 Called by :compute"
(let* ( (ri (* (send self :pweights) (send (send self :varfun) :pearson-resid (send self :y) (send self :fit-means))))
	
	(corrthings (send (send self :corr-structure) :estimate ri (send self :group) (send self :times) (send self :allow-missing) ))
	)
  (send self :alpha corrthings)
  (send self :scale (/ (second corrthings)))
))


;;;
;;; DISPLAY method 
;;;

(defmeth gee-proto :display ()
"Prints the model summary"
   (if (send self :needs-computing) (send self :compute))
  (let ((coefs (coerce (send self :coef-estimates) 'list))
	(se-s (send self :coef-standard-errors))
	(n-se-s (send self :naive-coef-standard-errors))
	(x (send self :x))
	(p-names (send self :predictor-names)))
    (format t "~%GEE Estimates:~2%")
    (format t "~25tCoefficient~40t Std Error~60t Naive Std Error~%")
    (when (send self :intercept)
	  (format t "Constant~25t~13,6g~40t(~,6g)~60t(~,6g)~%" (car coefs) (car se-s) (car n-se-s))
	  (setf coefs (cdr coefs))
	  (setf n-se-s (cdr n-se-s))
	  (setf se-s (cdr se-s)))
    (dotimes (i (array-dimension x 1)) 
	     (format t "~a~25t~13,6g~40t(~,6g)~60t(~,6g)~%"
		     (select p-names i) (car coefs) (car se-s) (car n-se-s))
	     (setf coefs (cdr coefs) se-s (cdr se-s) n-se-s (cdr n-se-s)))
    (format t "~%")
    (format t "Scale Estimate:~25t~13,6g~%" (send self :scale))
    (format t "Independence model deviance:~25t~13,6g~%" (send self :quasideviance))
    (format t "Number of cases:~25t~9d~%" (send self :num-cases))
    (format t "Link:~25t") (send (send self :link) :print)
    (send (send self :varfun) :prettyprint)
    (send (send self :corr-structure) :prettyprint (send self :alpha))
 ))
 
(defmeth gee-proto :display-with-formula (&key (block-only *gee-display-block-only*))
"Arguments: &key (block-only *gee-display-block-only*)
Prints the model summary based on a model formula"
(cond ( (null (send self :formula)) (error "No formula present"))
      (t 
       (if (send self :needs-computing) (send self :compute))
       (let* ((beta (coerce (send self :coef-estimates) 'list))
	     (formula (send self :formula))
	     (cov-mat (send self :covariance-matrix))
	     (xintercept (if (send self :intercept) 1 0))
	     (nblocks (length (send formula :block-indices)))
	     (block-indices (+ xintercept (send formula :block-indices)))
	     (varnames (if (send self :intercept) (cons "Intercept" (send self :predictor-names)) (send self :predictor-names)))
	     )
	 (format t "~%GEE Estimates:~2%")
	 (format t "~a~20t~a~35t~a~%" "Block" "Wald Chisq" "p-value")
	 (block-test (list 0) beta cov-mat :names varnames :blockname "Intercept" :block-only block-only) 
;;	 (dolist (i (iseq (- nblocks 1) 0))
	 (dolist (i (iseq 0 (- nblocks 1)))
		 (block-test (elt block-indices i) beta cov-mat :names varnames :blockname  (elt (send formula :block-names) i) :block-only block-only))
	 (format t "~%")
	 (format t "Scale Estimate:~25t~13,6g~%" (send self :scale))
	 (format t "Independence model deviance:~25t~13,6g~%" (send self :quasideviance))
	 (format t "Number of cases:~25t~9d~%" (send self :num-cases))
	 (format t "Link:~25t") (send (send self :link) :print)
	 (send (send self :varfun) :prettyprint)
	 (send (send self :corr-structure) :prettyprint (send self :alpha))
	 )
       )
))

;;;;
;;;; :ISNEW method
;;;;

(defmeth gee-proto :isnew (&key x 
				 y
				 group
				 (times nil)
				 link
				 varfun
				 (corr-structure *gee-default-correlation*)
				 (offset 0)
				 (intercept t)
				 (formula nil)
				 pweights
				 (allow-missing nil)
				 (verbose t)
				 predictor-names
				 response-name
				 (recycle nil)
				 case-labels
				 init-beta
				 init-alpha
				 count-limit
				 (print-summary nil)
				 (print t)
				 (finite-sample-correction *gee-correction*))
  (send self :x x)
  (send self :y y)
  (send self :group (coerce group 'list))
  (send self :varfun varfun)
  (send self :corr-structure corr-structure)
  (send self :link link)
  (send self :offset offset)
  (send self :times times)
  (send self :intercept intercept)
  (if pweights (send self :pweights pweights) (send self :pweights (repeat 1 (length y))))
  (send self :allow-missing allow-missing)
  (send self :formula formula)
  (send self :recycle recycle)
  (send self :verbose verbose)
  (send self :set-beta init-beta)
  (send self :finite-sample-correction finite-sample-correction)
;;  (if included (send self :included included))
  (if predictor-names (send self :predictor-names predictor-names))
  (if response-name (send self :response-name response-name))
  (if (or y case-labels) (send self :case-labels case-labels)) ; needs fixing
(send self :epsilon *gee-tolerance*) 
(send self :epsilon-dev  *gee-tolerance*)
(if (null count-limit) (send self :count-limit *gee-count-limit*) (send self :count-limit count-limit))
 (if print (send self :display))
 (if print-summary (send self :display-with-formula))
)

;;;;;
;;;;; Nice friendly function to construct these models
;;;;;


(defun gee-model (&key x y g (times nil)
		             (error nil error-given) 
		             (link nil link-given)
			     (correlation *gee-default-correlation* corr-given)
			     (offset 0)
			     (intercept t) 
			     (verbose *gee-verbose*)
			     (print nil has-print)
			     (print-summary t)
			     (allow-missing *gee-allow-missing*)
			     (init-beta nil)
			     (init-alpha nil)
			     (count-limit *gee-count-limit*)
			     (prior-weights nil)
			     (predictor-names nil has-names)
			     (finite-sample *gee-correction*)
			     )
  (if (null (and x y g (or link-given error-given))) (error "You must supply :x :y :g and at least one of :link and :error~%"))
  (if (and verbose (null corr-given)) (format t "No correlation specified, using *gee-default-correlation*~%"))
  (if (and verbose (null link-given)) (format t "No link specified - canonical link used~%"))
  (if (and verbose (null error-given)) (format t "No error function specified - using default for this link~%"))
  (cond ((objectp x) 
	 (send gee-proto :new :x (send x :design-matrix) :y y :group (coerce g 'list)
	       :times times :corr-structure correlation 
	       :link (if link-given link  (eval (cdr (assoc error default-link)))) 
	       :varfun (if error-given error (eval (cdr (assoc link default-error))))
	       :intercept intercept :verbose verbose  
	       :offset offset 
	       :allow-missing allow-missing
	       :init-beta init-beta
	       :count-limit count-limit
	       :pweights prior-weights
	       :predictor-names (if has-names predictor-names (send x :name-list))
	       :formula x
	       :print print
	       :print-summary print-summary
	       ))
	(t  (send gee-proto :new :x x :y y :group (coerce g 'list) 
	       :times times :corr-structure correlation 
	       :link (if link-given link  (eval (cdr (assoc error default-link)))) 
	       :varfun (if error-given error (eval (cdr (assoc link default-error))))
	       :intercept intercept :verbose verbose  
	       :offset offset 
	       :allow-missing allow-missing
	       :init-beta init-beta
	       :count-limit count-limit
	       :pweights prior-weights
	       :predictor-names predictor-names
	       :print (if has-print print print-summary)
	       :finite-sample-correction finite-sample
	       )
	    )
	)
)

;;;
;;; control variables. They need to be after the definition of the correlation
;;; structures
;;;

;default value for :verbose
(defvar *gee-verbose* t)

;default value for :count-limit (maximum number of iterations)
(defvar *gee-count-limit* 30)

;t to give only Wald tests for blocks by default
;nil to display all parameter estimates
(defvar *gee-display-block-only* nil)

;default value of :allow-missing
(defvar *gee-allow-missing* nil)

;default working correlation 
;(eg for compatability with SPIDA which uses exchangeable-corr
; as the default)
(defvar *gee-default-correlation* independence-corr)

; default convergence tolerance
(defvar *gee-tolerance* (sqrt (sqrt machine-epsilon)))

; save diagnostic information while fitting the model
; (slower fitting, but faster diagnostics, if t)
(defvar *gee-diagnostics-while-fitting* nil)

;number of digits for printing out correlations
(defvar *gee-corr-digits* 2)


; use finite-sample correction to variance estimate
; (recommended -- only nil for compatibility)
(defvar *gee-correction* nil)