pdf("g10.pdf");options(width=64)
#setwd("C:\\Users\\kolassa\\Class553")
setwd("~/Taught1/960-553/Data")
<pre>
#<a name=block1> Block 1</a>
#**********************************************************/
# Stuart (1953 Biometrika 40, 105-110), quoted by Stuart */
# (1955 Biometrika 42, 412-416), and Agresti (1983 Bio-  */
# metrics 39(2) 505-510), present an ordinal categoriza- */
# tion of right and left eye vision of Royal Ordinance   */
# female factory workers 1943-46.                        */
#*********************************************************/
eyes<-read.table("eyes.dat",col.names=c("r","l","count"))
#*****************************************************/
# Tests of Symmetry                                  */
#*****************************************************/
cat('\n Eyesight example: test symmetry  \n')
#R credits McNemar here, and not Bowker.
(eyetab<-xtabs(count~r+l,data=eyes))
mcnemar.test(eyetab)
#***********************************************************/
# 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<-read.table("elections.dat",
   col.names=c("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)
elections
cat('\n Electon example: test symmetry  \n')
#*************************************************************/
# With zero denominators, we have a problem.                 */
#*************************************************************/
(etab1<-xtabs(~w97+w05,data=elections))
mcnemar.test(etab1)
(etab2<-xtabs(~w97+w05c,data=elections))
mcnemar.test(etab2)
#*************************************************************/
# 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;
for(i in seq(n)) for(j in i:n){ 
  eyes$sym[((eyes$r==levs[i])& 
     (eyes$l==levs[j]))|((eyes$r==levs[j])& 
     (eyes$l==levs[i]))]<-paste(levs[i],levs[j],sep="")
} 
eyes
#*************************************************************/
#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.  */
#*************************************************************/
symmetricfit<-glm(count~-1+factor(sym), data=eyes,
  family=poisson) 
symmetricfit
fitted(symmetricfit) 
#*************************************************************/
#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.                      */
#*************************************************************/
saturated<-glm(count~factor(paste(r,l))-1,data=eyes,
   family=poisson)
summary(saturated) 
anova(symmetricfit,saturated)
#Use the function estimable from the gmodels library to gen-
#erate a Wald test for the contrasts.                       
library(gmodels)#For estimable; Fits linear constraints to generalized linear models.
#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.  First, examing contrasts we use:
t(contr.sym(16,full=F))
#I get a warning here that seems safe to ignore.
estimable(saturated,t(contr.sym(16,full=F)),joint.test=T)
#<a name=block2> Block 2</a>
#*****************************************************/
# Quasisymmetry                                      */
#*****************************************************/
#Test quasi-symmetry:
estimable(saturated,t(contr.qsym(16,full=F)),joint.test=T)
#<a name=block3> Block 3</a>
# Marginal symmetry/homogeneity test.
library(cta)#For mph.fit; Fits nonlinear constraints to contingency tables.
# h.fct is a function of cell counts or proportions, and 
# defines contrasts among proportions corresponding to 
# marginal symmetry.  M.fct also comes from library cta.
# It examines the data frame and determines which entires
# corresond to which explanatory variable levels.  Marginal 
# symmetry then is equivalent to difference in sums over 
# similar categories being zero.
h.fct<-function(p){
    p.r<-M.fct(eyes$r)%*%p
    p.l<-M.fct(eyes$l)%*%p
    as.matrix(c(p.r[-4]-p.l[-4]))
}
# Note that this calculation really can't be written in 
# format where input data is written in a data= format.
# Maybe this could have been done using attach, or with.
mph.summary(mph.fit(y=eyes$count,h.fct=h.fct))
#<a name=block4> Block 4</a>
#**********************************************************/
#Test the quasi-independence model.  Create a variable d  */
# giving the row and column number for diagonal elements. */
# Use this variable as a factor.  The baseline is 0, corr.*/
# to off-diagonal elements.     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 poisson regression  \n')
eyes$d<-eyes$r*(eyes$r==eyes$l)
quasiindfit<-glm(count~factor(r)+factor(l)+factor(d),
   data=eyes, family=poisson) 
summary(quasiindfit)
#Test quasiindependence
anova(saturated,quasiindfit)
#<a name=block5> Block 5</a>
#*****************************************************/
# Polychoric Correlation.                            */
#*****************************************************/
cat('\n Polychoric Correlation  \n')
# Note subtle spelling difference between package and
# function name.  Estimator is not MLE unless you add
# flag ML=TRUE.
library(polycor)#For polychor; Calculates polychoric correlation.
polychor(xtabs(count~r+l,data=eyes),std.err=TRUE,
   ML=TRUE)
</pre>