;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;
;;;;;;  Permutation tests for any method in the inverse-regression-model-proto.
;;;;;;  Written by Dennis Cook.
;;;;;;  Version 0.98, May 15, 2000
;;;;;;
;;;;;;  After loading the code and setting up an inverse regression, sir, save,
;;;;;;  phd, ..., a new Permutation Tests item appears in the model menu.
;;;;;;  Selecting that item and supplying a number of permutations results in
;;;;;;  permutation tests for dimension as output.  The tests are used in the
;;;;;;  same way as the tests for sir or phd.  Computations may take several minutes
;;;;;;  if the data set or the number of permutations is large.
;;;;;;
;;;;;;  Seems to work OK with deleted observations and/or weights.  At present
;;;;;;  the response is always permuted even if using pHd with residuals.
;;;;;;
;;;;;;  Reference:  Cook, R. D. and Yin, X. (2000).  Dimension Reduction and
;;;;;;  Visualization in Discriminant Analysis.  Australia/New Zealand Journal
;;;;;;  of Statistics, to appear.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;



(rc-menu-item 'permutation-test-menu-item  "Permutation tests"
        :perm-test)

(defparameter *inverse-regression-menu-items*
  (list 'display-sir-menu-item
   'plot-sir-directions-menu-item
   'permutation-test-menu-item
   'remove-model-menu-item))



(defmeth inverse-regression-model-proto :perm-test-pvals (B &key (num-lc
*ir-num-lc*))
"Method args: ((B &key (num-lc *ir-num-lc*))
This is a general permutation test for partial sums of eigenvalues.  This
version
records the counts in the slot 'count for further study. B is the number of
permutations. The pvalues are stored in slot 'perm-pvals."

	(let* (
        (method (send self :method))
        (numpred (length (send self :basis)))
        (ntemp (min numpred num-lc))
;;;
;;; In sir numtests, the number of permutation tests, cannot exceed
;;;  the number of slices minus 1
;;;
        (numtests (if (equal method :sir)
                      (min ntemp (- (send self :actual-nslices) 1))
                      ntemp))
;;;
;;; Store original data set
;;;
        (inc (which (send self :included)))
        (n (send self :num-cases))
        (ninc (send self :num-included))
        (evs (send self :eigenvalues))
        (ux0 (send self :x))
        (ux (send self :dirs))
        (uy (send self :y))
        (wts (send self :weights))
;;;
;;; Initialize lists of counts and p-values.  Power is the power to which the eigenvalues
;;; are raised before adding. Power=1 for SIR and SAVE but power=2 for pHd and perhaps
;;; future methods.  power=1 by default.  To set power=2 for a method and the method in
;;; the case function below.
;;;
        (count (repeat (list (repeat nil 1)) numtests))
        (pvals (repeat nil numtests))
        (power 1)
        pan
        )
   (case method
      ((:PHD :PHDY)
       (setf power 2))
      (t))
;;;
;;; A friendly message
;;;
   (format t "Calculating, please be patient...~%")
;;;
;;;  Set weights to 1 for use in testing
;;;
;;;  (send self :weights nil)
;;;
;;; Start the main permutation loop. B = no of permutations.  Only the included
;;; cases are permutated in the second statement of the main loop.  Pan is the list of
;;; permuted indices.
;;;
   (dotimes (i B)
            (setf pan (iseq n))
            (setf (select pan inc) (sample inc ninc))
            (send self :x ux)
            (send self :y (select uy pan))
            (setf (nth 0 count) (append (nth 0 count)
                                        (list (sum (^ (send self
:eigenvalues) power)))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Weights if pesent are assumed to be a function of X only, as would be
the case when
;;; using weights from the Reweight-for-Symmetry command.  Under the null
;;; hypothesis of d dimensions, the weights stay with the permuted
;;; d predictors (X_d) since the weighting doesn't matter for the last p-d
;;; unpermuted perdictors that are independent of Y given X_d under the hypothesis.
;;;
          (if wts (send self :weights (select wts pan)))
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Next, Loop thru the predictors (dirs)
;;;
            (dotimes (j (- numtests 1))
                     (let* (
                            (k (+ j 1))
                            (u (copy-list ux))
                            )
                       (setf (nth j u) (select (nth j u) pan))
                       (send self :x u)
                       (setf (nth k count)
                             (append (nth k count)
                                     (list (sum
                                            (select
                                             (^ (send self :eigenvalues) power)
                                             (iseq k (- numpred 1)))
                                            ))))
                       )))
;;;
;;; Remove the initial nils in count
;;;
   (setf count (mapcar #'rest count))
;;;
;;; Return data set to original state.
;;;
   (send self :x ux0)
   (send self :y uy)
   (send self :weights wts)
;;;
;;; Calculate p-values
;;;
   (dotimes (i numtests)
            (setf (nth i pvals) (/ (length
                                    (which
                                     (< (sum (^ (select evs (iseq i (-
numpred 1))) power))
                                        (nth i count)))) B))
            )
;;;
;;; Add slots
;;;
   (send self :add-slot 'perm-pvals pvals)
   (send self :add-slot 'count count)
   (list pvals power)))

(defmeth inverse-regression-model-proto :perm-test ()
"Method args: ()
Retruns the permutation tests after calling :perm-test-pvals"
  (let* (
         (B (let ((ans (get-value-dialog "Number of permutations" :initial 50)))
              (if ans (first ans) (return-from :perm-test))))
         (nobs (send self :num-included))
         (var (^ (standard-deviation (send self :sir-y)) 2))
         (g 1) ; needed if phd
         (pvalspower (send self :perm-test-pvals B :num-lc *ir-num-lc*))
         (pvals (first pvalspower))
         (power (nth 1 pvalspower))
         (e (^ (reverse (first (send self :eigen-structure))) power))
         (p (length e))
         (method (send self :method))
         )
    (case method
      ((:PHD :PHDY)
       (setf g (* .5 (/ var))))
      (t))
    (send self :description1)
    (send self :description2)
    (format t "~%Permutation pvalues for")
    (format t "~%sums of eigenvalues times ~a" nobs)
    (format t "~%Number of permutations: ~a~%~%" B)
    (format t "Number of   ~13tTest~30tPermutation~%")
    (format t "Components  ~13tStatistic~32tp-value~%")
    (dotimes (i (length pvals))
            (let* ((a (* nobs (sum (select e (iseq (- p i))))))
                  (a (* a g))
                  )
              (format t "~4d~12t~11a~32t~5,3f~%" (+ i 1) (optfmt1 a) (nth i
pvals))))
    ))

(defmeth inverse-regression-model-proto :perm-pvalues ()
"Method args: ()
Returns the pvalues for the last permutation test; used by sir and save."
  (if (send self :needs-computing) (send self :compute))
  (slot-value 'perm-pvals))

(defmeth inverse-regression-model-proto :perm-test-stats ()
"Method args: ()
Returns the permutation test statistics as list of lists ordered on the
hypothesis starting with dimension=0. The statistics are the partial sums
of eigenvalues
eigenvalues squared and are not normalized"
  (if (send self :needs-computing) (send self :compute))
  (slot-value 'count))