postscript("l02.ps");options(width=68)
#*******************************************************/
#          Lotto Example                               */
# CA lotto data from Friedman Pisani and Purvis (1998).*/
# Only one set of balls is used during each lottery    */
# drawing.  Columns represent sets of balls.  Rows rep-*/
# resent numbers 0-9.  Entries in table represent count*/
# of how often each ball appears.  Each set of balls   */
# was used 120 times, and so the expected number of    */
# times each ball comes up is 12.                      */
#*******************************************************/
lotto<-read.table("lotto.dat")                       
names(lotto)<-c("seta","setb","setc","setd")         
#***********************************************************/
# Mark.A Perform the chi-square test of the null hypothesis*/
# that each ball is equally likely on each draw.           */ 
#***********************************************************/
# The data step calculates the contribution for each ball.*/
lotto$setad<-(lotto$seta-12)^2/12  
lotto$setbd<-(lotto$setb-12)^2/12 
lotto$setcd<-(lotto$setc-12)^2/12 
lotto$setdd<-(lotto$setd-12)^2/12 
lotto        
# apply adds the contributions by set, only for those vari-
# ables containing squared differences.
chisq<-apply(as.matrix(lotto[,
     c("setad","setbd","setcd","setdd")]),2,"sum") 
# Use the chi-square cdf function to get the p value.  Note
# that lower.tail=F gives the upper tail.  Lower tail is
# the default.  Add this column of p-values to the existing
# vector of statistic values.
chisq<-rbind(chisq,pchisq(chisq,9,lower.tail=F)) 
dimnames(chisq)[[1]]<-c("Statistic value","p-value") 
print(chisq)
cat('\n Using built-in routines for goodness of fit test  \n')
#The univarite goodness of fit tests takes equally-likely
#categories as the default, and so the two following con-
#structions are equivalent.
chisq.test(lotto$seta) 
chisq.test(lotto$seta,p=rep(.1,10)) 
#******************************************************/
# 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"))            
#*******************************************************/
# Test to see whether expected number of falls is prop-*/
# ortional to the number of days in the month.  In this*/
# case the expected number of deaths from falls in a   */
# given month, under the null hypothesis, is the total */
# number of deaths times the number of days in the     */
# month divided by the number of days in the year.     */
#*******************************************************/
falls$exp<-falls$days*sum(falls$falls)/sum(falls$days)
tests<-list(                                          
   pears=sum((falls$falls-falls$exp)^2/falls$exp))    
tests$pearspv<-pchisq(tests$pears,                    
   length(falls$falls)-1,lower.tail=F)                
tests                                                 
chisq.test(falls$falls,p=falls$days/sum(falls$days))  
#***********************************************************/
# Mark.B Perform the test of trend on falls.               */
#***********************************************************/
tests<-falls; tests$score<-c(1,1,2,0,0,0,0,0,0,0,2,2)
trendtests<-list(                                     
   trend=sum(tests$score*(tests$falls-tests$exp)),    
   v1=sum(tests$score^2*tests$exp),                   
   v2=sum(tests$score*tests$exp))                     
trendtests$trend<-trendtests$trend/                   
  sqrt(trendtests$v1-trendtests$v2^2/sum(falls$falls))
trendtests$trendpv<-pchisq(trendtests$trend^2,1,      
                                      lower.tail=F)   
trendtests