DiggleChetwynd.test<-function(case,control,polygon,s,Nrep){ # ----------------------------------------------------------------- # An example of execution # # See section 6.5 for details. # # --------------------------------------------------------------------- # ARGUEMENTS # case$x : a vector of x-coordinates of case # case$y : a vector of y-coordinates of case # cont$x : a vector of x-coordinates of control # cont$y : a vector of y-coordinates of control # s : a vector of theta's values at which to calculate the # test statistic # Nrep : The number of Monte Carlo replications, e.g., 999, 9999 # # VALUES # Freq.Dsum : a vector of simulated values of Dsum # Freq.Dmax : a vector of simulated values of Dmax # Sum.p.value : simulated p-value for Dsum # Max.p.value : simulated p-value for Dmax #----------------------------------------------------------------------- nc<-length(case)/2 n<-(length(case)+length(control))/2 ind<-1:n xdat<-c(case[,1],control[,1]) ydat<-c(case[,2],control[,2]) Kcase<- khat(case, polygon, s) Kcont<- khat(control, polygon, s) SEcc <- secal(case, control, polygon, s) matplot(s, 2*cbind(SEcc, -SEcc), type='l', lty=3, col=1) lines(s, Kcase-Kcont, lwd=2) D1<- sum( (Kcase-Kcont)/SEcc ) D2<- max( (Kcase-Kcont)/SEcc ) for (i in 2:Nrep+1){ Sind<-sample(ind) Scases<-as.points(xdat[Sind<=nc],ydat[Sind<=nc]) Sconts<-as.points(xdat[Sind>=nc+1],ydat[Sind>=nc+1]) Kcase<- khat(Scases, polygon, s) Kcont<- khat(Sconts, polygon, s) SEcc <-secal(Scases, Sconts, polygon, s) sim1<- sum( (Kcase-Kcont)/SEcc ) sim2<- max( (Kcase-Kcont)/SEcc ) D1<-c(D1,sim1) D2<-c(D2,sim2) } p1<-rank(-D1)[1]/(Nrep+1) p2<-rank(-D2)[1]/(Nrep+1) list(Freq.Dsum=D1, Freq.Dmax=D2, Sum.p.value=p1, Max.p.value=p2) }