(load "expsurv1.lsp")
(defun term (x &key namebase)
  "A metric predictor variable" 
  (if namebase 
      (list (bind-columns x)  namebase))
  (bind-columns x)
  )

(defun factor (x &key namebase ) 
"Defines treatment contrast matrix for x and optionally 
a set of names based on namebase"
  (let* (
         (xlist (sort-data (remove-duplicates (coerce x 'list) :test #'equalp)))
         (p (- (length xlist) 1))
         (rows (row-list (select (identity-matrix (+ p 1)) 
                                 (iseq 0 p) (iseq 1 p))))
         (decoder (mapcar #'cons xlist rows))
         (xmatrix (apply #'bind-rows 
                         (coerce (map-elements #'(lambda (xx) 
                                     (cdr (assoc xx decoder
                                                 :test #'equalp))) 
                                 x) 'list)))
         (xnames (if (null  namebase) 
                     nil 
                   (cdr (mapcar #'(lambda (xx) 
                                    (write-to-string (list namebase xx)
                                                     :escape nil)) 
                                xlist))))
    )
    (if (null namebase) 
        xmatrix 
      (list xmatrix xnames))
    )
)

(defun interaction ( xs &key (is-factor (repeat t (length xs))) namebases ) 
"Interactions of any order for categorical and continuous variables"
  (let*  ( (n (length (first xs)))
          (zlist nil)
          (znames nil)
          (p (length xs))
          (names (if (null namebases) 
                     (repeat "." p) 
                   namebases))
          (junk (dotimes (i p)
                        (let* (
                               (thisz  (if (elt is-factor i) 
                                           (factor (elt xs i) 
                                                   :namebase (elt names i)) 
                                         (list (bind-columns (elt xs i)) 
                                               (list (elt names i)))))
                               )
                          (if (null zlist) 
                              (setf zlist (column-list (first thisz))) 
                            (setf zlist 
                                  (apply #'append 
                                         (mapcar #'(lambda (xx) 
                                                     (mapcar #'(lambda (yy) 
                                                                 (*  xx yy)) 
                                                             (column-list 
                                                              (first thisz))))
                                                 zlist))))
                          (if (null znames) 
                              (setf znames (second thisz))  
                            (setf znames 
                                  (apply #'append 
                                         (mapcar #'(lambda (xx) 
                                                     (mapcar #'(lambda (yy) 
                                                                 (write-to-string 
                                                                  (list xx yy) 
                                                                  :escape nil)) 
                                                             (second thisz))) 
                                                 znames))))
                          )
                        )
                )
          (zmat (apply #'bind-columns zlist))
          )
    (if (null namebases) 
        zmat 
      (list zmat znames))
))
          
(defun design (varlist) (apply #'bind-columns (mapcar #'first varlist)))

(defun names (varlist) (apply #'append (mapcar #'second varlist)))

(defun col-count (matlist) 
  "adds up columns"  
  (apply #'sum (mapcar #'(lambda (x) 
                           (cond ((matrixp x) 
                                  (second (array-dimensions x))) 
                                 (t 1))) 
                       matlist)))

(defun formula (model-formula) 
"Args: model-formula
Calculates a design matrix, variable names and information for running 
(block-test) on the fitted model. The argument is  a quoted list of metric 
variables, factor variables and interactions. Metric variables are 
specified by (term x), factor variables by (factor a) and interactions
 by (interaction (list a b x) :is-factor (list t t nil)). "
(let* (
       (design-matrix nil)
       (name-list nil)
       (block-names nil)
       (block-indices nil)
       (i 0)
       (junk (dolist (x model-formula) 
                     (setf i (if (null design-matrix) 
                                 0 
                               (col-count design-matrix)))
                     (cond
                      ((eql (car x) 'term)  
                       (let* ((namepos (position ':namebase x))
                              )
                         (setf 
                          design-matrix (cons (eval (second x)) design-matrix)
                          name-list (cons (if (null (find ':namebase x)) 
                                              (list (write-to-string (second x))) 
                                            (list (fourth x))) name-list)
                          block-names (cons (write-to-string 
                                             (if (null namepos) 
                                                 (second x) 
                                               (elt x (+ 1 namepos)))) 
                                            block-names)
                          block-indices (cons (list i) block-indices)
                          )
                         )
                       )
                       ((eql (car x) 'factor)
                        (let* ((namepos (position ':namebase x))
                               (factorx (eval 
                                         (if (null namepos) 
                                             (append x 
                                                     (list :namebase 
                                                           (write-to-string (second x)))) 
                                           x)))
                               )
                          (setf 
                           design-matrix (cons (first factorx) design-matrix)
                           name-list (cons (second factorx) name-list)
                           block-names (cons (write-to-string 
                                              (if (null namepos)
                                                  (second x) 
                                                (elt x (+ 1 namepos)))) 
                                             block-names)
                           block-indices (cons 
                                          (iseq i (- (col-count design-matrix) 1)) 
                                          block-indices)
                           )
                          )
                        )
                       ((eql (car x) 'interaction)
                        (let* ( 
                               (namepos (position ':namebases x))
                               (interx (eval 
                                        (if (null namepos)
                                            (append x 
                                                    (list :namebases 
                                                          (cons 'list 
                                                                (mapcar 
                                                                 #'write-to-string 
                                                                 (cdr (second x))))))
                                          x)))
                               )
                          (setf 
                           design-matrix (cons (first interx) design-matrix)
                           name-list (cons (second interx) name-list)
                           block-names (cons (write-to-string 
                                              (cons 'interaction 
                                                    (if (null namepos) 
                                                        (cdr (second x)) 
                                                      (cdr (elt x (+ 1 namepos)))) 
                                                    ) 
                                              :escape nil)  
                                             block-names)
                           block-indices (cons 
                                          (iseq i (- (col-count design-matrix) 1))
                                          block-indices)
                           )
                          )
                        )
                       (t (error "I don't understand ~s" x))
                       )
                     )
             )
       )
  (list (apply #'bind-columns design-matrix) 
        (apply #'append name-list) 
        block-names 
        (- (- (col-count design-matrix) 1)
           block-indices))
  )
)

(defun block-test (index beta covmat &key blockname names (block-only nil))
  (let* (
         (subbeta (select beta index))
         (subcov  (select covmat index index))
         (waldchisq (inner-product subbeta (matmult (inverse subcov) subbeta)))
         (waldp (- 1 (chisq-cdf waldchisq (length index))))
         (subse (sqrt (diagonal subcov)))
         (subz (/ subbeta subse))
         (subp (* 2 (- 1 (normal-cdf (abs subz)))))
         (nn (if (null names) 
                 (repeat "" (length index)) 
               (select names index)))
         (blockn (if (null blockname) 
                     "block" 
                   blockname))
         )
    (format t "~a~20t~13,5g~35t~,4f~%" blockn waldchisq waldp)
    (unless block-only
            (format t "~5t Variable~25t Estimate~40t Std.Err.~%")
            (dolist (i (iseq 0 (- (length index) 1)))
                    (format t "~5t~a~25t~13,5g~40t(~,6g)~%" 
                            (select nn i) 
                            (select subbeta i) 
                            (select subse i) 
                            ))
            )
    )
  )     

(defun formula-display 
  (beta cov-mat model-formula &key (block-only t) (intercept t)) 
"Display a model summary for a model with coefficient estimates 
beta and covariance estimates cov-mat and design matrix from 
model-formula.  This should really be a method for the relevant 
object. Little checking is performed: garbage in; garbage out"
(let* ( 
       (xintercept (if intercept 1 0))
       (nblocks (length (fourth model-formula)))
       (varnames (if intercept (cons "Intercept" (second model-formula)) 
                   (second model-formula)))
       (block-indices (+ xintercept (fourth model-formula)))
       )
  (format t "~a~20t~a~35t~a~%" "Block" "Wald Chisq" "p-value")
  (if (equal (length beta) 
             (+ xintercept (second (array-dimensions (first model-formula)))))
      (dolist (i (iseq (- nblocks 1) 0))
              (block-test (elt block-indices i) 
                          beta 
                          cov-mat 
                          :names varnames 
                          :blockname (elt (third model-formula) i) 
                          :block-only block-only))
    (error "length of beta doesn't match design matrix")
    )
  )
)

(defproto model-formula-proto 
  '(formula design-matrix name-list block-names block-indices) 
  '() 
  *object* 
  "Model formula")

(defmeth model-formula-proto :formula (&optional (new nil set))
"Set or return model formula"
  (when set
        (setf (slot-value 'formula) new)
        (send self :compute))
  (slot-value 'formula))

(defmeth model-formula-proto :isnew (&key formula)
  (send self :formula formula)
  )

(defmeth model-formula-proto :design-matrix (&optional (new nil set))
"Sets or returns design-matrix"
  (when set
        (setf (slot-value 'design-matrix) new))
  (slot-value 'design-matrix))

(defmeth model-formula-proto :name-list (&optional (new nil set) )
"Sets or returns list of column names"
  (when set
        (setf (slot-value 'name-list) new))
  (slot-value 'name-list))

(defmeth model-formula-proto :block-names (&optional (new nil set))
"Sets or returns list of block names"
  (when set
        (setf (slot-value 'block-names) new))
  (slot-value 'block-names))

(defmeth model-formula-proto :block-indices (&optional (new nil set) )
"Sets or returns list of column names"
  (when set
        (setf (slot-value 'block-indices) new))
  (slot-value 'block-indices))

(defmeth model-formula-proto :compute ()
"Internal use"
(let* ( 
       (result (formula (send self :formula)))
       )
  (send self :design-matrix (first result))
  (send self :name-list (second result))
  (send self :block-names (third result))
  (send self :block-indices (fourth result))
))

(defun as-formula (x) (send model-formula-proto :new :formula x))


(provide "modelformula")
(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)))


(defun plot-steps (x y &rest graphargs &key (right t))
"Arguments: x y &rest graphargs &key (right t)
Plots a step function with the given step points. 
Default is right-continuous, optionally left continuous"
(let* (
       (plot (apply #'plot-points (append (list x y) graphargs)))
       (xstep (if right 1 0))
       (ystep (if right 0 1))
       )
  (dotimes (i (- (length x) 1))
           (send plot :add-lines (select x (list i (+ i xstep) (+ i 1)))
                 (select y (list i (+ i ystep) (+ i 1))))
           )
  plot
  )
)

(defun add-steps (x y plot &key (right t))
"Args: x y plot &key (right t)
Adds a step function to plot, which must inherit from
scatterplot-proto"
(when (null (kind-of-p plot scatterplot-proto)) (error "Not a plot"))
(let* (
       (xstep (if right 1 0))
       (ystep (if right 0 1))
      )
  (send plot :add-points x y)
  (dotimes (i (- (length x) 1))
           (send plot :add-lines (select x (list i (+ i xstep) (+ i 1)))
                 (select y (list i (+ i ystep) (+ i 1))))
           )
  plot
  )
)

(defun boxplot-factor (x y)
"Args:x y
Plot boxplots of y for each level of x"
  (boxplot-x (unique x) (groups y x))
)

(defun cumprod (s) (accumulate #'* s))

(defun unique (x) (remove-duplicates x :test #'equalp))

(defun argselect (xs x) "returns indices where xs=s" 
        (if (numberp x) 
            (which (= x xs))
          (which (mapcar #'(lambda (xi) (equalp x xi)) xs))
          ))

(defun groups (y x) 
"Args: y x
Break y into groups based on values in x"
  (let* ((xi (unique x))
         )
    (mapcar #'(lambda (i) (select y (argselect x i) )) xi)
    )
) 

(defun zero-matrix (k) (* 0 (identity-matrix k)))

(defun one-matrix (k) (outer-product (repeat 1 k) (repeat 1 k)))

(defun firsts (x) (map-elements #'(lambda (xi) (elt xi 0)) x))
(defun seconds (x) (map-elements #'(lambda (xi) (elt xi 1)) x))
(defun thirds (x) (map-elements #'(lambda (xi) (elt xi 2)) x))
(defun fourths (x) (map-elements #'(lambda (xi) (elt xi 3)) x))
(defun nths (x n) (map-elements #'(lambda (xi) (elt xi (- n 1))) x))

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

(defun col2row (x) 
"Args: x
Converts a list of columns to a list of rows"
(row-list (apply #'bind-columns x)))

(defun row2col (x) "Args: x
Converts a list of rows to a list of columns"
(column-list (apply #'bind-rows x)))

(defun append-seq (&rest seqs) 
  (apply #'append 
         (mapcar #'(lambda (x) (coerce x 'list)) seqs)))

(defun flatten (listlist) (apply #'append listlist))

(defun list-of-seqs-p (x) 
     (when (listp x) (when (sequencep (car x)) (numberp (car (car x))))))

(defun gtrp (x y) (if (numberp x) (> x y) (STRING> x y)))

;; used in the examples.
(defun wait () 
  (format t "--Press enter to continue--") 
  (read-line))

;  the first four elements of a and b are time status entry stratum
; stratum could be a string

(defun after (a b) 
  (if (equalp (elt a 3) (elt b 3)) 
        (if (= (elt a 0) (elt b 0))
            (if (= (elt a 2) (elt b 2)) 
                (< (elt a 1) (elt b 1))
              (> (elt a 2) (elt b 2))
              )
          (> (elt a 0) (elt b 0))
          )
    (gtrp (elt a 3) (elt b 3))
    )
)

(defun tied (a b) (if (and (= (elt a 0) (elt b 0))
                           (= 1 (elt a 1) (elt b 1))
                           (= (elt a 2) (elt b 2)) 
                           (equalp (elt a 3) (elt b 3)))
                      1 0))


(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 correlation (x y) 
  (let* ((m (covariance-matrix x y)))
    (/ (aref m 0 1) (sqrt (prod (diagonal m)))) 
    )
)

(defun n+ (i j n) (if (null i) n (+ i j)))

(defun interpolate (xnew x y &key (right t))
"Interpolates the step function with corners at (x y) at the values 
 in xnew. Default is a right-continuous step function like the
 Kaplan-Meier and Nelson-Aalen estimators"
(let* ((fudge (if right -1 0))
       (n (- (length y) 1)))
  (map-elements 
   #'(lambda (xnewi) (elt y (n+ (position xnewi x :test #'<)  fudge n)))
   xnew)
  )
)
(defproto cox-tie-proto) 
(defproto breslow-ties () () cox-tie-proto)
(defproto counting-process-ties () () cox-tie-proto)
(defproto peto-ties () () cox-tie-proto)
(defproto efron-ties () () cox-tie-proto)
(defmeth cox-tie-proto  :increment (zlist betahat ns0 ns1 ns2)
(let* (      
       (expbtz (map-elements #'(lambda (z) (exp (inner-product z betahat))) zlist))
       (zbtz (* expbtz zlist))
       )
   (list (- (apply #'+ zlist) (* (length zlist) (/ ns1 ns0)))
         (* (length zlist) (- (/ ns2 ns0) (outer-product (/ ns1 ns0) (/ ns1 ns0)))))
  )
)
(defmeth efron-ties :increment (zlist betahat ns0 ns1 ns2)
(let* (      
       (expbtz (map-elements #'(lambda (z) (exp (inner-product z betahat))) zlist))
       (zbtz (* expbtz zlist))
       (z2btz (* expbtz (map-elements #'outer-product zlist zlist)))
       (n (length zlist))
       (w (/ (iseq  n) n))
       (zsum (apply #'+ zbtz))
       (z2sum (apply #'+ z2btz))
       (ebtzsum (apply #'+ expbtz))
       (zbars (/ (mapcar #'(lambda (wi) (- ns1 (* wi zsum))) w) 
                 (- ns0 (* w ebtzsum))))
       (z2bars (/ (mapcar #'(lambda (wi) (- ns2 (* wi  z2sum))) w) 
                  (- ns0 (* w ebtzsum))))
       )
   (list (- (apply #'+ zlist) (apply #'+ zbars))  
         (- (apply #'+ z2bars) (apply #'+ (mapcar #'outer-product zbars zbars))))
   )
)

(defmeth cox-tie-proto :dLambda (zlist betahat ns0)
"Args: zlist betahat ns0 
Computes the increment of the baseline cumulative hazard for the tied failures.
equiv-at-risk is the sum of exp(b'z) just before this time, so
(/ nS0) is the increment if there is only one failure"
  (/ (length zlist) ns0)
)
(defmeth efron-ties :dLambda (zlist betahat ns0)
"Args: zlist betahat ns0 
Computes the increment of the baseline cumulative hazard for the tied failures.
equiv-at-risk is the sum of exp(b'z) just before this time, so
(/ nS0) is the increment if there is only one failure"
(let*  (      
        (expbtz (map-elements #'(lambda (z) (exp (inner-product z betahat))) 
                              zlist))
        (n (length zlist))
        (w (/ (iseq  n) n))
        (ebtzsum (apply #'+ expbtz))
        )
  (sum (/ (- ns0 (* w ebtzsum))))
  )
)
(defmeth cox-tie-proto :dloglike (zlist betahat ns0)
"Args: zlist betahat ns0
Computes the increment in the loglikelihood"
(let* ((btz (map-elements #'(lambda (z) (inner-product z betahat)) zlist))
       )
  (- (apply #'+ btz) (* (length zlist) (log ns0)))
  )
)
(defmeth efron-ties :dloglike (zlist betahat ns0)
"Args: zlist betahat ns0
Computes the increment in the loglikelihood"
(let*  (      
        (btz (map-elements #'(lambda (z) (inner-product z betahat)) zlist))
        (expbtz (exp btz))
        (n (length zlist))
        (w (/ (iseq  n) n))
        (ebtzsum (apply #'+ expbtz))
        (btzsum (apply #'+ btz))
        )
  (- btzsum (sum (log  (- ns0 (* w ebtzsum)))))
  )
)
(defmeth cox-tie-proto :schoenfeld (zlist betahat ns0 ns1)
"Args: zlist betahat
Computes the appropriate Schoenfeld residual(s) for the tied failures"
(let* ((expbtz (map-elements #'(lambda (z) (exp (inner-product z betahat))) zlist))
       (n (length zlist))
       )
   (- zlist (repeat (list (/ ns1 ns0)) n))
  )
)
(defmeth efron-ties :schoenfeld (zlist betahat ns0 ns1)
"Args: zlist betahat
Computes the appropriate Schoenfeld residual(s) for the tied failures"
(let*  (      
        (expbtz (map-elements #'(lambda (z) (exp (inner-product z betahat))) 
                              zlist))
        (zbtz (* expbtz zlist))
        (n (length zlist))
        (w (/ (iseq  n) n))
        (zsum (apply #'+ zbtz))
        (ebtzsum (apply #'+ expbtz))
        (zbars (/ (mapcar #'(lambda (wi) (- ns1 (* wi zsum))) w) 
                  (- ns0 (* w ebtzsum))))
        )
     (- zlist zbars)
  )
)
(defmeth cox-tie-proto :prettyprint (&optional (stream t))
  (format stream "Tie-handling prototype~%")
)
(defmeth peto-ties :prettyprint (&optional (stream t))
  (format stream "Peto's (counting process) approximation~%")
)
(defmeth breslow-ties :prettyprint (&optional (stream t))
  (format stream "Breslow's (counting process) approximation~%")
)
(defmeth counting-process-ties :prettyprint (&optional (stream t))
  (format stream "Counting process (Breslow's approximation)~%")
)
(defmeth efron-ties :prettyprint (&optional (stream t))
  (format stream "Efron's approximation~%")
)
(defproto cox-regression-proto 
  '(times entry status x formula beta id strata weights
         epsilon epsilon-dev count-limit verbose recycle
         eta loglike agnostic-covariance-matrix  
         model-covariance-matrix tie-method recycle agnostic
         needs-computing cumulative-hazard
         martingales
         schoenfeld
         scores suppress-warnings
         predictor-names case-labels response-name
         statusint keyint entryint tieint comments time-dependent zbars
         paranoia
)
  '() 
  *object*
  "Cox regression model")
(defmeth cox-regression-proto :times (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns failure/censoring times"
  (when set
        (setf (slot-value 'times) new)
        (send self :needs-computing t))
  (slot-value 'times))


(defmeth cox-regression-proto :count-limit (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns iteration limit"
  (when set
        (setf (slot-value 'count-limit) new)
        (send self :needs-computing t))
  (slot-value 'count-limit))

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

(defmeth cox-regression-proto :intercept () nil)

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

(defmeth cox-regression-proto :id (&optional (new nil set))
"Message args: &optional (new nil set)
Sets or returns case id"
(when set
      (setf (slot-value 'id) id))
(slot-value 'id))

(defmeth cox-regression-proto :beta () 
  "Returns coefficient estimates"
  (slot-value 'beta))

(defmeth cox-regression-proto :keyint (&optional (new nil set))
"INTERNAL [ Message args: (&optional new)
  Sets or returns index from munged times into covariates]"
  (when set
        (setf (slot-value 'keyint) new)
        (send self :needs-computing t))
  (slot-value 'keyint))


(defmeth cox-regression-proto :statusint (&optional (new nil set))
  "INTERNAL [ Message args: (&optional new)
 Sets or returns failure status at munged times]"
  (when set
        (setf (slot-value 'statusint) new)
        (send self :needs-computing t))
  (slot-value 'statusint))

(defmeth cox-regression-proto :entryint (&optional (new nil set))
  "INTERNAL [Message args: (&optional new)
 Sets or returns entry status at munged times]"
  (when set
        (setf (slot-value 'entryint) new)
        (send self :needs-computing t))
  (slot-value 'entryint))

(defmeth cox-regression-proto :tieint (&optional (new nil set))
  "INTERNAL [Message args: (&optional new)
 Sets or returns tie status at munged times]"
  (when set
        (setf (slot-value 'tieint) new)
        (send self :needs-computing t))
  (slot-value 'tieint))

(defmeth cox-regression-proto :weights (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns case weights. These have no effect whatsoever."
  (when set
        (setf (slot-value 'weights) new)
        (send self :needs-computing t))
  (slot-value 'weights))

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

(defmeth cox-regression-proto :agnostic (&optional (new nil set))
  "Message args: (&optional new)
 Use agnostic covariance estimator"
  (when set
        (setf (slot-value 'agnostic) new)
        (send self :needs-computing t))
  (slot-value 'agnostic))


(defmeth cox-regression-proto :tie-method (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns tie-method"
  (when set
        (setf (slot-value 'tie-method) new)
        (send self :needs-computing t))
  (slot-value 'tie-method))

(defmeth cox-regression-proto :paranoia (&optional (new nil set))
   "Message args: (&optional new)"
   (when set
         (setf (slot-value 'paranoia) new)
         )
   (slot-value 'paranoia))
 

(defmeth cox-regression-proto :num-cases ()
  (length (unique (send self :id))))

(defmeth cox-regression-proto :x (&optional x)
  (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 cox-regression-proto :formula (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns model formula"
  (when set
        (setf (slot-value 'formula) new)
        (send self :needs-computing t))
  (slot-value 'formula))

(defmeth cox-regression-proto :comments (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns model description text"
  (when set
        (setf (slot-value 'comments) new)
        )
  (slot-value 'comments))

(defmeth cox-regression-proto :epsilon (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns convergence criterion for beta"
  (when set
        (setf (slot-value 'epsilon) new)
        (send self :needs-computing t))
  (slot-value 'epsilon))

(defmeth cox-regression-proto :epsilon-dev (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns convergence criterion for loglikelihood"
  (when set
        (setf (slot-value 'epsilon-dev) new)
        (send self :needs-computing t))
  (slot-value 'epsilon-dev))

(defmeth cox-regression-proto :recycle (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns indicator for reusing old values when refitting"
  (when set
        (setf (slot-value 'recycle) new)
        (send self :needs-computing t))
  (slot-value 'recycle))

(defmeth cox-regression-proto :id (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns case identifier (important for agnostic
covariance estimates or to get single residuals per case
with time-dependent covariates)"
  (when set
        (setf (slot-value 'id) new)
        (send self :needs-computing t))
  (slot-value 'id))

(defmeth cox-regression-proto :covariance-matrix (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns covariance matrix"
  (if (send self :agnostic) 
      (send self :agnostic-covariance-matrix)
    (send self :model-covariance-matrix))
  )

(defmeth cox-regression-proto :model-covariance-matrix (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns model-based covariance matrix"
  (when set
        (setf (slot-value 'model-covariance-matrix) new)
        )
  (slot-value 'model-covariance-matrix))

(defmeth cox-regression-proto :agnostic-covariance-matrix (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns agnostic (model-robust) covariance matrix"
  (when set
        (setf (slot-value 'agnostic-covariance-matrix) new)
        )
  (slot-value 'agnostic-covariance-matrix))

(defmeth cox-regression-proto :verbose (&optional (new nil set))
  "Message args: (&optional new)
 Sets or returns verboseness control (default is t to 
report iteration information)"
  (when set
        (setf (slot-value 'verbose) new)
        (send self :needs-computing t))
  (slot-value 'verbose))

(defmeth cox-regression-proto :suppress-warnings (&optional (new nil set))
  (when set
        (setf (slot-value 'suppress-warnings) new)
        (send self :needs-computing t))
  (slot-value 'suppress-warnings))

(defmeth cox-regression-proto :loglike (&optional (new nil set))
  "Message args: (&optional new)
Sets or returns loglikelihood"
  (when set
        (setf (slot-value 'loglike) new)
        (send self :needs-computing t))
  (slot-value 'loglike))

(defmeth cox-regression-proto :needs-computing (&optional (new nil set))
  "Message args: (&optional new)
Flag: t if model needs updating"
  (when set
        (setf (slot-value 'needs-computing) new)
        )
  (slot-value 'needs-computing))

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

(defmeth cox-regression-proto :model-coef-standard-errors ()
(sqrt (diagonal (send self :model-covariance-matrix))))

(defmeth cox-regression-proto :agnostic-coef-standard-errors ()
 (sqrt (diagonal (send self :agnostic-covariance-matrix))))

(defmeth cox-regression-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* ((pfixed (array-dimension (send self :x) 1))
        (ptime (length (first (send self :time-dependent))))
        (p (+ pfixed ptime))
        (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 cox-regression-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 cox-regression-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))

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

(defmeth cox-regression-proto :fit-values ()
"Message args: ()
Returns the fitted linear predictor values for the model."
  (matmult (send self :x) (send self :beta)))

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

(defmeth cox-regression-proto :set-eta (&optional val)
  (if val
      (setf (slot-value 'eta) val)
      (setf (slot-value 'eta)
            (matmult (send self :x) (send self :beta))))
  (slot-value 'eta))

(defmeth cox-regression-proto :time-dependent (&optional (new nil set))
"Message args: &optional (new nil set)
Sets or returns the information on continuously time-dependent
covariates, consisting of a list containing a list of functions and a covariate matrix"
     (when set (setf (slot-value 'time-dependent) new))
     (slot-value 'time-dependent)
)

(defmeth cox-regression-proto :fit-hazard-ratio (&optional (eta (send self :eta)))
"Message args: (&optional (eta (send self :eta)))
Returns mean values for current or supplied ETA."
(if (null eta) (exp (send self :eta))  (exp eta)))

(defmeth cox-regression-proto :cumulative-hazard (&optional (new nil set))
 (when set
        (setf (slot-value 'cumulative-hazard) new)
        )
  (slot-value 'cumulative-hazard))

(defmeth cox-regression-proto :schoenfeld-residuals (&optional (new nil set))
 (when set
        (setf (slot-value 'schoenfeld) new)
        )
  (slot-value 'schoenfeld))

(defmeth cox-regression-proto :zbars (&optional (new nil set))
"Args: (new nil set)
Sets or returns a list whose elements are the time, stratum number, number at risk 
in the stratum and the weighted mean covariate for each failure"
 (when set
        (setf (slot-value 'zbars) new)
        )
  (slot-value 'zbars))

(defmeth cox-regression-proto :scores (&optional (new nil set))
"Message ars: (&optional new)
Sets or returns martingale score residuals"
  (when set
        (setf (slot-value 'scores) new)
        )
  (slot-value 'scores))
(defmeth cox-regression-proto :initialize-search () 
  (let*
      ((np (+ (length (column-list (send self :x))) 
              (length (first (send self :time-dependent)))))) 
    (send self :set-beta (repeat 0 np))) )

(defmeth cox-regression-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 :beta) (send self :recycle))
            (send self :initialize-search))
    (do ((count 1 (+ count 1))
         (beta (send self :beta)  (send self :beta))
         (last-beta -1 beta)
         (dev  (- (/ machine-epsilon)) (send self :loglike))
         (last-dev  0 dev))
        ((or (> count maxcount) 
             (< (max (abs (/ (- beta last-beta)
                             (pmax (abs last-beta) low-lim))))
                epsilon)
             (< (abs (- dev last-dev)) epsilon-dev)) 
         (when (and  (> count maxcount) (null (send self :suppress-warnings))) 
               (format t 
                       "~%~%Failed to converge in ~d iterations~%Treat results with caution ~%" 
                       maxcount)))
        (send self :paranoia nil)
        (if (car (send self :time-dependent))
            (send self :compute-tstep) 
          (send self :compute-step))
         (if (and (> count 1) (or (send self :paranoia) (< (- dev last-dev) (- epsilon-dev))))
          (progn 
           (send self :set-beta (/ (+ beta last-beta) 2))
           (format t "Possible numerical instability: step-halving~%")))

        (if verbose 
            (format t "Iteration ~d: log-likelihood = ~,6g~%" 
                    count (send self :loglike)))
        )
    )
  (if (car (send self :time-dependent)) 
      (send self :compute-tmartingale)
    (send self :compute-martingale))
  (when  (send self :agnostic) 
         (send self :compute-agnostic-covariance)
         )
  (send self :needs-computing nil)
)
(defmeth cox-regression-proto :compute-tstep 
  ()
"Does one step of Newton-Raphson algorithm for models with 
continuously time-dependent covariates"
(let* (
       (keyint (send self :keyint))
       (statusint (send self :statusint))
       (entryint (send self :entryint))
       (tieint (send self :tieint))
       (times (send self :times))
       (strata (send self :strata))
       (stratlist (unique strata))
       (x-list (row-list (send self :x)))
       (riskmin 1) 
       (riskmax 1) 
       (oldbeta (send self :beta))
       (timecovs (send self :time-dependent))
       (timefuns (first timecovs))
       (timexs (second timecovs))
       (pfixed (length (first x-list)))
       (ptime (if timecovs (length (first timecovs)) 0))
       (p (+ pfixed ptime))
       (nz (length x-list))
       (nt (length keyint))
       (dcumhaz nil)
       (schoenfelds nil)
       (loglike 0)
       (score (repeat 0 p))
       (varscore (matrix (list p p) (repeat 0 (* p p))))
       (current-ties nil)
       (tie-method (send self :tie-method))
       (last-stratum (elt strata (elt keyint 0)))
       (this-stratum (elt strata (elt keyint 0)))
       (zbars nil)
       (riskset nil)
       (junk (dotimes (i nt)
                      (setf this-stratum (elt strata (elt keyint i)))
                      (cond ((equalp this-stratum last-stratum) 
                             (setf last-stratum this-stratum 
                                   ))
                            (t 
                             (setf riskset nil 
                                   last-stratum this-stratum 
                                   ))
                            )
                      (cond ((= (elt tieint i) 1) ;;one of a set of ties
                             (setf current-ties 
                                   (cons (elt keyint i) current-ties))
                             )
                            (current-ties ;; last of a tied set of failures
                             (unless (= (elt statusint i) 1) 
                                     (error "Can't happen"))
                             (setf current-ties 
                                   (cons (elt keyint i) current-ties))
                               ;;; recalculate nS0, nS1, nS2 first      
                             (setf nS0 0)
                             (setf nS1 (repeat 0 p))
                             (setf nS2 (zero-matrix p))
                             (setf riskset (append current-ties riskset))
                             (dolist (this-obs riskset)
                                     (let* ((timez (mapcar 
                                                    #'(lambda (f) 
                                                        (funcall f  
                                                                 (elt timexs this-obs) 
                                                                 (elt times (elt keyint i)))) 
                                                    timefuns))
                                            (z (append-seq 
                                                (select x-list this-obs) 
                                                timez))
                                            (ebtz (exp (inner-product z oldbeta)))
                                            (setf riskmin (min riskmin ebtz))
                                            (setf riskmax (max riskmax ebtz))
                                            )
                                       (setf nS0 (+ nS0 ebtz))
                                       (setf nS1 (+ nS1 (* z ebtz)))
                                       (setf nS2 (+ nS2 (* ebtz (outer-product z z))))
                                       )                                       
                                     )
                             (let* ((timexlist (select timexs current-ties))
                                    (timezlist (mapcar 
                                                #'(lambda (keyi) 
                                                    (mapcar 
                                                     #'(lambda (f) (funcall f  
                                                                            (elt timexs keyi) 
                                                                            (elt times keyi))) 
                                                     timefuns)) 
                                                current-ties))
                                    (zlist (row-list 
                                            (bind-columns 
                                             (apply #'bind-rows (select x-list current-ties)) 
                                             (apply #'bind-rows timezlist))))
                                    (temp (send tie-method 
                                                :increment zlist  oldbeta nS0 nS1 nS2))
                                    (thistime (select times (select keyint i)))
                                    )
                               (setf score (+ score (first temp)))
                               (setf varscore (+ varscore (second temp)))
                               (setf loglike (+ loglike 
                                                (send tie-method 
                                                      :dloglike zlist oldbeta ns0)))
                               (setf dcumhaz
                                     (cons 
                                      (list thistime  
                                            (send tie-method 
                                                  :dLambda  zlist oldbeta nS0)
                                            this-stratum) 
                                      dcumhaz))
                               (setf schoenfelds 
                                     (append 
                                      (mapcar #'(lambda (si)  
                                                  (list thistime 
                                                        si 
                                                        this-stratum)) 
                                              (send tie-method 
                                                    :schoenfeld zlist oldbeta ns0 ns1)) 
                                      schoenfelds))
                               (setf zbars (cons 
                                            (list thistime 
                                                  this-stratum 
                                                  (length riskset) 
                                                  (repeat (/ nS1 nS0) 
                                                          (length current-ties))
                                                  ) 
                                            zbars))
                               (setf current-ties nil)
                               )
                             )
                            ((= (elt statusint i) 1)  ;; untied failure
                             (setf riskset (cons (elt keyint i) riskset))
                               ;;; recalculate nS0, nS1, nS2 first      
                             (setf nS0 0)
                             (setf nS1 (repeat 0 p))
                             (setf nS2 (zero-matrix p))
                             (dolist (this-obs  riskset)
                                     (let* ((timez (mapcar 
                                                    #'(lambda (f) 
                                                        (funcall f  
                                                                 (elt timexs this-obs) 
                                                                 (elt times (elt keyint i)))) 
                                                    timefuns))
                                            (z (append-seq 
                                                (select x-list this-obs) 
                                                timez))
                                            (ebtz (exp (inner-product z oldbeta)))
                                            (setf riskmin (min riskmin ebtz))
                                            (setf riskmax (max riskmax ebtz))
                                            )
                                       (setf nS0 (+ nS0 ebtz))
                                       (setf nS1 (+ nS1 (* z ebtz)))
                                       (setf nS2 (+ nS2 (* ebtz (outer-product z z))))
                                       )                                       
                                     )
                             
                             (let* ((timez (mapcar 
                                            #'(lambda (f) 
                                                (funcall f  
                                                         (elt timexs (elt keyint i)) 
                                                         (elt times (elt keyint i)))) 
                                            timefuns))
                                    (z (append-seq 
                                        (select x-list (select keyint i)) 
                                        timez))
                                    (btz (inner-product z oldbeta))
                                    (ebtz (exp btz))
                                    (thistime (select times (select keyint i)))
                                    )
                               (setf current-ties nil)
                               (setf score (+ score (- z (/ nS1 nS0))))
                               (setf varscore (+ varscore 
                                                 (/ nS2 nS0) 
                                                 (/ (outer-product (- nS1) nS1) 
                                                    (* nS0 nS0))))
                               (setf dcumhaz (cons 
                                              (list thistime (/ ns0) this-stratum) 
                                              dcumhaz))
                               (setf schoenfelds (cons 
                                                  (list thistime 
                                                        (- z (/ nS1 nS0)) 
                                                        this-stratum) 
                                                  schoenfelds))
                               (setf zbars (cons 
                                            (list thistime 
                                                  this-stratum 
                                                  (length riskset) 
                                                  (/ nS1 nS0)) 
                                            zbars))
                               (setf loglike (+ loglike btz (- (log nS0))))
                               )
                             )
                            ((= (elt entryint i) 1) ;; censoring
                             (setf riskset 
                                   (cons (elt keyint i) riskset))
                             )
                            ((= (elt entryint i) -1) ;; entry
                             (if (position (elt keyint i) riskset)
                                 (setf riskset (remove (elt keyint i) riskset))
                               (error "Entry time with no corresponding exit time"))
                             )
                            (t (error "Can't happen in :compute-step"))
                            )
                      ))
       (invL (inverse (first (chol-decomp varscore))))
       (inv (matmult (transpose invL) invL))
       (cumhaz (map-elements 
                #'(lambda (stratum) 
                    (let* ((these (select zbars 
                                          (argselect (seconds zbars) stratum) )))
                      (if these
                          (list stratum 
                                (append-seq #(0) (firsts these))  
                                (append-seq #(0) (cumsum (seconds these))))
                        (list stratum '(0 0) '(0 0)))))
                    stratlist))
       (rotzbars (map-elements 
                  #'(lambda (stratum) 
                      (let* ((these (select zbars 
                                            (argselect (seconds zbars) stratum) ))) 
                        (if these
                            (list stratum 
                                  (firsts these) 
                                  (fourths these) 
                                  (thirds these))
                          (list stratum '(0) (list (repeat 0 p)) '(0)))))
                      stratlist))
       )
  (send self :set-beta (+ oldbeta (matmult score inv)))
  (send self :model-covariance-matrix inv )
  (send self :loglike loglike)
  (send self :cumulative-hazard cumhaz)
  (send self :schoenfeld-residuals schoenfelds)
  (send self :zbars rotzbars) 
  (when (and *cox-paranoia* ( < (/ riskmin riskmax) (sqrt machine-epsilon)))
  (send self :paranoia t))
  )
)
(defmeth cox-regression-proto :compute-step ()
"Does one step of Newton-Raphson algorithm"
(let* (
       (keyint (send self :keyint))
       (statusint (send self :statusint))
       (entryint (send self :entryint))
       (tieint (send self :tieint))
       (times (send self :times))
       (strata (send self :strata))
       (stratlist (unique strata))
       (x-list (row-list (send self :x)))
       (riskmin 1)
       (riskmax 1)
       (oldbeta (send self :beta))
       (p (length (first x-list)))
       (nz (length x-list))
       (nt (length keyint))
       (dcumhaz nil)
       (schoenfelds nil)
       (nS0 0)
       (nrisk 0)
       (loglike 0)
       (nS1 (repeat 0 p))
       (nS2 (matrix (list p p) (repeat 0 (* p p))))
       (score (copy-seq nS1))
       (varscore (copy-array nS2))
       (current-ties 0)
       (zbars nil)
       (tie-method (send self :tie-method))
       (last-stratum (elt strata (elt keyint 0)))
       (this-stratum (elt strata (elt keyint 0)))
       (junk (dotimes (i nt)
                      ;; work backwards in time, accumulating information as we go
                      (setf this-stratum (elt strata (elt keyint i)))
                      (cond ((equalp this-stratum last-stratum) 
                             (setf last-stratum this-stratum 
                                   ))
                            (t 
                             (setf nrisk 0 
                                   nS0 0 
                                   nS1 (repeat 0 p) 
                                   nS2 (matrix (list p p) (repeat 0 (* p p))) 
                                   last-stratum this-stratum 
                                   ))
                            )
                      (cond ((= (elt tieint i) 1)  ; tied failure
                             (let* (
                                    (z (select x-list (select keyint i)))
                                    (ebtz (exp (inner-product z oldbeta)))
                                    )
                               (setf riskmin (min riskmin ebtz))
                               (setf riskmax (max riskmax ebtz))
                               (setf nS0 (+ nS0 ebtz))
                               (setf nS1 (+ nS1 (* z ebtz)))
                               (setf nS2 (+ nS2 (* ebtz (outer-product z z))))
                               (setf current-ties (+ current-ties 1))
                               (setf nrisk (+ nrisk 1))
                               )
                             )
                            ((> current-ties 0)  ; last of tied set of failures
                             (let* (
                                    (z (select x-list (select keyint i)))
                                    (ebtz (exp (inner-product z oldbeta)))
                                    )
                               (setf riskmin (min riskmin ebtz))
                               (setf riskmax (max riskmax ebtz))
                               (setf nS0 (+ nS0 ebtz))
                               (setf nS1 (+ nS1 (* z ebtz)))
                               (setf nS2 (+ nS2 (* ebtz (outer-product z z))))
                               (setf nrisk (+ 1 nrisk))
                               (setf current-ties (+ current-ties 1)))
                             (let* ((keylist (select keyint (- i (iseq  current-ties))))
                                    (zlist (select x-list keylist))
                                    (temp (send tie-method 
                                                :increment zlist  oldbeta nS0 nS1 nS2))
                                    (thistime (select times (select keyint i)))
                                    )
                               (setf score (+ score (first temp)))
                               (setf varscore (+ varscore (second temp)))
                               (setf loglike (+ loglike 
                                                (send tie-method :dloglike 
                                                      (select x-list keylist) 
                                                      oldbeta ns0)))
                               (setf dcumhaz (cons 
                                              (list thistime  
                                                    (send tie-method :dLambda  
                                                          zlist 
                                                          oldbeta 
                                                          nS0) 
                                                    this-stratum) 
                                              dcumhaz))
                               (setf schoenfelds (append 
                                                  (mapcar #'(lambda (si)  
                                                              (list thistime 
                                                                    si 
                                                                    this-stratum)) 
                                                          (send tie-method :schoenfeld 
                                                                zlist 
                                                                oldbeta 
                                                                ns0 
                                                                ns1)) 
                                                  schoenfelds))
                               (setf zbars (cons 
                                            (repeat (list thistime 
                                                          this-stratum 
                                                          nrisk 
                                                          (/ nS1 nS0)) 
                                                    current-ties) 
                                            zbars))
                               (setf current-ties 0)
                               )
                             )
                            ((= (elt statusint i) 1)  ; untied failure
                             (let* ((z (select x-list (select keyint i)))
                                    (btz (inner-product z oldbeta))
                                    (ebtz (exp btz))
                                    (thistime (select times (select keyint i)))
                                    )
                               (setf riskmin (min riskmin ebtz))
                               (setf riskmax (max riskmax ebtz))
                               (setf nS0 (+ nS0 ebtz))
                               (setf nS1 (+ nS1 (* z ebtz)))
                               (setf nS2 (+ nS2 (* ebtz (outer-product z z))))
                               (setf current-ties 0)
                               (setf nrisk (+ nrisk 1))
                               (setf score (+ score (- z (/ nS1 nS0))))
                               (setf varscore (+ varscore 
                                                 (/ nS2 nS0) 
                                                 (/ (outer-product (- nS1) nS1) 
                                                    (* nS0 nS0))))
                               (setf dcumhaz (cons (list thistime 
                                                         (/ ns0) 
                                                         this-stratum) 
                                                   dcumhaz))
                               (setf schoenfelds (cons (list thistime 
                                                             (- z (/ nS1 nS0)) 
                                                             this-stratum) 
                                                       schoenfelds))
                               (setf zbars (cons (list thistime 
                                                       this-stratum 
                                                       nrisk 
                                                       (/ nS1 nS0)) 
                                                 zbars))
                               (setf loglike (+ loglike btz (- (log nS0))))
                               )
                            )
                            ((= (elt entryint i) 1) ;censoring
                             (let* ((z (select x-list (select keyint i)))
                                    (ebtz (exp (inner-product z oldbeta)))
                                    )
                               (setf riskmin (min riskmin ebtz))
                               (setf riskmax (max riskmax ebtz))
                               (setf nS0 (+ nS0 ebtz))
                               (setf nS1 (+ nS1 (* z ebtz)))
                               (setf nS2 (+ nS2 (* ebtz (outer-product z z))))
                               (setf nrisk (+ nrisk 1))
                               )
                             )
                            ((= (elt entryint i) -1) ;left truncation
                             (let* ((z (select x-list (select keyint i)))
                                    (ebtz (exp (inner-product z oldbeta)))
                                   )
                               (setf riskmin (min riskmin ebtz))
                               (setf riskmax (max riskmax ebtz))
                               (setf nrisk (- nrisk 1))
                               (setf nS0 (- nS0 ebtz))
                               (setf nS1 (- nS1 (* z ebtz)))
                               (setf nS2 (- nS2 (* ebtz (outer-product z z))))
                               )
                             )
                            (t (error "Can't happen in :compute-step")) 
                            )
                      ))
       (invL (inverse (first (chol-decomp varscore))))
       (inv (matmult (transpose invL) invL))
       (cumhaz (map-elements 
                #'(lambda (stratum) 
                    (let* ((these (select dcumhaz 
                                          (argselect (thirds dcumhaz) stratum)))
                           )
                      (if these
                          (list stratum 
                                (append-seq  #(0) (firsts these)) 
                                (append-seq  #(0) (cumsum (seconds these))))
                        (list stratum '(0 0) '(0 0))
                        )
                      ))
                stratlist))
       (rotzbars (map-elements 
                  #'(lambda (stratum) 
                      (let* ((these (select zbars 
                                            (argselect (seconds zbars) stratum) )))
                        (if these
                            (list stratum 
                                  (firsts these) 
                                  (fourths these) 
                                  (thirds these))
                          (list stratum '(0) (list (repeat 0 p))  '(0))
                          )
                        ))
                  stratlist))
       )
  (send self :set-beta (+ oldbeta (matmult score inv)))
  (send self :model-covariance-matrix inv )
  (send self :loglike loglike)
  (send self :cumulative-hazard cumhaz)
  (send self :schoenfeld-residuals schoenfelds)
  (send self :zbars rotzbars)
  (when (and *cox-paranoia* ( < (/ riskmin riskmax) (sqrt machine-epsilon)))
    (send self :paranoia t))
  )
)
(provide "martingale")

(defun outer-square (seq) (outer-product seq seq))

(defmeth cox-regression-proto :compute-martingale ()
"Martingale and score residuals collapsed over entry/exit but not id"
(let* (
       (keyint (send self :keyint))
       (statusint (send self :statusint))
       (entryint (send self :entryint))
       (tieint (send self :tieint))
       (times (send self :times))
       (entrytimes (send self :entry))
       (strata (send self :strata))
       (stratlist (unique strata))
       (x-list (row-list (send self :x)))
       (p (length (first x-list)))
       (nz (length x-list))
       (nt (length keyint))
       (cumhaz (send self :cumulative-hazard))
       (zbars (send self :zbars))
       (id (send self :id))
       (beta (send self :beta))
       (dcumhaz  (map-elements 
                  #'(lambda (stratum)  (list (first stratum) 
                                             (second stratum) 
                                             (difference (third stratum)))) 
                  cumhaz))
       (zbardL (map-elements  #'(lambda (hazs zbars)  
                                  (list (first hazs) 
                                        (second hazs) 
                                        (* (third hazs) (third zbars)))) 
                              dcumhaz 
                              zbars))
       (intzbardL (map-elements  #'(lambda (stratum) 
                                     (list (first stratum) 
                                           (second stratum) 
                                           (cumsum (third stratum)))) 
                                 zbardL))
       (rval (row-list (bind-columns (repeat 0 nz) (* 0 (send self :x)))))
       (fn #'(lambda (caseno) 
               (let* ((key (elt keyint caseno))
                      (idno (elt id key))
                      (stratno (elt strata key))
                      (status (elt statusint caseno))
                      (entry (elt entryint caseno))
                      (thistime (if (= entry 1) 
                                    (elt times key) 
                                  (elt entrytimes key)))
                      (zi (elt x-list key))
                      (risk (exp (inner-product zi beta)))
                      (whichtminus (position 
                                    thistime 
                                    (second (find stratno 
                                                  zbars 
                                                  :test #'(lambda (x y) 
                                                            (equalp  x (elt y 0))))) 
                                    :test #'>= 
                                    :from-end t))
                      (zbarti (if whichtminus 
                                  (elt (third (find 
                                               stratno 
                                               zbars 
                                               :test #'(lambda (x y) 
                                                         (equalp x (elt y 0))))) 
                                       whichtminus)
                                (* 0 zi)))
                      (Lti (* risk (if whichtminus 
                                       (elt (third (find 
                                                    stratno 
                                                    cumhaz 
                                                    :test #'(lambda (x y) 
                                                              (equalp  x (elt y 0))))) 
                                            (+ 1 whichtminus)) 0)))
                      (intzbardLti (* risk 
                                      (if whichtminus 
                                          (elt (third (find 
                                                       stratno 
                                                       intzbardL 
                                                       :test #'(lambda (x y) 
                                                                 (equalp x (elt y 0))))) 
                                               whichtminus) 
                                        0)))
                      (Mti (cons 
                            (- status Lti) 
                            (coerce (+ 
                                     (* status (- zi zbarti)) 
                                     intzbardLti 
                                     (- (* zi Lti))) 
                                    'list)))
                      )
                 (setf (select rval key) (+ (select rval key) (* entry Mti)))
                 )
               t))
       (junk (map-elements fn (iseq nt)))
       (mgales (column-list (apply #'bind-rows rval)))
       )
  (send self :martingales (first mgales))
  (send self :scores (apply #'bind-columns (cdr mgales)))
  )          
t
)


(defmeth cox-regression-proto :compute-tmartingale ()
"Martingale and score residuals collapsed over id for
 general time-dependence"
(format t 
        "Warning: martingale residuals not yet available for general time-dependence\n")
)


(defmeth cox-regression-proto :compute-agnostic-covariance ()
"INTERNAL: compute and set appropriate agnostic covariance matrix"
(let* ((id (send self :id))
       (unique-id (coerce (unique id) 'list))
       (scores (or 
                (row-list (send self :scores)) 
                (error "Can't do agnostic covariance - no martingale residuals")))
       (sumsq (cond ( (= (length id) (length unique-id))
                      (apply #'+ 
                             (mapcar #'outer-square scores)))
                    (t 
                     (apply #'+ 
                            (mapcar #'outer-square 
                                    (mapcar #'(lambda (i) 
                                                (apply #'+ 
                                                       (select scores 
                                                               (argselect id i)))) 
                                            unique-id))))
                    ))
       (modelv (send self :model-covariance-matrix))
       )
  (send self :agnostic-covariance-matrix (matmult modelv sumsq modelv))
))

(defmeth cox-regression-proto :display ()
"Prints the model summary"
   (if (send self :needs-computing) (send self :compute))
  (let ((coefs (coerce (send self :beta) 'list))
        (se-s (if (send self :agnostic) 
                  (send self :agnostic-coef-standard-errors) 
                (send self :coef-standard-errors)))
        (x (send self :x))
        (p-names (send self :predictor-names)))
    (format t (send self :comments))
    (if (> (length (unique (send self :strata))) 1) 
        (format t "~% Stratified ") 
      (format t "~%"))
    (format t "Maximum Partial Likelihood Estimates:~2%")
    (format t "~25tCoefficient~40t Std Error~%")
    (when (send self :intercept)
          (format t "Constant~25t~13,6g~40t(~,6g)~%" (car coefs) (car se-s))
          (setf coefs (cdr coefs))
          (setf se-s (cdr se-s)))
    (dotimes (i (length coefs)) 
             (format t "~a~25t~13,6g~40t(~,6g)~%"
                     (select p-names i) (car coefs) (car se-s))
             (setf coefs (cdr coefs) se-s (cdr se-s)))
    (format t "~%")
    (format t "-2log-likelihood:~25t~13,6g~%" (* -2 (send self :loglike)))
    (format t "Number of cases:~25t~9d~%" (send self :num-cases))
    (format t "Number of records: ~25t~9d~%" (length (send self :keyint)))
    (format t "Tie handling: ~25t")
    (if (> (sum (send self :tieint)) 0) 
        (send (send self :tie-method) :prettyprint) 
      (format t "No ties present~%"))
 ))
 
(defmeth cox-regression-proto :display-with-formula 
  (&key (block-only *cox-display-block-only*))
"Arguments: &key (block-only *cox-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 :beta) 'list))
             (formula (send self :formula))
             (cov-mat (if (send self :agnostic) 
                          (send self :agnostic-covariance-matrix) 
                        (send self :covariance-matrix)))
             (nblocks (length (send formula :block-indices)))
             (block-indices (send formula :block-indices))
             (varnames (send self :predictor-names))
             )
         (format t (send self :comments))
         (if (> (length (unique (send self :strata))) 1) 
             (format t "~% Stratified ") 
           (format t "~%"))
         (format t "Maximum Partial Likelihood Estimates:~2%")
         (format t "~a~20t~a~35t~a~%" "Block" "Wald Chisq" "p-value")
         (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 "-2log-likelihood:~25t~13,6g~%" 
                 (* -2 (send self :loglike)))
         (format t "Number of cases:~25t~9d~%" 
                 (send self :num-cases))
         (format t "Number of records: ~25t~9d~%" 
                 (length (send self :keyint)))
         (format t "Tie handling: ~25t")
         (if (> (sum (send self :tieint)) 0) 
             (send (send self :tie-method) :prettyprint) 
           (format t "No ties present~%"))
         )
       )
      )
)

(defmeth cox-regression-proto :conf-interval 
  (&key (coverage 0.95) (names t))
"Arguments: &key  (coverage 0.95) (names t)
Produces a matrix of Wald confidence intervals for beta with predictor
names if requested  and  transforms them to the hazard ratio 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 (exp baseci))
       (name-list (send self :predictor-names))
       )
  (if names (bind-columns name-list xform) xform)
  )
)

(defmeth cox-regression-proto :cumulative-hazard-function (&key (xvec nil))
"Cumulative hazards evaluated at given covariate vector"
(let* (
       (hazards (send self :cumulative-hazard))
       (beta (send self :beta))
       (rr (exp (if xvec (inner-product xvec beta) 0)))
       )
  (map-elements #'(lambda (stratum) (list (first stratum) 
                                          (second stratum) 
                                          (* rr (third stratum))))
                hazards)
  )
)

(defmeth cox-regression-proto :survival-function 
  (&key (xvec nil) (product-limit t))
"Args: &key (xvec nil) (product-limit t)
Kaplan-Meier or near offer
Evaluated at specified covariate vector"
(let* (
       (hazards (send self :cumulative-hazard-function ))
       (beta (send self :beta))
       (rr (exp (if xvec (inner-product xvec beta) 0)))
       )
  (if product-limit
      (map-elements 
       #'(lambda (stratum)
           (let* (
                  (dcumhaz (difference (third stratum)))
                  (survival (^ (append-seq #(1) (cumprod (- 1 dcumhaz))) rr))
                  )
             (list (first stratum) (second stratum) survival)
             )
           )
       hazards
       )
    (map-elements
     #'(lambda (stratum) 
         (list (first stratum) (second stratum) 
               (exp (- (* rr (third stratum))))))
     hazards
     )
    )
  )
)

(defmeth cox-regression-proto :survival-function-closure 
  (&key (stratum 0 has-stratum) (xvec nil) (product-limit t))
"Arguments: &key (stratum 0 has-stratum)  (xvec nil) (product-limit t)
Function closure giving survival function for specified stratum and 
optionally covariate pattern"
(unless (or has-stratum (= 1 (length (unique (send self :strata))))) 
        (error "Must specify stratum"))
(let* (
       (kms (assoc stratum (coerce (send self :survival-function 
                                         :xvec xvec 
                                         :product-limit product-limit) 
                                   'list)
                   :test #'equalp
                   ))
       (ff (if kms  #'(lambda (ti) (interpolate ti (second kms) (third kms)))
             (error "Stratum ~a not found" stratum)))
       )
  ff
  )
)


(defmeth cox-regression-proto :cumulative-hazard-function-closure 
  (&key (stratum 0 has-stratum)  (xvec nil))
"Arguments: &key (stratum 0 has-stratum)  (xvec nil)
Function closure giving cumulative hazard function for specified 
stratum and optionally covariate pattern"
(unless (or has-stratum (= 1 (length (unique (send self :strata))))) 
        (error "Must specify stratum"))
(let* (
       (haz (assoc stratum 
                   (coerce (send self :cumulative-hazard-function 
                                 :xvec xvec) 
                           'list)
                   :test #'equalp
                   ))
       (ff (if haz  
               #'(lambda (ti) (interpolate ti (second haz) (third haz)))
             (error "Stratum ~a not found" stratum)))
       )
  ff
  )
)

(defmeth cox-regression-proto :plot-cumulative-hazard 
  (&key (one-plot t)(xvec nil))
"Arguments: (&key (one-plot t) (xvec nil)
Plots cumulative hazard functions, for baseline or for specified
covariate pattern"
(let* (
       (hazards (send self :cumulative-hazard-function :xvec xvec))
       )
  (cond ((and one-plot (> (length hazards) 1)) 
         (let* ((firstone (car hazards))
                (rest (cdr hazards))
                (plot (plot-steps (second firstone) (third firstone)  
                                  :variable-labels (list "time" "Hazard")))
                )
           (map-elements #'(lambda (stratum)
                     (add-steps (second stratum) (third stratum) plot)) rest)
           (send plot :adjust-to-data)
           (list plot)
           ))
        (t (map-elements #'(lambda (stratum)
                    (plot-steps (second stratum) (third stratum) 
                                :title (format nil "Stratum ~a" 
                                               (first stratum)) 
                                :variable-labels (list "time" "Hazard"))) 
                         hazards)
           )
        )
  )
)

(defmeth cox-regression-proto :plot-survival 
  (&key (inverted nil) (one-plot t) (xvec nil) (product-limit t))
"Arguments: &key (inverted nil) (one-plot t) (xvec nil) (product-limit t)
Plots survival function (or mortality function)
for each stratum at baseline or specified covariate pattern."
(let* (
       (survivals (send self :survival-function 
                        :xvec xvec 
                        :product-limit product-limit))
       )
  (cond ((and one-plot (> (length survivals) 1))
         (let* ((firstone (car survivals))
                (rest (cdr survivals))
                (plot (plot-steps (second firstone) 
                                  (if inverted 
                                      (- 1 (third firstone)) 
                                    (third firstone))  
                                  :variable-labels 
                                  (list "time"  
                                        (if inverted 
                                            "Proportion failing" 
                                          "Proportion surviving"))))
                )
           (map-elements #'(lambda (stratum)
                             (add-steps (second stratum) 
                                  (if inverted 
                                      (- 1 (third stratum)) 
                                    (third stratum)) 
                                  plot)) 
                         rest)
           (send plot :adjust-to-data)
           (if (null inverted) (send plot :range 1 0 1))
           (list plot)
           ))
        (t 
         (let (
               (plots (map-elements 
                       #'(lambda (stratum)
                           (plot-steps (second stratum) 
                                       (if inverted 
                                           (- 1 (third stratum)) 
                                         (third stratum)) 
                                       :title (format nil 
                                                      "Stratum ~a" 
                                                      (first stratum))
                                       :variable-labels 
                                       (list "time" (if inverted 
                                                        "Proportion failing" 
                                                      "Proportion surviving")))) 
                       survivals))
               )
           (unless inverted 
                   (map-elements 
                    #'(lambda (plot) (send plot :range 1 0 1))
                    plots))
           plots
           )
         )
        )
  )
)

(require "sliders")

(defmeth cox-regression-proto :slider-plot ()
"Arguments: 
Draws a survival plot with sliders"
(let* (
       (x (column-list (send self :x)))
       (names (send self :predictor-names))
       (ranges (mapcar #'(lambda (xi) 
                            (if (< (length (unique xi)) 10)
                                (sort-data (unique xi))
                                (list (min xi) (max xi)))
                          )
                x))
       (current-xvec (mapcar #'min x))
       (current-plot (elt (send self :plot-survival :one-plot t
                                :xvec current-xvec) 0))       
       (actions (mapcar 
                 #'(lambda (i) 
                     (format nil 
                             "(lambda (x) (adjust-to-sliders x ~s))" i))
                 (iseq (length x))))
       (lspec (mapcar #'(lambda (i) (list (elt names i) 
                                          (elt ranges i) 
                                          (elt actions i))) 
                      (iseq (length x))))
       )
  (defun adjust-to-sliders (x i) 
    (setf (select current-xvec i) x)
    (send current-plot :clear-points :draw nil)
    (send current-plot :clear-lines :draw nil)
    (map-elements #'(lambda (stratum)
                      (add-steps (second stratum) 
                                 (third stratum) 
                                 current-plot))  
                  (send self :survival-function :xvec current-xvec) )
    (send current-plot :range 1 0 1)
    )
  (make-and-run-sliders lspec)
  )
)

(defmeth cox-regression-proto :plot-group-survival-function 
  (&key new-x new-strata (latent-times nil) (product-limit t))
"Method arguments: &key new-x (latent-times nil) (product-limit t)
Plots expected survival for cohort with covariates given by :new-x, using 
the direct-adjusted method if :latent-times is nil and Hakulinen's method
otherwise."
(let* ((curve (send self :group-survival-function :new-x new-x 
                    :new-strata new-strata :latent-times latent-times 
                    :product-limit product-limit))
       (plot (plot-steps (first curve) 
                         (second curve) 
                         :title "Group Expected Survival" 
                         :variable-labels (list "time" "proportion surviving")
                         ))
       )
  (send plot :range 1 0 1)
  plot
)
)



(defmeth cox-regression-proto :group-survival-function 
  (&key new-x new-strata (latent-times nil) (product-limit t))
"Method arguments: &key new-x (product-limit t)
Returns expected survival for cohort with covariates given by :new-x, using 
the direct-adjusted method "
  (let* ((x (cond ((matrixp  new-x) (row-list  new-x))
                  ((objectp  new-x) (send new-x :design-matrix))
                  ((and (sequencep  new-x) (numberp (elt  new-x 0))) 
                   (row-list (bind-columns  new-x)))
                  ((and (sequencep  new-x) (sequencep (elt  new-x 0))) 
                   (row-list (apply #'bind-columns  new-x)))
             )
        ))
    (send self :exact-group-surv x new-strata product-limit))
  )



(defmeth cox-regression-proto :exact-group-surv (x strata product-limit)
  (let* ((beta (send self :beta))
         (stratlist (unique strata))
         (alltimes (sort-data (unique (send self :times))))
         (survs  (if strata 
                     (mapcar #'(lambda (stratum) 
                                 (send self :survival-function-closure 
                                       :stratum stratum 
                                       :product-limit product-limit)
                                 )
                             stratlist
                             )
                   (list (send self :survival-function-closure 
                               :product-limit product-limit))
                   ))
         (rr (mapcar #'(lambda (xi) (exp (inner-product xi beta))) x))
         (rrstrat (if strata 
                      (map-elements #'(lambda (stratum)  
                                        (select rr (argselect strata stratum)))
                                    stratlist)
                    (list rr)))
         (ff (lambda (rrs sts) 
               (map-elements #'(lambda (sti) (sum (^ sti rrs))) sts)))
         (basesurvs (mapcar 
                     #'(lambda (fstrat) (funcall fstrat alltimes)) 
                     survs))
         (survcurvs (mapcar 
                     #'(lambda (rrs sts) (funcall ff rrs sts))  
                     rrstrat basesurvs))
         )
    (list alltimes  (/ (apply #'+ survcurvs) (length x)))
    )
  )

(defmeth cox-regression-proto :martingales (&optional (new nil set))
"Message ars: (&optional new)
Sets or returns martingale residuals"
  (when set
        (setf (slot-value 'martingales) new)
        )
  (slot-value 'martingales))


(defmeth cox-regression-proto :scaled-schoenfeld-residuals ()
  (let* ((ri (send self :schoenfeld-residuals))
         (n (length ri))
         (betas (matmult (bind-columns (repeat 1 n)) 
                         (bind-rows (send self :beta))))
         (scaled  (+ betas 
                     (matmult  (apply #'bind-rows (mapcar #'second ri)) 
                               (* n (send self :covariance-matrix)))))
         )
    (cons  (mapcar #'first ri) (column-list scaled))
    )
  )


(defmeth cox-regression-proto :plot-ph-test (&optional vars)
"Arguments: vars
Plots smoothed scaled Schoenfeld residuals"
(let* (
       (r (send self :scaled-schoenfeld-residuals))
       (vns (send self :predictor-names))
       (varlist (or vars (iseq (length (send self :beta)))))
       )
  (map-elements 
   #'(lambda (v) 
       (scatter-smooth (first r) (elt r (+ 1 v)) 
                       :title 
                       (format nil "Proportional hazards: ~a" (elt vns v)) 
                       :variable-labels 
                       (list "Time" "Scaled Schoenfeld Residuals"))) varlist)
  )
)

(defmeth cox-regression-proto :harrell-lee-test ()
"Args:
Harrell & Lee's test for proportional hazards 
based on the Schoenfeld residuals"
  (second (send self :ph-test :scaled nil :unscaled #'rank))
) 

(defmeth cox-regression-proto :ph-test (&key (scaled #'identity) (unscaled nil))
"Args: &key (scaled #'identity) (unscaled nil)
Tests for proportional hazards based on scaled or raw Schoenfeld residuals"
(list (cond ((null scaled) nil)
            ((functionp scaled)
             (let* ((sr (send self :scaled-schoenfeld-residuals))
                    (srs (cons (funcall scaled (first sr)) (cdr sr)))
                    (cmat (apply #'covariance-matrix srs))
                    (sds (/ (sqrt (diagonal cmat))))
                    (rmat (matmult (diagonal sds) cmat (diagonal sds)))
                    (rhos (cdr (coerce (first (row-list rmat)) 'list)))
                    (n2 (- (length (first srs)) 2))
                    (ts (* rhos (sqrt n2)))
                    (names (send self :predictor-names))
                    (ps (* 2 (t-cdf (- (abs ts)) n2)))
                    )
               (bind-columns names rhos ts ps)
               )
             )
            )
      (cond ((null unscaled) nil)
            ((functionp unscaled)
             (let* ((sr (send self :schoenfeld-residuals))
                    (srs (cons (funcall unscaled (firsts sr)) 
                               (column-list (apply #'bind-rows (seconds sr)))))
                    (cmat (apply #'covariance-matrix srs))
                    (sds (/ (sqrt (diagonal cmat))))
                    (rmat (matmult (diagonal sds) cmat (diagonal sds)))
                    (rhos (cdr (coerce (first (row-list rmat)) 'list)))
                    (n2 (- (length (first srs)) 2))
                    (ts (* rhos (sqrt n2)))
                    (names (send self :predictor-names))
                    (ps (* 2 (t-cdf (- (abs ts)) n2))) 
                    )
               (bind-columns names rhos ts ps)
               )
             )
            )
      )
)
(defmeth cox-regression-proto :delta-betas (&key (total t))
"Args: &key (total t)
Returns a list of sequences containing delta-betas.  If total is t the delta-betas are
collapsed over records with the same id and the first sequence is the unique ids; if 
total is nil the ordering is the same as in the data"
(let* ((id (send self :id))
       (unique-id (unique id))
       (scores (or (send self :scores) 
                   (error "Can't calculate delta-betas - no martingale residuals")))
       (invinfomat (send self :model-covariance-matrix))
       (partial-db (row-list (matmult scores invinfomat )))
       (db  (cond 
             (total 
              (mapcar #'(lambda (i) 
                          (apply #'+ (select partial-db (argselect id i)))) 
                      unique-id)) 
             (t 
              partial-db)))
       )
  (if total  
      (cons unique-id (column-list (apply #'bind-rows db)))
    (column-list (apply #'bind-rows db))
    )
  )
)

(defmeth cox-regression-proto :plot-delta-betas (&key (total t) (link t) (vars nil))
"Args: &key (total t) (link t) (vars nil)
Plots delta-betas in linked plots labelled by id. A subset of the variables 
can be chosen with vars.
The default is to sum the delta-betas over all records with 
the same id, set total to nil to get separate delta-betas for different records."

(let* ((db (send self :delta-betas :total total))
       (idlabel (map-elements #'write-to-string 
                              (if total 
                                  (first db) 
                                (send self :id))))
       (titles (send self :predictor-names))
       (ylabels (mapcar #'(lambda (title) 
                            (format nil "dBeta for ~a" title)) titles))
       (rval (mapcar #'(lambda (db title ylabel) 
                         (index-plot db 
                                     :title title 
                                     :point-labels idlabel 
                                     :variable-labels (list "index" ylabel))) 
                     (if total 
                         (cdr db) 
                       db) 
                     titles 
                     ylabels))
       )
  (when link (apply #'link-views rval))
  rval
))


(defmeth cox-regression-proto :plot-new-variable (newx &optional label)
"Args: newx &optional label
 Smooth plot of variable newx against martingale residuals, where newx should be a
 variable not in the model.  The smooth estimates the correct functional form of newx"
(let* ((rmart (send self :martingales))
       (xlabel (or label "New covariate"))
       )
  (scatter-smooth newx 
                  rmart 
                  :variable-labels (list xlabel "Martingale residuals") 
                  :point-labels (mapcar #'write-to-string (send self :id)))
  )
)


(require "expsurv1")
(defmeth  cox-regression-proto :scatterplot-expected 
  (&key covariates covariate-labels)
"Args: &key covariates covariate-labels
Scatterplot matrix of covariates linked to expected & K-M survival graph"
  (let* ((scat (scatterplot-matrix covariates 
                                   :variable-labels covariate-labels))
         (time (send self :times))
         (status (send self :status))
         (coxobject self)
         (data (km-plot time status (kmest time status)))
         (km (plot-lines (first data) 
                         (second data) 
                         :title "Survival Plot"))
         (current-points (iseq (length time))))
    (send km :range 1 0 1)
    (send scat :point-state (iseq (length time)) 'hilited)
    (defmeth scat :unselect-all-points ()
      (setf current-points nil)
      (send km :clear)
      (call-next-method))
    (defmeth scat :adjust-points-in-rect (x1 y1 w h s)
      (let* ((p-i-r (send self :points-in-rect x1 y1 w h)))
        (setf current-points
              (case (send self :mouse-mode)
                (brushing p-i-r) 
                (selecting (union p-i-r current-points))))
        (if current-points
            (prog* 
             ((points (sort-data current-points))
              (time-sel (select time points))
              (status-sel (select status points))
              (expected (send coxobject 
                              :group-survival-function 
                              :new-x  (mapcar 
                                       #'(lambda (x) 
                                           (select x points))
                                       (column-list (send coxobject :x))) 
                              :new-strata 
                              (select (send coxobject :strata) points)))
              (new-plot (km-plot time-sel status-sel 
                                 (kmest time-sel status-sel))))
             (send km :clear :draw nil)
             (send km :add-lines (first expected) (second expected))
             (send km :add-lines 
                   (first new-plot) (second new-plot))))
        (call-next-method x1 y1 w h s)))
    scat))

(defmeth  cox-regression-proto :scatterplot-predicted 
  (&key (covariates nil) covariate-labels  x times status (strata nil))
"Args: &key (covariates nil)  x times status (strata nil)
Scatterplot matrix of covariates linked to predicted & K-M survival graph"
  (let* ((scat (scatterplot-matrix (or covariates x) 
                                   :variable-labels covariate-labels))
         (time times)
         (coxobject self)
         (strata (or strata (* times 0)))
         (data (km-plot time status (kmest time status)))
         (km (plot-lines (first data) 
                         (second data) 
                         :title "Survival Plot"))
         (current-points (iseq (length time))))
    (send km :range 1 0 1)
    (send scat :point-state (iseq (length time)) 'hilited)
    (defmeth scat :unselect-all-points ()
      (setf current-points nil)
      (send km :clear)
      (call-next-method))
    (defmeth scat :adjust-points-in-rect (x1 y1 w h s)
      (let* ((p-i-r (send self :points-in-rect x1 y1 w h)))
        (setf current-points
              (case (send self :mouse-mode)
                (brushing p-i-r) 
                (selecting (union p-i-r current-points))))
        (if current-points
            (prog* 
             ((points (sort-data current-points))
              (time-sel (select time points))
              (status-sel (select status points))
              (expected (send coxobject :group-survival-function 
                              :new-x   (mapcar 
                                        #'(lambda (xi) 
                                            (select xi points))
                                        x) 
                              :new-strata (select strata points)))
              (new-plot (km-plot time-sel status-sel 
                                 (kmest time-sel status-sel))))
             (send km :clear :draw nil)
             (send km :add-lines (first expected) 
                   (second expected))
             (send km :add-lines 
                   (first new-plot) (second new-plot))))
        (call-next-method x1 y1 w h s)))
    scat))



(defmeth cox-regression-proto :isnew (&key x times status
                                           entry
                                           id
                                           (formula nil)
                                           (ties *cox-default-ties*)
                                           pweights
                                           (verbose t)
                                           predictor-names
                                           response-name
                                           (recycle nil)
                                           case-labels
                                           (init-beta nil)
                                           count-limit
                                           (print t) 
                                           print-summary
                                           statusint 
                                           keyint
                                           agnostic
                                           entryint
                                           strata
                                           tieint
                                           comments
                                           (suppress-warnings nil)
                                           (compute t)
                                           (time-dependent nil)
                                           )
  (send self :x x)
  (send self :time-dependent time-dependent)
  (send self :times times)
  (send self :status status)
  (send self :entry entry)
  (send self :tie-method ties)
  (send self :weights (or pweights (repeat 1 (length times))))
  (send self :formula formula)
  (send self :strata strata)
  (send self :recycle recycle)
  (cond (init-beta (send self :set-beta init-beta) (send self :recycle t)))
  (send self :verbose verbose)
  (send self :agnostic agnostic)
  (send self :suppress-warnings suppress-warnings)
  (when predictor-names (send self :predictor-names predictor-names))
  (when response-name (send self :response-name response-name))
  (send self :epsilon *cox-tolerance*) 
  (send self :epsilon-dev  *cox-tolerance*)
  (send self :count-limit (or count-limit *cox-count-limit*))
  (send self :statusint statusint)
  (send self :keyint keyint)
  (send self :tieint tieint)
  (send self :entryint entryint)
  (send self :comments comments)
  (send self :id id)
  (cond (print (send self :display))
        (print-summary (send self :display-with-formula))
        (compute (send self :compute))
        )
)
(provide "coxreg")
(defun cox-model (&key x times status
                       (id (iseq (length times)) has-id)
                       (entry nil)
                       (strata (* 0 status) has-strata)
                       (agnostic nil)
                       (verbose *cox-verbose*)
                       (print t)
                       (init-beta nil)
                       (count-limit *cox-count-limit*)
                       (weights nil)
                       (predictor-names nil has-names)
                       (response-name "Failure time" has-response-name)
                       (ties *cox-default-ties*)
                       (suppress-warnings nil)
                       time-x
                       time-funs
                       )
"Constructs a Cox regression modelling with right censoring and 
optionally left truncation and strata."
(unless (and x times status) 
        (error "You must supply :x :times and :status.~%"))
(let* ((n (length times))
       (has-unique-id (= (length (remove-duplicates id)) (length id)))
       (comments (cond 
                  ((and agnostic (not has-unique-id)) 
                   "Cox model with LWA/WLW agnostic covariance estimator")
                  ((and agnostic has-unique-id) 
                   "Cox model with LW agnostic covariance estimator")
                  (t 
                   "Standard Cox regression")
                  )
                 )
       (alldata (if entry  
                    (col2row (list (append-seq entry times)  
                                   (append-seq (repeat 0 n) status) 
                                   (append-seq (repeat -1 n) (repeat 1 n)) 
                                   (append-seq strata strata) (repeat (iseq n) 2)))
                  (col2row (list times  status  (repeat 1 n) strata (iseq n)))
                  ))
       (rtemp  (sort (copy-list alldata) #'after))
       (temp (row2col rtemp))
       (keyint (select temp 4))
       (statusint (select temp 1))
       (entryint (select temp 2))
       (ntimes (length keyint))
       (timexlist (cond
                   ((null time-x) 
                    nil)
                   ((matrixp time-x) 
                       (row-list time-x))
                   ((and (sequencep time-x) (numberp (elt time-x 0))) 
                    (row-list (bind-columns time-x)))
                   ((and (sequencep time-x) (sequencep (elt time-x 0))) 
                    (row-list (apply #'bind-columns time-x)))
                   (t 
                    (error "Wrong sort of :time-x"))
                   ))
       ;; tieint =1 if this entry is tied with the previous one, else 0
       (tieint (append  
                (mapcar #'tied (select rtemp (iseq 0 (- ntimes 2))) 
                        (select rtemp (iseq 1 (- ntimes 1)))) 
                   '(0)))
       )
  (cond ((objectp x)
         (send cox-regression-proto :new :x (send x :design-matrix)
               :times times :status status :entry entry
               :predictor-names 
               (if has-names predictor-names (send x :name-list))
               :init-beta init-beta
               :count-limit count-limit
               :weights weights
               :formula x
               :print nil
               :print-summary print
               :agnostic agnostic
               :verbose verbose
               :keyint keyint
               :strata strata
               :statusint statusint
               :entryint entryint
               :tieint tieint
               :ties ties
               :reponse-name response-name
               :comments comments
               :suppress-warnings suppress-warnings
               :time-dependent (list time-funs timexlist)
               :id id
               ))
        (t 
         (send cox-regression-proto :new :x x
               :times times :status status :entry entry
               :predictor-names predictor-names
               :init-beta init-beta
               :count-limit count-limit
               :weights weights
               :formula x
               :print print
               :print-summary nil
               :agnostic agnostic
               :verbose verbose
               :strata strata
               :keyint keyint
               :statusint statusint
               :entryint entryint
               :tieint tieint
               :ties ties
               :reponse-name response-name
               :comments comments
               :suppress-warnings suppress-warnings
               :time-dependent (list time-funs timexlist)
                  :id id
                  )  
         )
        )
  )

)
(defun kaplan-meier 
(&key times status (entry nil) (groups (repeat 0 (length times))))
"Args: &key times status groups entry"
(let* ((fakex (uniform-rand (length times)))
       (nullcox (cox-model :x fakex :times times :entry entry :status status 
                           :strata groups :verbose nil :print nil :print-summary nil 
                           :count-limit 1 :suppress-warnings t))
       )
   (send nullcox :compute)
   (send nullcox :survival-function)
  )
)

(defun plot-kaplan-meier 
  (&key times status (entry nil) (groups  (repeat 0 (length times))) inverted (one-plot t))
"Args: &key times status groups entry inverted (one-plot t)"
(let* ((fakex (uniform-rand (length times)))
       (nullcox (cox-model :x fakex :times times :entry entry :status status 
                           :strata groups :verbose nil :print nil :print-summary nil 
                           :count-limit 1 :suppress-warnings t))
       )
   (send nullcox :compute)
   (send nullcox :plot-survival :inverted inverted :one-plot one-plot)
   )
)

(defun nelson-aalen 
  (&key times status (entry nil)(groups  (repeat 0 (length times))))
"Args: &key times status groups entry "
(let* ((fakex (uniform-rand (length times)))
       (nullcox (cox-model :x fakex :times times :entry entry :status status 
                           :strata groups :verbose nil :print nil :print-summary nil 
                           :count-limit 1 :suppress-warnings t))
       )
  (send nullcox :compute)
  (send nullcox :cumulative-hazard)
  )
)

(defun plot-nelson-aalen 
  (&key times status (entry nil) (groups  (repeat 0 (length times))) (one-plot t))
"Args: &key times status groups entry (one-plot t)"
(let* ((fakex (uniform-rand (length times)))
       (nullcox (cox-model :x fakex :times times :entry entry :status status 
                           :strata groups :verbose nil :print nil :print-summary nil 
                           :count-limit 1 :suppress-warnings t))
       )
  (send nullcox :compute)
  (send nullcox :plot-cumulative-hazard :one-plot one-plot)
  )
)

(defun cox-score-test 
  (&key groups times status (entry nil) (trend nil) (strata (repeat 0 (length times))))
"Args: &key groups times status entry (trend nil)
Peto's approximate Logrank test for differences between groups.
Requires modelformula extensions."
(let* (
       (gmodel (factor groups :namebase "Group "))
       (gmat (first gmodel))
       (gnames (second gmodel))
       (p (second (array-dimensions gmat)))
       (x (if trend (matmult gmat (iseq 1 p)) gmat))
       (nullcox (cox-model :x x :times times :entry entry :status status :strata strata 
                           :verbose nil :print nil :print-summary nil :count-limit 1
                           :suppress-warnings t :ties peto-ties))
       )
  (send nullcox :compute)
  (let* (
         (beta (send nullcox :beta))
         (infomat (inverse (send nullcox :covariance-matrix)))
         (score (matmult beta infomat))
         (grouptests (/ score (sqrt (diagonal infomat))))
         (overalltest (inner-product score beta))
         )
    (format t (if trend "Logrank Trend Test~%" "Logrank test~%"))
    (when (and (> p 1) (null trend))
          (dotimes (i p)
                   (format t "Group ~a.  Z=~6,4g   p=~6,4g~%" 
                           (elt gnames i) 
                           (elt grouptests i) 
                           (* 2 (- 1 (normal-cdf (abs (elt grouptests i))))))
                   )
          )
    (format t "Overall test: X^2=~8,6g   p=~6,4g~%" 
            overalltest 
            (- 1 (chisq-cdf overalltest p)))
    (if (< (length (unique times)) (length times)) 
        (format t "Note: Tied times present -- results slightly conservative~%"))
    (list grouptests overalltest)
    )
  )
)
(def *cox-count-limit* 20)
(def *cox-default-ties* efron-ties)
(def *cox-verbose* t)
(def *cox-tolerance* (sqrt (sqrt machine-epsilon)))
(def *cox-display-block-only* nil)
(def *cox-paranoia* t)