1# Copyright 20014 by Roger Bivand, Virgilio Gómez-Rubio
2#
3
4
5lee.mc <- function(x, y, listw, nsim, zero.policy=NULL,
6	alternative="greater", na.action=na.fail, spChk=NULL,
7        return_boot=FALSE) {
8	alternative <- match.arg(alternative, c("greater", "less"))
9	if(!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
10		"is not a listw object"))
11	if(!is.numeric(x) | !is.numeric(y)) stop(paste(deparse(substitute(x)),
12		"is not a numeric vector"))
13        if (is.null(zero.policy))
14            zero.policy <- get("zeroPolicy", envir = .spdepOptions)
15        stopifnot(is.logical(zero.policy))
16	if(missing(nsim)) stop("nsim must be given")
17	if (is.null(spChk)) spChk <- get.spChkOption()
18	if (spChk && !chkIDs(x, listw) && !chkIDs(y, listw))
19		stop("Check of data and weights ID integrity failed")
20	cards <- card(listw$neighbours)
21	if (!zero.policy && any(cards == 0))
22		stop("regions with no neighbours found")
23#	if (any(is.na(x))) stop("NA in X")
24#	if (any(is.na(y))) stop("NA in Y")
25	xname <- deparse(substitute(x))
26	yname <- deparse(substitute(y))
27	wname <- deparse(substitute(listw))
28	if (deparse(substitute(na.action)) == "na.pass")
29	    stop("na.pass not permitted")
30
31	#Check NA's in both vectors
32	na.act <- attr(na.action(cbind(x,y)), "na.action")
33
34	x[na.act]<-NA
35	y[na.act]<-NA
36
37	x<-na.action(x)
38	y<-na.action(y)
39
40	if (!is.null(na.act)) {
41	    subset <- !(1:length(listw$neighbours) %in% na.act)
42	    listw <- subset(listw, subset, zero.policy=zero.policy)
43	}
44	n <- length(listw$neighbours)
45	if ((n != length(x)) |(n != length(y))) stop("objects of different length")
46        gamres <- suppressWarnings(nsim > gamma(n + 1))
47        if (gamres) stop("nsim too large for this number of observations")
48	if (nsim < 1) stop("nsim too small")
49
50#	S0 <- Szero(listw)
51	S2<-sum ( (unlist(lapply(listw$weights, sum)))^2 )
52
53	#Data frame with x, y
54	xy<-data.frame(x,y)
55        if (return_boot) {
56            lee_boot <- function(var, i, ...) {
57#                var <- var[i]
58#                return(moran(x=var, ...)$I)
59		return(lee(x=var[i,1], y=var[i,2], ...)$L)
60            }
61            cores <- get.coresOption()
62            if (is.null(cores)) {
63            parallel <- "no"
64            } else {
65                parallel <- ifelse (get.mcOption(), "multicore", "snow")
66            }
67            ncpus <- ifelse(is.null(cores), 1L, cores)
68            cl <- NULL
69            if (parallel == "snow") {
70                cl <- get.ClusterOption()
71                if (is.null(cl)) {
72                    parallel <- "no"
73                    warning("no cluster in ClusterOption, parallel set to no")
74                }
75            }
76            res <- boot(xy, statistic=lee_boot, R=nsim,
77                sim="permutation", listw=listw, n=n, S2=S2,
78                zero.policy=zero.policy, parallel=parallel, ncpus=ncpus, cl=cl)
79            return(res)
80        }
81	res <- numeric(length=nsim+1)
82	for (i in 1:nsim)
83	{
84		idx<-sample(1:n)
85		res[i] <- lee(x[idx], y[idx], listw, n, S2,
86	    zero.policy)$L
87	}
88	res[nsim+1] <- lee(x, y, listw, n, S2, zero.policy)$L
89	rankres <- rank(res)
90	xrank <- rankres[length(res)]
91	diff <- nsim - xrank
92	diff <- ifelse(diff > 0, diff, 0)
93	if (alternative == "less")
94        	pval <- punif((diff + 1)/(nsim + 1), lower.tail=FALSE)
95    	else if (alternative == "greater")
96        	pval <- punif((diff + 1)/(nsim + 1))
97	if (!is.finite(pval) || pval < 0 || pval > 1)
98		warning("Out-of-range p-value: reconsider test arguments")
99	statistic <- res[nsim+1]
100	names(statistic) <- "statistic"
101	parameter <- xrank
102	names(parameter) <- "observed rank"
103	method <- "Monte-Carlo simulation of Lee's L"
104	data.name <- paste(
105	    xname, ", ", yname, "\nweights:",
106	    wname, ifelse(is.null(na.act), "", paste("\nomitted:",
107	    paste(na.act, collapse=", "))),
108"\nnumber of simulations + 1:",
109	    nsim+1, "\n")
110	lres <- list(statistic=statistic, parameter=parameter,
111	    p.value=pval, alternative=alternative, method=method,
112	    data.name=data.name, res=res)
113	if (!is.null(na.act) ) attr(lres, "na.action") <- na.act
114	class(lres) <- c("htest", "mc.sim")
115	lres
116}
117
118