DiggleETAL.test<-function(pts,time,polygon,range,s,t,Nrep){ # ----------------------------------------------------------------- # An example of execution # # dat<-scan("d:/demo/McHardy.dat",list(x=0,y=0,time=0)) # pts<- as.points(dat$x, dat$y) # kap.poly<-spoints(scan("d:/book/DiseaseCluster/kappoly.txt")) # (The file "kappoly.txt" contains a set of points for a polygon) # range<-c(1955, 1980); s<- seq(1,5); t<- seq(0, 4)) # Nrep<-999 # out<-DiggleETAL.test(pts,dat$t,polygon,range,s,t,Nrep) # # --------------------------------------------------------------------- # ARGUEMENTS # pts : a sets of points (x, y) as defined in Splancs # time : a vector of times, the same length of as the number of points in pts # polygon : a polygon enclosing the points # range : a vector of length 2 specifying the lower and # upper temporal domain # s : a vector of equally spaced spatial distances for the analysis # t : a vector of equally spaced temporal distances for the analysis # Nrep : The number of Monte Carlo replications, e.g., 999, 9999 # # VALUES # Sum.p.value : simulated p-value for the test statistic Usum # Max.p.value : simulated p-value for the test statistic Umax # Freq.Usum : a vector of simulated test statistics Usum under the null # Freq.Umax : a vector of simulated test statistics Umax under the null # #----------------------------------------------------------------------- kap1 <- stkhat(pts, time, polygon, range, s, t) se <- stsecal(pts, time, polygon, range, s, t) Res <- turD / kap1se D1<- sum( Res ) D2<- max( Res ) for (i in 2:Nrep+1){ Stime <-sample(time) Skap1 <- stkhat(pts, Stime, polygon, range,s,t) f1<-matrix(Skap1$ks) f2<-matrix(Skap1$kt) f1f2 <- f1 %*% t(f2) SturD<-Skap1$kst - f1f2 Sse <- stsecal(pts, Stime, polygon, range, s, t) SRes <- SturD / Sse sim1<- sum( SRes ) sim2<- max( SRes ) D1<-c(D1,sim1) D2<-c(D2,sim2) } p1<-rank(-D1)[1]/(Nrep+1) p2<-rank(-D2)[1]/(Nrep+1) list(Sum.p.value=p1, Max.p.value=p2, Freq.Usum=D1, Freq.Umax=D2) }