(require "coxreg")
(load "pbc.lsp")
(wait)
(def mayo-model-formula 
     (as-formula '((term age) (term (log bili)) (term (log protime)) 
                   (term (log alb)) (term edtrt))))

(def mayo-model (cox-model :times time :status status :x mayo-model-formula))
(wait)
(def mayo-model-without-formula 
     (cox-model :times time :status status 
                :x (list age (log bili) (log protime) (log alb) edtrt)))

(wait)
(def stage-model (cox-model :times time :status status 
                            :x (as-formula '((factor stage)))))

(wait)
(def mayo-model-stage 
     (cox-model :times time :status status 
                :x (as-formula '((term age) (term (log bili)) 
                                 (term (log protime)) (term (log alb)) 
                                 (term edtrt) (factor stage)
                                 ))))

(- 1 (chisq-cdf 
      (* 2 (- (send mayo-model-stage :loglike) 
              (send mayo-model :loglike)))
      3
      )
   )

(def mayo-model-stratified 
     (cox-model :times time :status status :strata edtrt
                :x (as-formula '((term age) (term (log bili)) 
                                 (term (log protime)) (term (log alb)) 
                                 ))))
(wait)
(print-matrix (send mayo-model :conf-interval))
(wait)
(send mayo-model :plot-survival)
(wait)
(def typical (list 0 (log 3.5) (log 10) (log 1.75 ) 50)) 
(send mayo-model :plot-survival :xvec typical)
(wait)
(def typical2 (list (log 3.5) (log 10) (log 1.75 ) 50)) 
(send mayo-model-stratified :plot-survival :xvec typical2 :one-plot t)
(wait)

(def sfun (send mayo-model :survival-function-closure :xvec typical))
(funcall sfun (* 365.25 (iseq 1 10)))
(wait)

(send mayo-model :slider-plot)
(wait)
(load "pbc-new.lsp")
(def risks (matmult 
            (bind-columns edtrt1 (log alb1) (log protime1)
                          (log bili1) age1) 
            (send mayo-model :beta)))

(def low (which (> 5.5 risks)))
(def high (which (< 6.5 risks)))
(def med (which 
          (map-elements #'(lambda (x y) 
                            (and x y)) 
                        (> 6.5 risks) (< 5.5 risks)
                        )
          )
     )
(def lowrisk (send mayo-model :group-survival-function 
                   :new-x (mapcar #'(lambda (l) 
                                      (select l low)) 
                                  (list edtrt1 (log alb1) 
                                        (log protime1) (log bili1) 
                                        age1
                                        )
                                  )
                   )
     )
(def medrisk (send mayo-model :group-survival-function
                   :new-x (mapcar #'(lambda (l) 
                                      (select l med)) 
                                  (list edtrt1 (log alb1)
                                        (log protime1) (log bili1)
                                        age1
                                        )
                                  )
                   )
     )
(def hirisk (send mayo-model :group-survival-function 
                  :new-x (mapcar #'(lambda (l)
                                     (select l high)) 
                                 (list edtrt1 (log alb1)
                                       (log protime1) (log bili1)
                                       age1
                                       )
                                 )
                  )
     )

(def gtime1 (mapcar #'(lambda (i)
                        (select time1 i)) 
                    (list low med high)
                    )
     )
(def gstat1 (mapcar #'(lambda (i) 
                        (select status1 i))
                    (list low med high)
                    )
     )
(def kmplots (mapcar #'(lambda (tt ss)
                         (km-plot tt ss (kmest tt ss)))
                     gtime1 
                     gstat1
                     )
     )

(def validation 
     (plot-lines (elt kmplots 0)
                 :type 'D 
                 :title "Predictive accuracy of Mayo Model" 
                 :variable-labels (list "Days on study"
                                        "Proportion surviving")
                 )
     )
(send validation :add-lines (elt kmplots 1) :type 'D)
(send validation :add-lines (elt kmplots 2) :type 'D)
(send validation :range 1 0 1)
(send validation :add-lines lowrisk :color (elt *colors* 2))
(send validation :add-lines medrisk :color (elt *colors* 2))
(send validation :add-lines hirisk :color (elt *colors* 2))
(wait)
(def mayo-ph (send mayo-model :ph-test :scaled #'log))
(print-matrix (first mayo-ph))

(print-matrix (send mayo-model :harrell-lee-test))
(wait)

(send mayo-model :plot-ph-test)
(wait)

(def mayo-formula-no-bili (as-formula '((term age) (term (log protime)) 
                                        (term (log alb)) (term edtrt))))

(def mayo-no-bili (cox-model  :times time :status status 
                              :x mayo-formula-no-bili))

(send mayo-no-bili :plot-new-variable bili "Bilirubin")
(wait)

(send mayo-no-bili :plot-new-variable (log bili) "log(Bilirubin)")
(wait)

(send mayo-model :plot-delta-betas)
(wait)

(send mayo-model :scatterplot-expected :covariates (list stage copper sgot)
                :covariate-labels (list "Stage" "urine copper" "SGOT"))
(wait)

(load "pbc-new")
(send mayo-model :scatterplot-predicted 
      :x (list edtrt1 (log alb1) (log protime1) (log bili1) age1)
      :times time1
      :status status1
      :covariate-labels (list "edema" "log alb" "log protime" "log bili" "age")
      )



(wait)
(cox-model :times time :status status :x mayo-model-formula :ties efron-ties)
(cox-model :times time :status status :x mayo-model-formula :ties breslow-ties)
(def tt (- time (* 0.01 status (uniform-rand 312))))
(cox-model :times tt :status status :x mayo-model-formula)
(wait)
(cox-score-test :groups trt :times time :status status)
(wait)
(plot-kaplan-meier :groups trt :times time :status status)
(def id (iseq 312))
(def after2 (select id (which (>= time 1096))))
(def protime2 (append (repeat 1 312) (select protime after2)))
(def protime1 (append protime (repeat 1 (length after2))))
(def edtrt12 (append edtrt (select edtrt after2)))
(def alb12 (append alb (select alb after2)))
(def age12 (append age (select age after2)))
(def bili12 (append bili (select bili after2)))
(def id12 (append id after2))
(def time12 (append (mapcar #'(lambda (ti) (min ti 1096)) time) 
                   (select time after2)))
(def entry12 (append (repeat 0 312) (repeat 1096 (length after2))))
(def status12 (append (mapcar #'(lambda (ti si) 
                                  (if (< ti 1096) si 0)) time status)
                      (select status after2)))


(def mayo-model-3-years 
     (cox-model :times time12 :status status12 :entry entry12
                :x (as-formula '((term (log protime1)) (term (log protime2))
                                 (term edtrt12)
                                 (term (log alb12)) (term age12)
                                 (term (log bili12))))
                ))

(def mayo-model-timedep
     (cox-model :times time :status status 
                :x (list (log alb) (log bili) edtrt (log protime))
                :time-x (list protime age) 
                :time-funs (list #'(lambda (xi ti)
                                     (* (elt xi 0)  (/ ti 1000)))
                                 #'(lambda (xi ti)
                                     (+ (elt xi 1) ti))
                                 )
                :predictor-names (list "log alb" "log bili" "edema" 
                                       "log protime" "(log protime)(t)"
                                       "attained age")
                ))

(wait)
(load "cgd")

(def cgd-formula-1 (as-formula '((factor cgd-trt))))
(def cgd-formula-2 (as-formula '((factor cgd-trt) (term cgd-age) 
                                 (factor cgd-type) (factor cgd-sex) 
                                 )))

(def cgd-first (pmin cgd-seqnum 2))

(def cgd-formula-1-strat (as-formula '((factor cgd-trt)
                                       (interaction (list cgd-trt cgd-first))
                                       )))

(def cgd-formula-2-strat (as-formula '((factor cgd-trt) (term cgd-age)
                                       (factor cgd-type) (factor cgd-sex)
                                       (interaction (list cgd-trt cgd-first))
                                       (interaction (list cgd-type cgd-first))
                                       )))



(wait)
(def andersen-gill-1 
     (cox-model :times cgd-stop :entry cgd-start :status cgd-status
                :x cgd-formula-1
                :id cgd-id :agnostic t))

(def andersen-gill-2
     (cox-model :times cgd-stop :entry cgd-start :status cgd-status
                :x cgd-formula-2
                :id cgd-id :agnostic t))
(wait)

(def conditional-1
     (cox-model :times cgd-stop :entry cgd-start :status cgd-status
                :x cgd-formula-1 :strata cgd-seqnum
                :id cgd-id :agnostic t))
(def conditional-2
     (cox-model :times cgd-stop :entry cgd-start :status cgd-status
                :x cgd-formula-2 :strata cgd-seqnum
                :id cgd-id :agnostic t))