postscript("l05.ps");options(width=68)
#********************************************************/
# Mark.A Example of McNemar's test.                     */
#********************************************************/
#***********************************************************/
# Purity data, example H of Cox and Snell (1981) Applied   */
# Statistics.  Two manufacturing processes were applied to */
# raw materials from 22 batches, yeilding 22 matched pairs.*/
# First variable is covariate measuring purity of batch,   */
# and next two variables reflect whether the result was    */
# flawed.  Does one of the processes produce fewer flaws?  */
#***********************************************************/
purity<-as.data.frame(scan("purity.dat",
  what=list(pure=0,r1="",r2="")))
cat('\n McNemar test  \n')
mcnemar.test(purity$r1,purity$r2,F)
purity$batch<-seq(length(purity$r1))
puritya<-purity;puritya$proc<-1;puritya$r<-puritya$r1
purityb<-purity;purityb$proc<-2;purityb$r<-purityb$r2
purity2<-rbind(puritya,purityb)
mantelhaen.test(xtabs(~proc+r+batch,data=purity2),correct=F)
#*******************************************************/
# Mark.B Create three fake data sets, one consistent   */
# with H0 for CMH and McNemar, and two consistent with */
# one H0 but not with the other.                       */
#*******************************************************/
cat('\n Fake comparisons of McNemar and Pearson  \n')
fake<-as.data.frame(list(a=rep(c(0,0,1,1),3),
  b=rep(c(0,1),6),ex=rep(1:3,rep(4,3)),
  we=c(100,50,50,25,100,50,50,0,200,50,100,25)))
for(i in 1:3){
  tab<-xtabs(we~a+b,subset=(ex==i),data=fake)
  print(tab);print(chisq.test(tab,correct=F))
  print(mcnemar.test(tab,correct=F))
}
#*******************************************************/
# Mark.C Poisson regression example.                   */
#*******************************************************/
#******************************************************/
# Number of accidental deaths from falls, 1979        */
# World Almanac and Book of Facts 1984 (New York:     */
# Newspaper Enterprise Association, 1984), and quoted */
# in the Minitab Handbook, 3rd Edn. (Belmont, CA: Dux-*/
# bury, 1994).  Falls are by month.  Columns are num- */
# bers of falls, numbers of days in the month, and    */
# month abbreviation.  Deaths may be influence by the */
# weather; scorebelow scores months based on icyness. */
#******************************************************/
falls<-as.data.frame(scan("falls.dat",  
   what=list(falls=0,days=0,month=""))) 
falls$score<-(falls$month%in%c("Jan","Feb"))+         
   2*(falls$month%in%c("Nov","Dec","Mar"))            
#******************************************************/
# Fit a Poisson regression model, without and with    */
# accounting for number of days in each month.        */
#******************************************************/
cat('\n Poisson regression model  \n')
go1<-glm(falls~factor(month),data=falls,family=poisson)
cat('\n Poisson regression model with offset  \n')
go2<-glm(falls~offset(log(days))+factor(month),
   data=falls,family=poisson)
anova(go2)
cat('\n Chi-square approach  \n')
#********************************************************/
# Compare the Poisson regression approach to the Pearson*/
# chi-square approach.                                  */
#********************************************************/
chisq.test(falls$falls,p=falls$days/sum(falls$days))