postscript("l10.ps");options(width=68)
#*****************************************************/
# 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.     */
#*****************************************************/
eyes<-as.data.frame(scan("eyes.dat",
  what=list(r=0,l=0,count=0)))
#*************************************************************/
# Mark A: Tests of Symmetry                                  */
#*************************************************************/
cat('\n Eyesight example: test symmetry  \n')
#R credits McNemar here, and not Bowker.
mcnemar.test(xtabs(count~r+l,data=eyes))
#***********************************************************/
# 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.
#***********************************************************/
elections<-as.data.frame(scan("elections.dat",
  what=list(w97="",w01="",w05="")))
elections$w05c<-elections$w05
elections$w05c[elections$w05=="RU"]<-"L"
parties<-unique(sort(as.matrix(elections)))
elections$w97<-factor(elections$w97,levels=parties)
elections$w01<-factor(elections$w01,levels=parties)
elections$w05<-factor(elections$w05,levels=parties)
elections$w05c<-factor(elections$w05c,levels=parties)
cat('\n Electon example: test symmetry  \n')
#*************************************************************/
# With zero denominators, we have a problem.                 */
#*************************************************************/
mcnemar.test(xtabs(~w97+w05,data=elections))
mcnemar.test(xtabs(~w97+w05c,data=elections))
#*************************************************************/
# 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.            */
#*************************************************************/
cat('\n Eye Symmetry via Poisson Regression  \n')
levs<-unique(c(eyes$r,eyes$l)) 
n<-length(levs) 
eyes$sym<-eyes$count; count<-0 
for(i in seq(n)) for(j in i:n){ 
  count<-count+1 
  eyes$sym[((eyes$r==levs[i])& 
     (eyes$l==levs[j]))|((eyes$r==levs[j])& 
     (eyes$l==levs[i]))]<-paste(levs[i],levs[j],sep="")
} 
#*************************************************************/
#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.  */
#*************************************************************/
glmout1<-glm(count~-1+factor(sym), data=eyes,
  family=poisson) 
fitted(glmout1) 
cat('\n Eyesight example: fit symmetry via genmod  \n')
#*************************************************************/
#Fit the saturated model, and compare the log likelihoods.   */
#Verify that the resulting test statistic is similar to      */
#Bowker's test.  Both have the same DF.                      */
#*************************************************************/
glmout2<-glm(count~factor(paste(r,l))-1,data=eyes,
   family=poisson); summary(glmout2) 
library(gmodels)
#Use the function estimable from the gmodels library to gen-
#erate a Wald test for the contrasts.                       
#R has some built-in constrast matrices.  Do ?contr.sum to  
#get a list of these; these do not include the ones we need.
#I think I wrote the ones below based on contr.sum. This is 
#more than we need at present; it is of the general form    
#that R requires for fitting the contrasts internal to glm. 
contr.both<-function(n,scores=1:n,contrasts=FALSE,
       full=TRUE,sym=TRUE){
   if (is.numeric(n) && length(n) == 1)
       levs <- 1:n
   else {
       levs <- n
       n <- length(levs)
   }
   rtn<-sqrt(n)
   contr<-if(!full){
      array(0,c(n,(if(sym){rtn}else{rtn-2})*(rtn-1)/2))
   }else{
      array(0,c(n,n-1))
   }
   count<-0
   a<-if(sym){1}else{2}
   for(i in a:(rtn-1)) for(j in (i+1):rtn){
      contr[(i-1)*rtn+j,count<-count+1]<-1
      contr[(j-1)*rtn+i,count]<- -1
      if(!sym){ 
         contr[(1-1)*rtn+j,count]<- -1 
         contr[(i-1)*rtn+1,count]<- -1 
         contr[(1-1)*rtn+i,count]<- 1 
         contr[(j-1)*rtn+1,count]<- 1 
      } 
   }
   if(full){
     for(i in a:(rtn-1)) for(j in (i+1):rtn){
        contr[(i-1)*rtn+j,count<-count+1]<-1
        contr[(j-1)*rtn+i,count]<- 1
        contr[(j-1)*rtn+j,count]<- -1
        contr[(i-1)*rtn+i,count]<- -1
     }
     for(i in 1:(rtn-1)){
        contr[(i-1)*rtn+i,count<-count+1]<-1
        contr[1,count]<- -1
     }
  }
  return(contr)
}
contr.sym<-function(n,scores=1:n,contrasts=FALSE,full=TRUE){
  return(contr.both(n,scores=scores,contrasts=contrasts,
  full=full,sym=TRUE))}
contr.qsym<-function(n,scores=1:n,contrasts=FALSE,full=TRUE){
  return(contr.both(n,scores=scores,contrasts=contrasts,
  full=full,sym=FALSE))}
#Test symmetry.  I get a warning here that seems safe to
#ignore.                                                
estimable(glmout2,t(contr.sym(16,full=F)),joint.test=T)
cat('\n Eyesight example: fit sat. model to test symmetry  \n')
#*************************************************************/
# Mark B: Quasisymmetry                                      */
#*************************************************************/
#Test quasisymmetry:
estimable(glmout2,t(contr.qsym(16,full=F)),joint.test=T)
cat('\n Eyesight example: fit sat. model to test quasisym.  \n')
cat('\n Quasisymmetry via CATMOD title2;  \n')
cat('\n Marginal symmetry via CATMOD  \n')
#*************************************************************/
#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.                                    */
#*************************************************************/
cat('\n Quasiindependence via GENMOD  \n')
eyes$qi<-eyes$r*(eyes$r==eyes$l)
glmout3<-glm(count~factor(r)+factor(l)+factor(qi),
   data=eyes, family=poisson); summary(glmout3) 
#*************************************************************/
#Mark E:  Polychoric Correlation.                            */
#*************************************************************/
cat('\n Polychoric Correlation  \n')
library(polycor)#For polychor
polychor(xtabs(count~r+l,data=eyes))
library(psych)#For cohen.kappa
cohen.kappa(xtabs(count~r+l,data=eyes))