/**********************************/
/* CIRCULAR STATISTICS FUNCTIONS  */
/* AUTHOR: Mathias Kölliker       */
/*         Zoological Insitute    */
/*         University of Basel    */
/*         Evolutionary Biology   */
/*         Vesalgasse 1           */
/*         4051 Basel             */
/*         - Switzerland          */
/* mathias.koelliker-at-unibas.ch */
/**********************************/

/***********************************************************************/
/* Function: 		circular_moments                                   */
/* Purpose:			calculates mean direction and mean vector-length   */
/*              	for a sample of angular values                     */
/* Parameters:		data_in: 	Input Data Set                         */
/*              	data_out:	Output Data Set containing mean        */
/*								direction and vector length            */
/*              	circvar:	Name of the circular variable          */
/*                  arcs:       Number of arcs to which the circle     */
/*                              was divided during data collection     */
/*                              (anz_group has to be set to "no" when  */
/*                              measurements are on a continuous scale)*/
/*              	rad_deg:	Statement on whether the circular      */
/*                  			values are on a radian (rad)or         */
/*                              degree (deg) scale                     */
/*                  remarks:    Remarks which will be added to the     */
/*                              output (e.g. sample-information        */
/* Author:      	Mathias Kölliker                                   */
/* Created:			14.02.2002                                         */
/* Changed:         14.02.2002; MK; newly written                      */
/*                  27.03.2002; MK; inclusion of circular standard dev */
/***********************************************************************/

%MACRO circular_moments (data_in, data_out, circvar, arcs, rad_deg, remarks);

 %IF &rad_deg = rad %THEN %DO;
 DATA temp01;
   SET &data_in;
       x = COS(&circvar);
	   y = SIN(&circvar);
 RUN;
 %END;
 %ELSE %IF &rad_deg = deg %THEN %DO;
  DATA temp01;
    SET &data_in;
      x = COS(&circvar / 360 * 6.28318530718);
	  y = SIN(&circvar / 360 * 6.28318530718);
  RUN;
  %END;
 
 PROC SUMMARY DATA = temp01;
   VAR x y;
     OUTPUT OUT = temp02 (RENAME = (_freq_ = n) DROP = _type_)
	   MEAN(x y) = mx my;
 RUN;

 %IF &rad_deg = rad %THEN %DO;
 DATA temp03 (DROP = arcs);
   SET temp02;
     arcs = SYMGET("arcs");
     IF arcs = "no" THEN c = 1;
	 ELSE DO;
       c = (6.28318530718 / (2 * arcs)) / SIN(6.28318530718 / (2 * arcs));
	 END;
     r = c * SQRT((mx * mx)+(my * my));
	 IF r > 1 THEN r = 1;
	 PHI = .;
	   IF mx > 0 THEN PHI = ATAN(my / mx);
	   IF mx < 0 THEN PHI = (180 * 6.28318530718 / 360) + ATAN(my / mx);
	   IF mx = 0 AND my > 0 THEN PHI = (90 * 6.28318530718 / 360);
	   IF mx = 0 AND my < 0 THEN PHI = (270 * 6.28318530718 / 360);
	 s2 = 2 * (1 - r);
	 s  = SQRT(2 * (1 - r));
	 s0 = SQRT(-2 * LOG(r));
 RUN;
 %END;
 %ELSE %IF &rad_deg = deg %THEN %DO;
 DATA temp03 (DROP = arcs);
   SET temp02;
     arcs = SYMGET("arcs");
     IF arcs = "no" THEN c = 1;
	 ELSE DO;
       c = (6.28318530718 / (2 * arcs)) / SIN(6.28318530718 / (2 * arcs));
	 END;
     r = c * SQRT((mx * mx)+(my * my));
	 IF r > 1 THEN r = 1;
	 PHI = .;
	   IF mx > 0 THEN PHI = ATAN(my / mx) * 360 / 6.28318530718;
	   IF mx < 0 THEN PHI = (180 + (ATAN(my / mx) * 360 / 6.28318530718));
	   IF mx = 0 AND my > 0 THEN PHI = 90;
	   IF mx = 0 AND my < 0 THEN PHI = 270;
	 s2 = (360 / 6.28318530718)**2 * 2 * (1 - r);
	 s  = (360 / 6.28318530718) * SQRT(2 * (1 - r));
	 s0 = (360 / 6.28318530718) * SQRT(-2 * LOG(r));
 RUN;
 %END;

 PROC TRANSPOSE DATA = temp03
                 OUT = temp04;
 RUN;

 DATA temp04 (DROP = _name_ col1);
   SET temp04;
     FORMAT Statistics $35.
            Value     7.4;
	 IF _name_ EQ "mx" OR _name_ EQ "my" THEN DELETE;
     IF _name_ = "n"   THEN Statistics = "Sample Size";
	 IF _name_ = "r"   THEN Statistics = "Mean Vector Length";
	 IF _name_ = "PHI" THEN Statistics = "Mean Angle (&rad_deg)";
	 IF _name_ = "s2"  THEN Statistics = "Angular Variance";
	 IF _name_ = "s"   THEN Statistics = "Angular Deviation";
	 IF _name_ = "s0"  THEN Statistics = "Circular Standard Deviation";
     Value = COL1;
	 IF Statistics = " " THEN DELETE;
 RUN;
 
 DATA &data_out;
   SET temp04;
 RUN;

 PROC PRINT DATA = &data_out NOOBS;
   VAR Statistics Value;
   TITLE HEIGHT = 2.5 "Circular Moments and Measures of Location";
   TITLE2 HEIGHT = 2.0 "&remarks";
   TITLE3 HEIGHT = 5.0 " ";
   TITLE4 HEIGHT = 5.0 " ";
   TITLE5 HEIGHT = 2.0 "SAS-Macro written by:";
   TITLE6 HEIGHT = 1.5 "Mathias Kölliker";
   TITLE7 HEIGHT = 1.5 "Zoological Institute, Evolutionary Biology";
   TITLE8 HEIGHT = 1.5 "University of Basel, Switzerland";
   TITLE9 HEIGHT = 1.5 "mathias.koelliker-at-unibas.ch";
   TITLE10 HEIGHT = 5.0 " ";
 RUN;
 QUIT;
 
 GOPTIONS RESET = ALL;

%MEND circular_moments;


/***********************************************************************/
/* Function: 		raileigh_test                                      */
/* Purpose:			calculates the moments of the circular distribution*/
/*                  and runs a Raileigh test to assess if there is     */
/*                  evidence for a significant direction               */
/* Parameters:		data_in: 	Input Data Set                         */
/*              	data_out:	Output Data Set containing mean        */
/*								direction and vector length            */
/*              	circvar:	Name of the circular variable          */
/*                  arcs:       Number of arcs to which the circle     */
/*                              was divided during data collection     */
/*                              (anz_group has to be set to "no" when  */
/*                              measurements are on a continuous scale)*/
/*              	rad_deg:	Statement on whether the circular      */
/*                  			values are on a radian (rad)or         */
/*                              degree (deg) scale                     */ 
/*                  remarks:    any kind of text. Will be written      */
/*      			            to the list-output                     */
/*                  iterations: number of iterations that are performed*/
/*                              in the randomization test              */
/*                              (default = 1000)                       */ 
/* Warning:         During macro-execution, the system option          */
/*                  "COMPRESS = " is set to "NO". If compression shall */
/*                  be enabled during the rest of the SAS-session, you */
/*                  have to re-run "OPTIONS COMPRESS = YES" somewhere  */
/*                  on the program-editor                              */
/* Author:      	Mathias Kölliker                                   */
/* Created:			14.02.2002                                         */
/* Changed:         14.02.2002; MK; newly written                      */
/*                  20.03.2002; MK; suppression of log-output during   */
/*                                  randomization                      */
/*                  20.03.2002; MK; nr. of iterations can be defined   */
/*                                  by the user (default = 0)          */
/***********************************************************************/

%MACRO rayleigh_test(data_in, data_out, circvar, arcs, rad_deg, remarks, iterations = 0);

%circular_moments(&data_in, &data_out, &circvar, &arcs, &rad_deg, &remarks);

 DATA temp05;
   SET temp04;
     WHERE Statistics = "Mean Vector Length";
	 CALL SYMPUT ("rm",value);
 RUN;
 
 DATA temp06;
   SET &data_in;
 RUN;

  /* Suppression of log-output */
 OPTIONS COMPRESS = NO NONOTES NOSOURCE NOSOURCE2 ERRORS = 0;

 /* Randomization */
 %LET i = 1;
 %DO %WHILE (&i LE &iterations);
   %IF &rad_deg = rad %THEN %DO;
     DATA temp06 (DROP = rand&i);
       SET temp06;
         rand&i = RANUNI(0) * 6.28318530718;
		 x&i = COS(rand&i);
	     y&i = SIN(rand&i);
     RUN;

	 PROC SUMMARY DATA = temp06;
       VAR x&i
           y&i;
       OUTPUT OUT = t1
       MEAN(x&i y&i) = mx my;
     RUN;
   %END;
   %ELSE %IF &rad_deg = deg %THEN %DO;
     DATA temp06 (DROP = rand&i);
	   SET temp06;
	     rand&i = RANUNI(0) * 360;
		 x&i = COS(rand&i / 360 * 6.28318530718);
	     y&i = SIN(rand&i / 360 * 6.28318530718);
	 RUN;

	 PROC SUMMARY DATA = temp06;
       VAR x&i
           y&i;
       OUTPUT OUT = t1
       MEAN(x&i y&i) = mx my;
     RUN;
   %END;

   %IF &i = 1 %THEN %DO;
   DATA temp07;
     SET t1;
   RUN;
   %END;
   %ELSE %DO;
   DATA temp07;
     SET temp07
	     t1;
   RUN;
   %END;

   %PUT Permutation &i;

 %LET i = %EVAL(&i + 1);
 %END;

  /* Reactivation of log-output */
 OPTIONS NOTES SOURCE SOURCE2 ERRORS = 20;
 
 DATA temp03;
   SET temp03;
     CALL SYMPUT("c",c);
 RUN;

 DATA temp07;
   SET temp07;
     r = SYMGET("c") * SQRT((mx * mx)+(my * my));
     p = 0;
     IF r >= SYMGET("rm") THEN p = 1;
 RUN;

 PROC SUMMARY DATA = temp07;
   VAR p;
   OUTPUT OUT = temp08 (DROP = _type_ _freq_)
     MEAN(p) = p;
 RUN;

 DATA temp08 (DROP = p);
   SET temp08;
     Statistics = "Raileigh-Test: P (Rho>=r)";
	 Value     = p;
 RUN;

 DATA &data_out;
   SET &data_out
       temp08;
	 FORMAT Iterations  9.;
	 IF Statistics = "Raileigh-Test: P (Rho>=r)" THEN
	    Iterations = &iterations;
 RUN;

 PROC PRINT DATA = &data_out NOOBS;
   VAR Statistics Value Iterations;
   TITLE HEIGHT = 2.5 "Circular Moments and Measures of Location";
   TITLE2 HEIGHT = 2.5 "Raileigh-Test for Uniform Distribution"; 
   TITLE3 HEIGHT = 2.0 "&remarks";
   TITLE4 HEIGHT = 5.0 " ";
   TITLE5 HEIGHT = 5.0 " ";
   TITLE6 HEIGHT = 2.0 "SAS-Macro written by:";
   TITLE7 HEIGHT = 1.5 "Mathias Kölliker";
   TITLE8 HEIGHT = 1.5 "Zoological Institute, Evolutionary Biology";
   TITLE9 HEIGHT = 1.5 "University of Basel, Switzerland";
   TITLE10 HEIGHT = 1.5 "mathias.koelliker-at-unibas.ch";
   TITLE11 HEIGHT = 5.0 " ";
 RUN;
 QUIT;

 GOPTIONS RESET = ALL;

%MEND rayleigh_test;


/***********************************************************************/
/* Function: 		circular_corr                                      */
/* Purpose:			calculates various circular correlation            */
/*                  coefficients among two circular variables          */
/* Parameters:		data_in: 	Input Data Set                         */
/*              	data_out:	Output Data Set containing the         */
/*                              calculated correlation coefficient     */
/*              	circvar1:	Name of the first circular variable    */
/*                  curcvar2;   Name of the second circular variable   */
/*              	rad_deg:	Statement on whether the circular      */
/*                  			values are on a radian (rad)or         */
/*                              degree (deg) scale                     */
/*                  iterations: number of iterations that are performed*/
/*                              in the randomization test              */
/*                              (default = 0)                          */ 
/* Warning:         During macro-execution, the system option          */
/*                  "COMPRESS = " is set to "NO". If compression shall */
/*                  be enabled during the rest of the SAS-session, you */
/*                  have to re-run "OPTIONS COMPRESS = YES" somewhere  */
/*                  on the program-editor                              */
/* Author:      	Mathias Kölliker                                   */
/* Created:			07.03.2002                                         */
/* Changed:         07.03.2002; MK; newly written                      */
/*                  18.03.2002; MK; inclusion of non-parametric        */
/*                                  correlation                        */
/*                  20.03.2002; MK; suppression of log-output during   */
/*                                  randomization                      */
/*                  20.03.2002; MK; nr. of iterations can be defined   */
/*                                  by the user (default = 0)          */
/*                  24.03.2002; MK; rescaling of coefficient to range  */
/*   				                from 0 to 1                        */
/*                  27.03.2002; MK; inclusion of correlation coef.     */
/*                                  following Fisher and Lee           */
/*                  08.04.2002; MK; modular structure, inclusion of    */
/*                                  corr coef following Jammalamadaka  */
/*                                  and Sarma 1988                     */ 
/*                  25.04.2002; MK; controlable output, default is     */
/*                                  WITH list output                   */
/***********************************************************************/

%MACRO circular_corr(data_in, data_out, circvar1, circvar2, rad_deg, iterations = 0,
                     output = yes);

 /* Suppression of log-output */
 OPTIONS COMPRESS = NO NONOTES NOSOURCE NOSOURCE2 ERRORS = 0;

 %IF &rad_deg = rad %THEN %DO;
 DATA temp01;
   SET &data_in;
     circvar1 = &circvar1;
	 circvar2 = &circvar2;
	 cv1 = COS(circvar1);
	 sv1 = SIN(circvar1);
	 cv2 = COS(circvar2);
	 sv2 = sin(circvar2);
   IF &circvar1 NE . AND &circvar2 NE . THEN OUTPUT;
 RUN;
 %END;

 %ELSE %IF &rad_deg = deg %THEN %DO;
 DATA temp01;
   SET &data_in;
     circvar1 = &circvar1 / 360 * 6.28318530718;
	 circvar2 = &circvar2 / 360 * 6.28318530718;
	 cv1 = COS(circvar1);
	 sv1 = SIN(circvar1);
	 cv2 = COS(circvar2);
	 sv2 = sin(circvar2);
   IF &circvar1 NE . AND &circvar2 NE . THEN OUTPUT;
 RUN;
 %END;

 PROC SUMMARY DATA = temp01;
   OUTPUT OUT = temp_n;
 RUN;

 DATA temp_n;
   SET temp_n;
     CALL SYMPUT("n",_freq_);
 RUN;

 %jupp_mardia(temp01, temp02a);
 DATA temp02a;
   SET temp02a;
     CALL SYMPUT("ra", r);
 RUN;
     
 %fisher_lee_parm(temp01, temp02b);
 DATA temp02b;
   SET temp02b;
 	 CALL SYMPUT("rb", r);
 RUN;

 %fisher_lee_nparm(temp01, temp02c);
 DATA temp02c;
   SET temp02c;
     CALL SYMPUT("rc", r);
 RUN;    

 %jammalamadaka(temp01, temp02d);
 DATA temp02d;
   SET temp02d;
     CALL SYMPUT("rd", r);
 RUN;

 /*Randomization Test*/
 DATA temp01a_r;
   SET temp01 (KEEP = circvar1);
 DATA temp01b_r;
   SET temp01 (KEEP = circvar2);
 RUN;

 %LET i = 1;
 %DO %WHILE (&i LE &iterations);
 DATA temp02a_r;
   SET temp01a_r;
 RUN;

 DATA temp02b_r;
   SET temp01b_r;
     r_&i = ranuni(0);
 PROC SORT DATA = temp02b_r;
   BY r_&i;
 RUN;

 DATA temp03_r (DROP = r_&i);
   MERGE temp02a_r
         temp02b_r;
 RUN;

 DATA temp03_r;
   SET temp03_r;
	 cv1 = COS(circvar1);
	 sv1 = SIN(circvar1);
	 cv2 = COS(circvar2);
	 sv2 = sin(circvar2);
 RUN;

 %jupp_mardia(temp03_r, temp04a_r);
 %fisher_lee_parm(temp03_r, temp04b_r);
 %fisher_lee_nparm(temp03_r, temp04c_r);
 %jammalamadaka(temp03_r, temp04d_r);

 %IF &i = 1 %THEN %DO;
   DATA temp05a_r;
     SET temp04a_r;
   DATA temp05b_r;
     SET temp04b_r;
   DATA temp05c_r;
     SET temp04c_r;
   DATA temp05d_r;
     SET temp04d_r;
   RUN;
 %END;
 %ELSE %DO;
   DATA temp05a_r;
     SET temp04a_r
         temp05a_r;
   DATA temp05b_r;
     SET temp04b_r
         temp05b_r;
   DATA temp05c_r;
     SET temp04c_r
	     temp05c_r;
   DATA temp05d_r;
     SET temp04d_r
	     temp05d_r;
   RUN;
 %END;

 %PUT Permutation &i;

 %LET i = %EVAL(&i + 1);
 %END;

 /* P-value assessments / Output generation*/
 DATA temp05a_r;
   SET temp05a_r;
     IF r EQ . THEN p = .;
     ELSE IF ABS(r) GE ABS(SYMGET("ra")) THEN p = 1;
	 ELSE IF ABS(r) LT ABS(SYMGET("ra")) THEN p = 0;
 PROC SUMMARY DATA = temp05a_r;
   VAR p;
     OUTPUT OUT = temp06a_r
       MEAN(p) = p;
 DATA temp06a_r;
   SET temp06a_r;
     CALL SYMPUT("pa",p);
 DATA temp02a (KEEP = Coefficient Sample_Size Value P_Value Permutations Info Ref);
   SET temp02a;
     FORMAT Coefficient   $5.
			Sample_Size   5.0
	        Value         10.5
			P_Value       4.3
			Permutations  6.0
            Info          $20.
            Ref           $20.;
	   Coefficient = "r(2)";
	   Sample_Size = SYMGET("n");
	   Value       = r;
	   P_Value     = SYMGET("pa");
	   Permutations= &iterations;
	   Info        = "parametric";
	   Ref         = "Jupp-Mardia (1980)";
	IF Value = . THEN DELETE;
 RUN;

 DATA temp05b_r;
   SET temp05b_r;
     IF r EQ . THEN p = .;
     ELSE IF ABS(r) GE ABS(SYMGET("rb")) THEN p = 1;
	 ELSE IF ABS(r) LT ABS(SYMGET("rb")) THEN p = 0;
PROC SUMMARY DATA = temp05b_r;
   VAR p;
     OUTPUT OUT = temp06b_r
       MEAN(p) = p;
 DATA temp06b_r;
   SET temp06b_r;
     CALL SYMPUT("pb",p);
 DATA temp02b (KEEP = Coefficient Sample_Size Value P_Value Permutations Info Ref);
   SET temp02b;
     FORMAT Coefficient   $5.
			Sample_Size   5.0
	        Value         10.5
			P_Value       4.3
			Permutations  6.0
            Info          $20.
            Ref           $20.;
	   Coefficient = "r(t)";
	   Sample_Size = SYMGET("n");
	   Value       = r;
	   P_Value     = SYMGET("pb");
	   Permutations= &iterations;
	   Info        = "parametric";
	   Ref         = "Fisher-Lee (1983)";
	IF Value = . THEN DELETE;
 RUN;

 DATA temp05c_r;
   SET temp05c_r;
     IF r EQ . THEN p = .;
     ELSE IF ABS(r) GE ABS(SYMGET("rc")) THEN p = 1;
	 ELSE IF ABS(r) LT ABS(SYMGET("rc")) THEN p = 0;
PROC SUMMARY DATA = temp05c_r;
   VAR p;
     OUTPUT OUT = temp06c_r
       MEAN(p) = p;
 DATA temp06c_r;
   SET temp06c_r;
     CALL SYMPUT("pc",p);
 DATA temp02c (KEEP = Coefficient Sample_Size Value P_Value Permutations Info Ref);
   SET temp02c;
     FORMAT Coefficient   $5.
			Sample_Size   5.0
	        Value         10.5
			P_Value       4.3
			Permutations  6.0
            Info          $20.
            Ref           $20.;
	   Coefficient = "r(n)";
	   Sample_Size = SYMGET("n");
	   Value       = r;
	   P_Value     = SYMGET("pc");
	   Permutations= &iterations;
	   Info        = "non-parametric";
	   Ref         = "Fisher-Lee (1982)";
	IF Value = . THEN DELETE;
 RUN;

 DATA temp05d_r;
   SET temp05d_r;
     IF r EQ . THEN p = .;
     ELSE IF ABS(r) GE ABS(SYMGET("rd")) THEN p = 1;
	 ELSE IF ABS(r) LT ABS(SYMGET("rd")) THEN p = 0;
 PROC SUMMARY DATA = temp05d_r;
   VAR p;
     OUTPUT OUT = temp06d_r
       MEAN(p) = p;
 DATA temp06d_r;
   SET temp06d_r;
     CALL SYMPUT("pd",p);
 DATA temp02d (KEEP = Coefficient Sample_Size Value P_Value Permutations Info Ref);
   SET temp02d;
     FORMAT Coefficient   $5.
			Sample_Size   5.0
	        Value         10.5
			P_Value       4.3
			Permutations  6.0
            Info          $20.
            Ref           $20.;
	   Coefficient = "r(cn)";
	   Sample_Size = SYMGET("n");
	   Value       = r;
	   P_Value     = SYMGET("pd");
	   Permutations= &iterations;
	   Info        = "parametric";
	   Ref         = "Rao-SenGupta (2001)";
	IF Value = . THEN DELETE;
 RUN;

 DATA &data_out;
   SET temp02a
       temp02b
       temp02c
       temp02d;
   IF SYMGET("iterations") LE 50 THEN P_Value = " ";
 RUN;

 PROC TRANSPOSE DATA = &data_out
                 OUT = &data_out;
   VAR Coefficient Sample_Size Value P_Value Permutations Info Ref;
 RUN;

 DATA &data_out (KEEP = Statistics r_2 r_t r_n r_cn);
   SET &data_out;
     FORMAT Statistics             $11.
            r_2 r_t r_n r_cn       $20.;
     Statistics = _name_;
     r_2   = LEFT(TRIM(COL1));
	 r_t   = LEFT(TRIM(COL2));
	 r_n  = LEFT(TRIM(COL3));
	 r_cn  = LEFT(TRIM(COL4));
   IF _name_ = "Coefficient" THEN DELETE;
 RUN;

   /* Reactivation of log-output */
 OPTIONS NOTES SOURCE SOURCE2 ERRORS = 20;

 %IF &output = yes %THEN %DO;
 PROC PRINT DATA = &data_out NOOBS;
   VAR Statistics r_2 r_t r_n r_cn;
   TITLE HEIGHT  = 2.5 "Circular Correlation Coefficients";
   TITLE2 HEIGHT = 5.0 " ";
   TITLE3 HEIGHT = 5.0 " ";
   TITLE4 HEIGHT = 2.0 "SAS-Macro written by:";
   TITLE5 HEIGHT = 2.0 "Mathias Kölliker";
   TITLE6 HEIGHT = 1.5 "Zoological Institute, Evolutionary Biology";
   TITLE7 HEIGHT = 1.5 "University of Basel, Switzerland";
   TITLE8 HEIGHT = 1.5 "mathias.koelliker-at-unibas.ch";
   TITLE9 HEIGHT = 5.0 " ";
 RUN;
 QUIT;
 
 GOPTIONS RESET = ALL;
 %END;

%MEND circular_corr;

/*********************************************************************/
/* Function:		jupp_mardia                                      */
/* Purpose:			calculates a parametric circular correlation     */
/*                  coefficient following Jupp and Mardia (1980;     */
/*                  Biometrika 67, 163-173)                          */
/* Warning:         plug-in for function circular_corr               */
/*********************************************************************/

%MACRO jupp_mardia(data_in, data_out);

 PROC CORR DATA = &data_in NOPRINT NOSIMPLE OUTP = tempp_cc;
   VAR cv1 cv2;
 DATA tempp_cc (KEEP = rcc);
   SET tempp_cc;
     WHERE UPCASE(_type_) = "CORR" AND cv1 = 1;
       RENAME cv2 = rcc;
 
 PROC CORR DATA = &data_in NOPRINT NOSIMPLE OUTP = tempp_cs;
   VAR cv1 sv2;
 DATA tempp_cs (KEEP = rcs);
   SET tempp_cs;
     WHERE UPCASE(_type_) = "CORR" AND cv1 = 1;
       RENAME sv2 = rcs;

 PROC CORR DATA = &data_in NOPRINT NOSIMPLE OUTP = tempp_sc;
   VAR sv1 cv2;
 DATA tempp_sc (KEEP = rsc);
   SET tempp_sc;
     WHERE UPCASE(_type_) = "CORR" AND sv1 = 1;
	   RENAME cv2 = rsc;

 PROC CORR DATA = &data_in NOPRINT NOSIMPLE OUTP = tempp_ss;
   VAR sv1 sv2;
 DATA tempp_ss (KEEP = rss);
   SET tempp_ss;
     WHERE UPCASE(_type_) = "CORR" AND sv1 = 1;
	   RENAME sv2 = rss;

 PROC CORR DATA = &data_in NOPRINT NOSIMPLE OUTP = tempp_var1;
   VAR sv1 cv1;
 DATA tempp_var1 (KEEP = r1);
   SET tempp_var1;
     WHERE UPCASE(_type_) = "CORR" AND sv1 = 1;
	   RENAME cv1 = r1;

 PROC CORR DATA = &data_in NOPRINT NOSIMPLE OUTP = tempp_var2;
   VAR sv2 cv2;
 DATA tempp_var2 (KEEP = r2);
   SET tempp_var2;
     WHERE UPCASE(_type_) = "CORR" AND sv2 = 1;
	   RENAME cv2 = r2;

 DATA &data_out;
   MERGE tempp_cc tempp_cs tempp_sc tempp_ss tempp_var1 tempp_var2;
     r = ((rcc**2 + rcs**2 + rsc**2 + rss**2) + (2*r1*r2*(rcc*rss + rcs*rsc))
	     - (2*r2*(rcc*rcs + rsc*rss)) - (2*r1*(rcc*rsc + rcs*rss))) /
		 ((1 - r1**2)*(1 - r2**2));
 RUN;

%MEND jupp_mardia;



/*********************************************************************/
/* Function:		fisher_lee_nparm                                 */
/* Purpose:			calculates a non-parametric circular correlation */
/*                  coefficient following Fisher and Lee (1982;      */
/*                  Biometrika 69, 315-321) and Fisher (1993;        */
/*                  Statistical analysis of circular data.           */
/*                  Cambridge University Press)                      */
/* Warning:         plug-in for function circular_corr               */
/*********************************************************************/

%MACRO fisher_lee_nparm(data_in, data_out);

 PROC RANK FRACTION DATA = &data_in
                     OUT = temp02;
	VAR circvar1 circvar2;
	RANKS cr1 cr2;
 DATA temp02;
   SET temp02 (KEEP = circvar1 circvar2 cr1 cr2);
     C    = 6.28318530718;
     cr1 = C*cr1;
	 cr2 = C*cr2;
     tcc = COS(cr1) * COS(cr2);
	 tss = SIN(cr1) * SIN(cr2);
	 tsc = SIN(cr1) * COS(cr2);
	 tcs = COS(cr1) * SIN(cr2);
PROC SUMMARY DATA = temp02;
   VAR tcc tss tcs tsc;
     OUTPUT OUT = &data_out
	   MEAN() = ;
 DATA &data_out;
   SET &data_out;
     n  = SYMGET("n");
     R1 = (tcc + tss)**2 + (tsc - tcs)**2;
	 R2 = (tcc - tss)**2 + (tsc + tcs)**2;
	 r  = r1 - r2;
 RUN;

%MEND fisher_lee_nparm;



/*********************************************************************/
/* Function:		fisher_lee_parm                                  */
/* Purpose:			calculates a parametric circular correlation     */
/*                  coefficient following Fisher and Lee (1983;      */
/*                  Biometrika 70, 327-332) and Fisher (1993;        */
/*                  Statistical analysis of circular data.           */
/*                  Cambridge University Press)                      */
/* Warning:         plug-in for function circular_corr               */
/*********************************************************************/

%MACRO fisher_lee_parm(data_in, data_out);

/* Formula from Fisher (1993) page 151 */
 DATA temp02;
   SET &data_in;
     A = cv1 * cv2;
	 B = sv1 * sv2;
	 C = cv1 * sv2;
	 D = sv1 * cv2;
	 E = COS(2 * circvar1);
	 F = SIN(2 * circvar1);
	 G = COS(2 * circvar2);
	 H = SIN(2 * circvar2);
 PROC SUMMARY DATA = temp02;
   VAR A B C D E F G H;
     OUTPUT OUT = &data_out (DROP = _type_ RENAME = (_freq_ = n))
	   SUM() = ;
 DATA &data_out;
   SET &data_out;
     r = 4 * (A * B - C * D) / SQRT((n**2 - E**2 - F**2) * (n**2 - G**2 - H**2));
 RUN;

%MEND fisher_lee_parm;

/*********************************************************************/
/* Function:		jammalamadaka                                    */
/* Purpose:			calculates a parametric circular correlation     */
/*                  coefficient following Jammalamadaka and          */
/*                  SenGupta (2001; Topics in circular statistics.   */
/*                  Series on multivariate Analysis.                 */
/*                  World Scientific Publishing Co.), page 178       */
/* Warning:         plug-in for function circular_corr               */
/*********************************************************************/

%MACRO jammalamadaka(data_in, data_out);

 DATA temp02;
   SET &data_in (KEEP = circvar1 circvar2 cv1 sv1 cv2 sv2);
 PROC SUMMARY DATA = temp02;
   VAR cv1 sv1 cv2 sv2;
     OUTPUT OUT = temp03 (RENAME = (_freq_ = n) DROP = _type_)
	   MEAN(cv1 sv1 cv2 sv2) = mx1 my1 mx2 my2;
  DATA temp03;
   SET temp03;
	 PHI_1 = .;
	   IF mx1 > 0 THEN PHI_1 = ATAN(my1 / mx1);
	   IF mx1 < 0 THEN PHI_1 = (180 * 6.28318530718 / 360) + ATAN(my1 / mx1);
	   IF mx1 = 0 AND my1 > 0 THEN PHI_1 = (90 * 6.28318530718 / 360);
	   IF mx1 = 0 AND my1 < 0 THEN PHI_1 = (270 * 6.28318530718 / 360);
	 PHI_2 = .;
	   IF mx2 > 0 THEN PHI_2 = ATAN(my2 / mx2);
	   IF mx2 < 0 THEN PHI_2 = (180 * 6.28318530718 / 360) + ATAN(my2 / mx2);
	   IF mx2 = 0 AND my2 > 0 THEN PHI_2 = (90 * 6.28318530718 / 360);
	   IF mx2 = 0 AND my2 < 0 THEN PHI_2 = (270 * 6.28318530718 / 360);
	 CALL SYMPUT("phi_1", PHI_1);
	 CALL SYMPUT("phi_2", PHI_2);
 DATA temp02;
   SET temp02;
	 A    = SIN(circvar1 - SYMGET("phi_1")) * SIN(circvar2 - SYMGET("phi_2"));
	 B    = SIN(circvar1 - SYMGET("phi_1"))**2 * SIN(circvar2 - SYMGET("phi_2"))**2;
 PROC SUMMARY DATA = temp02;
   VAR A B;
     OUTPUT OUT = &data_out (DROP = _type_ RENAME = (_freq_ = n))
	   SUM(A B)  = A B;
 DATA &data_out;
   SET &data_out;
     r = A / (SQRT(n-1) * SQRT(B));
 RUN;

%MEND jammalamadaka;






/***********************************************************************/
/* Function: 		two_sample_test                                    */
/* Purpose:			tests the siginificance of the difference between  */
/*                  the mean angles of two independent samples         */
/*                  following the parametric method by                 */
/*                  Watson and Williams                                */
/* Parameters:		data_in: 	Input Data Set                         */
/*              	data_out:	Output Data Set containing the         */
/*                              test statistics and p-value            */
/*              	circvar:	Name of the circular variable          */
/*                  groupvar;   Name of the grouping variable          */
/*              	rad_deg:	Statement on whether the circular      */
/*                  			values are on a radian (rad)or         */
/*                              degree (deg) scale                     */
/*                  iterations: number of iterations that are performed*/
/*                              in the randomization test              */
/*                              (default = 0)                          */ 
/* Warnings:        -During macro-execution, the system option         */
/*                   "COMPRESS = " is set to "NO". If compression shall*/
/*                   be enabled during the rest of the SAS-session, you*/
/*                   have to re-run "OPTIONS COMPRESS = YES" somewhere */
/*                   on the program-editor                             */
/*                  -the grouping variable must have the values 1 and 2*/
/* Author:      	Mathias Kölliker                                   */
/* Created:			24.03.2002                                         */
/* Changed:         24.03.2002; MK; newly written                      */
/***********************************************************************/

%MACRO two_sample_test(data_in, data_out, circvar, groupvar, rad_deg, iterations = 0);

%IF &rad_deg = rad %THEN %DO;
 DATA temp01;
   SET &data_in (KEEP = &circvar &groupvar);
     WHERE &groupvar NE .;
       C = COS(&circvar);
	   S = SIN(&circvar);
 RUN;
 %END;
 %ELSE %IF &rad_deg = deg %THEN %DO;
 DATA temp01;
   SET &data_in (KEEP = &circvar &groupvar);
     WHERE &groupvar NE .;
       C = COS(&circvar / 360 * 6.28318530718);
	   S = SIN(&circvar / 360 * 6.28318530718);
 RUN;
 %END;

 PROC SORT DATA = temp01;
   BY &groupvar;
 RUN;

 DATA temp02;
   SET temp01 (OBS = 1);
     CALL SYMPUT("gr1", &groupvar);
 RUN;

 DATA temp01;
   SET temp01;
     IF &groupvar = SYMGET("gr1") THEN group = 1;
	 ELSE group = 2;
 RUN;

 DATA temp03a temp03b;
   SET temp01;
     IF group = 1 THEN OUTPUT temp03a;
	 IF group = 2 THEN OUTPUT temp03b;
 RUN;

 PROC SUMMARY DATA = temp03a;
   VAR C S;
     OUTPUT OUT = temp04a (DROP = _type_ RENAME = (_freq_ = n1))
	   SUM(C S) = C1 S1;
 RUN;

 PROC SUMMARY DATA = temp03b;
   VAR C S;
     OUTPUT OUT = temp04b (DROP = _type_ RENAME = (_freq_ = n2))
	   SUM(C S) = C2 S2;
 RUN;

 DATA temp05;
   MERGE temp04a
         temp04b;
     R1 = SQRT(C1**2 + S1**2);
	 R2 = SQRT(C2**2 + S2**2);
	 n  = n1 + n2;
	 C  = C1 + C2;
	 S  = S1 + S2;
	 R  = SQRT(C**2 + S**2);
	 d  = R1 + R2 - R;
	 CALL SYMPUT("d", d);
 RUN;

 /*Randomisation Test*/
 DATA temp06a (KEEP = group) temp06b (KEEP = C S);
   SET temp01;
 RUN;

  /* Suppression of log-output */
 OPTIONS COMPRESS = NO NONOTES NOSOURCE NOSOURCE2 ERRORS = 0;

 %LET i = 1;
 %DO %WHILE (&i LE &iterations);
   DATA temp07a;
     SET temp06a;
       r_&i = ranuni(0);
   PROC SORT DATA = temp07a;
     BY r_&i;
   RUN;

   DATA temp07b;
     SET temp06b;
      r_&i = ranuni(0);
   PROC SORT DATA = temp07b;
     BY r_&i;
   RUN;

   DATA temp08;
     MERGE temp07a
           temp07b;
   RUN;

   DATA temp09a temp09b;
     SET temp08;
       IF group = 1 THEN OUTPUT temp09a;
	   IF group = 2 THEN OUTPUT temp09b;
   RUN;

   PROC SUMMARY DATA = temp09a;
     VAR C S;
       OUTPUT OUT = temp10a (DROP = _type_ RENAME = (_freq_ = n1))
	     SUM(C S) = C1 S1;
   RUN;

   PROC SUMMARY DATA = temp09b;
     VAR C S;
       OUTPUT OUT = temp10b (DROP = _type_ RENAME = (_freq_ = n2))
	     SUM(C S) = C2 S2;
   RUN;

   DATA temp11;
     MERGE temp10a
           temp10b;
       R1 = SQRT(C1**2 + S1**2);
	   R2 = SQRT(C2**2 + S2**2);
	   n  = n1 + n2;
	   C  = C1 + C2;
	   S  = S1 + S2;
	   R  = SQRT(C**2 + S**2);
	   dr  = R1 + R2 - R;
   RUN;

   %IF &i = 1 %THEN %DO;
     DATA temp12;
       SET temp11 (KEEP = dr);
     RUN;
   %END;
   %ELSE %DO;
     DATA temp12;
       SET temp12
	       temp11;
     RUN;
   %END;

   %PUT Permutation &i;

   %LET i = %EVAL(&i + 1);
 %END;

  /* Reactivation of log-output */
 OPTIONS NOTES SOURCE SOURCE2 ERRORS = 20;
   
 DATA temp12;
   SET temp12;
     p = 0;
	 IF dr GE SYMGET("d") THEN p = 1;
 RUN;

 PROC SUMMARY DATA = temp12;
   VAR p;
     OUTPUT OUT = temp13
	   MEAN(p) = ;
 RUN;

 DATA temp13;
   SET temp13;
     CALL SYMPUT("p", p);
 RUN;

 DATA &data_out (DROP = n1 n2 d);
   SET temp05 (KEEP = n1 n2 d);
     FORMAT n_1 n_2 6.
	        Test_stats   $5.
			Value             10.5
			P_value           4.3;
     n_1        = n1;
	 n_2        = n2;
	 Test_stats = d;
	 P_value    = SYMGET("p");
 RUN;

 PROC PRINT DATA = &data_out NOOBS;
   VAR n_1 n_2 Test_stats P_value;
   TITLE HEIGHT  = 2.5 "Watson-Williams two-sample test";
   TITLE2 HEIGHT = 5.0 " ";
   TITLE3 HEIGHT = 5.0 " ";
   TITLE4 HEIGHT = 2.0 "SAS-Macro written by:";
   TITLE5 HEIGHT = 2.0 "Mathias Kölliker";
   TITLE6 HEIGHT = 1.5 "Zoological Institute, Evolutionary Biology";
   TITLE7 HEIGHT = 1.5 "University of Basel, Switzerland";
   TITLE8 HEIGHT = 1.5 "mathias.koelliker-at-unibas.ch";
   TITLE9 HEIGHT = 5.0 " ";
 RUN;
 QUIT;
 
 GOPTIONS RESET = ALL;

%MEND two_sample_test;