WallerLawson.test<-function(oo, ee, g, Nrep){ # ----------------------------------------------------------------- # An example of execution # # d<-1:10; g<- 1/d; Nrep<-999 # ob<- c(0, 1, 3, 2, 1, 8, 7, 4, 0, 2) # ex<- c(0.041,0.701,1.128,1.634,1.868,5.738,4.550,3.617,3.792,4.931) # out<-WallerLawson.test(oo,ee,g,Nrep) # # --------------------------------------------------------------------- # ARGUEMENTS # oo : a vector of observed number of cases sorted in the ascending # order of the distance to the point source # ee : a vector of conditional expected number of cases sorted # in the ascending order of the distance to the point source # g : a vector of exposure function of distance to the point source # Nrep : The number of Monte Carlo replications, e.g., 999, 9999 # # VALUES # Test.statistic : test statistic # p.value : simulated p-value # Freq : a vector of simulated test statistics under the null # #----------------------------------------------------------------------- p<-ee/sum(ee) q<-oo/sum(oo) nn<-sum(oo) nr<-length(oo) h1<- sum( g* (q-p) ) * nn h2<- ( sum( g*g * p ) - sum( g * p )^{2} ) * nn Twl<- h1 / sqrt( h2 ) # p<-ee/sum(ee) nloop<-Nrep+1 for (i in 2:nloop){ q<-rmultinom(1,nn,p) / nn h1<- sum( g* (q-p) ) * nn h2<- ( sum( g*g * p ) - sum( g * p )^{2} ) * nn Twl1<- h1/ sqrt( h2 ) Twl<- c(Twl, Twl1) } p1<- rank(-Twl)[1]/(Nrep+1) list(Test.statistic=Twl[1], p.value=p1, Freq=Twl) }