options pagesize=150 linesize=68;
filename grafout 'g10.ps';
goptions device=pscolor gsfname=grafout 
   rotate=landscape gsfmode=append;
/******************************************************/
/* Stuart (1953 Biometrika 40, 105-110), quoted by    */
/* Stuart (1955 Biometrika 42, 412-416), and Agresti  */
/* (1983 Biometrics 39(2) 505-510), present an ordinal*/
/* categorization of right and left eye vision of     */
/* Royal Ordnance female factory workers 1943-46.     */
/******************************************************/
data eyes; infile 'eyes.dat'; input r l count; run;
/**************************************************************/
/* Mark A: Tests of Symmetry                                  */
/**************************************************************/
* Bowker's test is labeled S.;
* Output from next block of code below will be used later;
title 'Eyesight example: test symmetry';
proc freq data=eyes; weight count; 
   table r*l/norow nocol nopct agree ; run;
/************************************************************/
/* Election data from 1997, 2001, 2005 for London boroughs  */
/* C=conservative, L=labor, LD=liberal democrat, RU=respect */
/* union.  From http://www.election.demon.co.uk/1997LB.html */
/* RU was a small recent party most ideologically related to*/
/* labor.
/************************************************************/
data elections; infile 'elections.dat';
   input w97 $ w01$ w05$; w05c=w05; if w05="RU" then w05c="L";
   run;
title 'Electon example: test symmetry';
/**************************************************************/
/* With zero denominators, we have a problem.                 */
/**************************************************************/
* Sas doesn't handle this correctly.;
proc freq data=elections; 
   table w97*w05/norow nocol nopct agree; run;
proc freq data=elections; 
   table w97*w05c/norow nocol nopct agree; run;
/**************************************************************/
/* Fit the symmetry model via Poisson reg.  Make covariates to*/
/* represent diagonal entries, and covariates for each sym-   */
/* metry pair.  If J=number of rows and columns, first J are  */
/* diagonal entries and next J*(J-1)/2 symmetry terms.  Neg-  */
/* ative multiplier makes 0 the baseline category.            */
/**************************************************************/
title 'Eye Symmetry via Poisson Regression';
data neweyes; set eyes;
   md=0;
   do i=1 to 4;
      md=md-i*(l=i)*(r=i);
   end;
   k=4;
   do i=1 to 3;
      do j=i+1 to 4;
         k=k+1;
         md=md-((r=i)*(l=j)+(r=j)*(l=i))*k;
      end;
   end;
   drop i j k; run;
/**************************************************************/
/*Fit the Poisson model.  Specify no main effects, so that the*/
/*probabilities really are symmetric.  Drop the intercept to  */
/*make the output cleaner.  Examining fitted values verifies  */
/*that off-diagonal fits are the averages of the two counts.  */
/**************************************************************/
*  Use the output statement, with out=        ;
* and predicted= to create data set   ;
* containing variable  with predictions .     ;
title 'Eyesight example: fit symmetry via genmod';
proc genmod data=neweyes; class md;
   model count=md /noint dist=poisson;
   output out=fitted predicted=outpred; run;
proc print data=fitted; run;
/**************************************************************/
/*Verify that the resulting test statistic is similar to      */
/*Bowker's test.  Both have the same DF.                      */
/**************************************************************/
title 'Eyesight example: fit sat. model to test symmetry';
proc genmod data=neweyes; class r l;
   model count=r*l/noint dist=poisson; 
   contrast 'r0' r*l 0  1  0  0
                    -1  0  0  0 
                     0  0  0  0 
                     0  0  0  0,
               r*l   0  0  1  0
                     0  0  0  0 
                    -1  0  0  0 
                     0  0  0  0,
               r*l   0  0  0  1
                     0  0  0  0 
                     0  0  0  0 
                    -1  0  0  0,
               r*l   0  0  0  0
                     0  0  1  0 
                     0 -1  0  0 
                     0  0  0  0,
               r*l   0  0  0  0
                     0  0  0  1 
                     0  0  0  0 
                     0 -1  0  0,
               r*l   0  0  0  0
                     0  0  0  0 
                     0  0  0  1 
                     0  0 -1  0; run;
* The contrast statement above forces calculation of LR stat-  ;
* istics for the null hypothesis of symmetry, defined by 6 con-;
* trasts on parameters.  First equates elements (2,1) and (1,2);
* by putting +1 and -1 and forcing linear combination to be 0. ;
/**************************************************************/
/* Mark B: Quasisymmetry                                      */
/**************************************************************/
* Test quasi-symmetry model using LR.                          ;
title 'Eyesight example: fit sat. model to test quasisym.';
proc genmod data=neweyes; class r l;
   model count=r*l/noint dist=poisson; 
   contrast 'r0' r*l 0  1 -1  0
                    -1  0  1  0
                     1 -1  0  0
                     0  0  0  0,
                 r*l 0  1  0 -1
                    -1  0  0  1
                     0  0  0  0
                     1 -1  0  0,
                 r*l 0  0  1 -1
                     0  0  0  0
                    -1  0  0  1
                     1  0 -1  0; run;
* You might hope that things would be easier than the above,; 
* since quasisymmetry really should involve just the inter- ; 
* actions in a two-way model with main effects included, but;
* I've done the calculations coding in by hand interactions ;
* in a model that SAS is treating as an otherwise unstruc-  ;
* tured multinomail.  As far as I can tell, there's no eas- ;
* ier way around this, because SAS won't let you add con-   ;
* trasts that don't sum to 0 within a given level of a      ;
* lower-order term.  You can also do this using PROC CATMOD,;
* but without any serious simplification.  Results differ   ;
* lightly from those of GENMOD, since CATMOD estimation is  ;
* via weighted least squares, and test is Wald. I show both.; 
title 'Quasisymmetry via CATMOD'; title2;
proc catmod data=neweyes; weight count;
   response 
      0  1 -1  0 -1  0  1  0  1 -1  0  0 0  0  0  0,
      0  1  0 -1 -1  0  0  1  0  0  0  0 1 -1  0  0,
      0  0  1 -1  0  0  0  0 -1  0  0  1 1  0 -1  0 log;
   model r*l=/noint; run;
*The above calculation contrasts the log of the expectations,;
*and show that quasisymmetry can't be ruled out.             ;
*Mark C: Marginal symmetry.  Lack of log in response makes   ;
*summation on the original scale.  R doesn't offer this.     ;
title 'Marginal symmetry via CATMOD';
proc catmod data=neweyes; weight count;
   response 0  1  1  1  
           -1  0  0  0  
           -1  0  0  0  
           -1  0  0  0,
            0 -1  0  0  
            1  0  1  1  
            0 -1  0  0  
            0 -1  0  0,
            0  0 -1  0  
            0  0 -1  0  
            1  1  0  1  
            0  0 -1  0;
   model r*l=/noint; run;
/**************************************************************/
/*Mark D: Test the quasi-independence model.  Create a vari-  */
/*able d that gives the row and column number for diagonal    */
/*elements.  Use this variable as a factor.  Make the baseline*/
/*be off-diagonal by setting d=5.  Test this alt. *hypothesis */
/*HA vs. H0 of independence using Type 3 analysis.  Test H0 of*/
/*quasi-independence vs. HA of the saturated model using  the */
/*likelihood for HA above.                                    */
/**************************************************************/
title 'Quasiindependence via GENMOD';
data neweyes; set neweyes; 
   d=(r=1)*(l=1)*1+(r=2)*(l=2)*2+(r=3)*(l=3)*3+(r=4)*(l=4)*4;
   if(d=0) then d=5; run;
proc genmod data=neweyes; class l r d;
   model count=r l d/dist=poisson type3; run;
/**************************************************************/
/*Mark E:  Polychoric Correlation.                            */
/**************************************************************/
title 'Polychoric Correlation';
proc freq data=eyes; weight count; 
   table r*l/noprint plcorr; run;