;;;
;;; Examples used in the documentation
;;;

;; load the data
(load "karim.lsp")

;; display the variables
(variables)

;; fit the first model
(def example1 (gee-model :y num-visits :x (list sex smoke age time2 time3 time4) :g id :error poisson-error :correlation (stat-m-dependence 3) :times times :predictor-names (list "sex" "smoker" "mother's age" "2nd vs 1st" "3rd vs 1st" "4th vs 1st")))

;; Print out a table with the hazard ratios and 90% confidence intervals
;; (the coefficients are log hazard ratios)

(print-matrix (send example1 :conf-interval :transform #'exp :coverage 0.9))

;;
;; Same model, using model formulae
(def karim-model (as-formula '((term age) (factor gender) (factor smoker) (factor times))))

(def example2 (gee-model :y num-visits :x karim-model :g id :error poisson-error :correlation (stat-m-dependence 3) :times times))


;; plot delta-betas for deletion of individual cases
;; Observations 15 and 12 highly influential on smoker
;; (and on third time variable)
(send example1 :plot-dbeta)

;; a plot-linking example for completeness
(def resid (send example1 :deviance-residuals)) 
(def qqnorm (quantile-plot resid))
(send example1 :bind-plot-to-gee qqnorm +observation+)

;; now plot delta-betas for deletion of observations, only for 
;; smoker and the third time variable
(send example1 :plot-dbeta :unit +group+ :variable (list 2 6))

;; Group 19 is influential: it's the same observations
(select id #(12 15))

;; Look at Pearson residuals against linear predictor
(send example1 :plot-pearson-residuals :variable (list 0))

;; and standardised Pearson residuals
(send example1 :plot-standardised-residuals :variable (list 0))



;; Now remove group 19 and see how the fit changes
;;
(def not19 (remove-if #'(lambda (i) (= (elt id i) 19)) (iseq 292)))
(def newy (select num-visits not19))
(def newsmoke (select smoke not19))
(def newsex (select sex not19))
(def newtime2 (select time2 not19))
(def newtime3 (select time3 not19))
(def newtime4 (select time4 not19))
(def newage (select age not19))
(def newid (select id not19))
(def newtimes (select times not19))

(def example3 (gee-model :y newy :x (list newsex newsmoke newage newtime2 newtime3 newtime4) :g newid :error poisson-error :correlation (stat-m-dependence 3) :times newtimes :predictor-names (list "sex" "smoker" "mother's age" "2nd vs 1st" "3rd vs 1st" "4th vs 1st")))

;; Perhaps there's something strange about the correlation structure
;; Try saturated correlation

(def example4 (gee-model :y newy :x (list newsex newsmoke newage newtime2 newtime3 newtime4) :g newid :error poisson-error :correlation saturated-corr :times newtimes :predictor-names (list "sex" "smoker" "mother's age" "2nd vs 1st" "3rd vs 1st" "4th vs 1st")))