Jacquez.test<-function(x,y,time,k,Nrep){ # ----------------------------------------------------------------- # An example of execution # # dat<-scan("d:/demo/McHardy.dat",list(x=0,y=0,time=0)) # k<-1; Nrep<-999 # out<-Jacquez.test(dat$x,dat$y,dat$t,k,Nrep) # # --------------------------------------------------------------------- # ARGUEMENTS # x : a vector of x-coordinates of case # y : a vector of y-coordinates of case # time : a vector of observed times # k : k of k nearest neighbors # Nrep : The number of Monte Carlo replications, e.g., 999, 9999 # # VALUES # Jacquez.T : test statistic # Expected : expected value of test statistic # Simulated.p.value : simulated p-value # Freq : a vector of simulated test statistics under the null # #----------------------------------------------------------------------- n<-length(x) time<-time+runif(n)/100 ktmon<-1:(Nrep+1) sdis<-matrix(0,n,n) tdis<-matrix(0,n,n) as<-sdis at<-tdis for (i in 1:n){ for (j in 1:n){ sdis[i,j]<-sqrt( (x[i]-x[j])^2 + (y[i]-y[j])^2 ) tdis[i,j]<- abs( time[i] - time[j]) } } for(i in 1:n){ ff<-order(sdis[i,]) gg<-order(tdis[i,]) for(v in 1:k){ m1<-ff[v+1] m2<-gg[v+1] as[i,m1]<-1 at[i,m2]<-1 } } s1<-0; s2<-0; s3<-0 for(i in 1:n){ for (j in 1:n){ s1<-s1+as[i,j] * at[i,j] s2<-s2+as[i,j] s3<-s3+at[i,j] } } obst <- s1/2 Ms <- s2/2 Mt <- s3/2 ee <- Ms*Mt/(n*(n-1)/2) # # Monte Carlo Hypothesis testing # for(ij in 1:Nrep){ timeR<-sample(time) for (i in 1:n){ for (j in 1:n){ tdis[i,j]<- abs( timeR[i] - timeR[j]) } } at[,]<-0 for(i in 1:n){ gg<-order(tdis[i,]) for(v in 1:k){ m2<-gg[v+1] at[i,m2]<-1 } } s1<-0 for(i in 1:n){ for (j in 1:n){ s1<-s1+as[i,j] * at[i,j] } } ktmon[ij] <- s1/2 } ktmon[Nrep+1]<-obst r<-length(ktmon[ktmon>=obst]) p<-r/(Nrep+1) list(Jacquez.T=obst, Expected=ee, Simulated.p.value=p, Freq=ktmon) }