#
# association analysis for categorical traits with covariates
#
eval (define (paste x) \
  (if (atom? x) \
    x \
    (string-append (car x) \
      (apply string-append \
        (map (lambda (x) (string-append " " x)) (cdr x))))))
#
# construct model terms
#
eval (define (twoway alist) \
  (define (pairoff x) \
   (let loop ((A (car x)) (B (cdr x)) (res '())) \
    (if (> (length B) 1) \
       (loop (car B) (cdr B) (append res (map (lambda (x) (list A x)) B))) \
       (append res (list (cons A B)))))) \
  (define (convert-pairoff x) \
    (map (lambda (y) (string-append (car y) "*" (cadr y))) x)) \
  (convert-pairoff (pairoff alist)))
#
# mul <trait> <gene> <cov1> <cov2> ...
#
eval (define (s-multinomial-1 trait marker covariates model-type) \
        (let ((args (append (list trait) (list marker) covariates))) \
          (run (string-append "llm " (paste args) " " \
                 (paste (twoway args)) " " model-type "; " \
                 "llm " (paste args) " " \
                 (paste (cdr (twoway args))) " " model-type)) \
          (let* ((s1 (stat-result "lrt")) \
                 (df1 (stat-result "df"))  \
                 (p (pchisq s1 df1))) \
            (locstat-set! marker p) \
            (format "~12A ~8,2F ~3D ~8,6F~\%" marker s1 df1 p))))
                   
#
# multinomial <trait> <cov1> <cov2>...
#
eval (define (s-multinomial trait covariates) \
  (if (string=? (loctyp trait) "c") \
     (let ((n (nloci))  \
           (cov (if (pair? covariates) covariates (string-split covariates)))) \
       (format "     Trait: ~S~\%" trait) \
       (format "Covariates: ~S ~\%~\%" (paste cov)) \
       (format "~A~\%" "Marker       Chi-sq    df  P") \
       (format "~A~\%" "------------ -------- --- -------- ") \
       (let loop ((i 1) (nmarkers 0)) \
         (if (< i n) \
            (begin \
              (if (string=? (loctyp i) "m") \
                (begin \
                  (s-multinomial-1 trait (loc i) cov "allelic") \
                  (loop (+ i 1) (+ nmarkers 1))) \
                (loop (+ i 1) nmarkers))) \
            (format "~\%Tested ~D markers~\%" nmarkers)))) \
     (format "~a~\%" "ERROR: Expected categorical trait and covariates."))) 
#
# front end
#
macro multinomial
 macro ple <- plevel
 set ple -2 silent 
 eval (s-multinomial "%1" "%+2")
 set ple %ple silent
;;;;
