1# Copyright 2002-2008 by Roger Bivand and Michael Tiefelsdorf,
2# with contributions by Danlin Yu
3#
4
5localmoran.sad <- function (model, select, nb, glist = NULL, style = "W",
6    zero.policy = NULL, alternative = "greater", spChk=NULL,
7    resfun=weighted.residuals,
8    save.Vi = FALSE, tol = .Machine$double.eps^0.5,
9    maxiter = 1000, tol.bounds=0.0001, save.M=FALSE, Omega=NULL) {
10# need to impose check on weights TODO!!
11# class to inherits Jari Oksanen 080603
12    if (!inherits(nb, "nb"))
13        stop(paste(deparse(substitute(nb)), "not an nb object"))
14        if (is.null(zero.policy))
15            zero.policy <- get("zeroPolicy", envir = .spdepOptions)
16        stopifnot(is.logical(zero.policy))
17    n <- length(nb)
18    dmc <- deparse(model$call)
19    if (!inherits(model, "lm"))
20     	stop(paste(deparse(substitute(model)), "not an lm object"))
21    u <- resfun(model)
22    if (n != length(u))
23        stop("objects of different length")
24    if (is.null(spChk)) spChk <- get.spChkOption()
25    if (spChk && !chkIDs(u, nb2listw(nb, zero.policy=zero.policy)))
26	stop("Check of data and weights ID integrity failed")
27    if (!(alternative %in% c("greater", "less", "two.sided")))
28	stop("alternative must be one of: \"greater\", \"less\", or \"two.sided\"")
29    if (missing(select)) select <- 1:n
30    u <- as.vector(u)
31    select <- unique(as.integer(select))
32    if (length(select) < 1L) stop("select too short")
33    if (any(select < 1) || any(select > n))
34        stop("select out of range")
35    utu <- c(t(u) %*% u)
36    p <- model$rank
37    p1 <- 1:p
38    nacoefs <- which(is.na(coefficients(model)))
39    m <- n - p - 2
40    XtXinv <- chol2inv(model$qr$qr[p1, p1, drop = FALSE])
41    X <- model.matrix(terms(model), model.frame(model))
42# fixed after looking at TOWN dummy in Boston data
43    if (length(nacoefs) > 0L) X <- X[,-nacoefs]
44    if (!is.null(wts <- weights(model))) {
45	X <- sqrt(diag(wts)) %*% X
46    }
47    cond.sad <- FALSE
48    if (!is.null(Omega)) {
49        Omega <- chol(Omega)
50        M <- diag(n) - X %*% tcrossprod(XtXinv, X)
51        M1 <- Omega %*% M
52        M2 <- M %*% t(Omega)
53        cond.sad <- TRUE
54    }
55    B <- listw2U(nb2listw(nb, glist=glist, style="B",
56	zero.policy=zero.policy))
57    D <- NULL
58    a <- NULL
59    if (style == "W") {
60        D <- 1/sapply(B$weights, sum)
61    } else if (style == "S") {
62        D <- 1 / sqrt(sapply(B$weights, function(x) sum(x^2)))
63#        a <- sum(unlist(B$weights))
64# correction by Danlin Yu, 25 March 2004
65	a <- sum(sapply(B$weights, function(x) sqrt(sum(x^2))))
66    } else if (style == "C") a <- sum(unlist(B$weights))
67    res <- vector(mode="list", length=length(select))
68    for (i in 1:length(select)) {
69        Vi <- listw2star(B, select[i], style=style, n, D, a,
70	    zero.policy=zero.policy)
71        Viu <- lag.listw(Vi, u, zero.policy=TRUE)
72	Ii <- c((t(u) %*% Viu) / utu)
73	if (cond.sad) {
74            obj <- sadLocalMoranAlt(Ii, Vi, M1, M2, n, tol.bounds,
75                tol, maxiter, ii=select[i], alternative=alternative)
76            sad.p <- obj$sad.p
77            sad.r <- obj$sad.r
78            sad.u <- obj$sad.u
79            omega <- obj$omega
80            p.sad <- obj$p.sad
81            gamma <- obj$gamma
82	} else {
83	    ViX <- lag.listw(Vi, X, zero.policy=TRUE)
84	    MViM <- t(X) %*% ViX %*% XtXinv
85	    t1 <- -sum(diag(MViM))
86	    sumsq.Vi <- function(x) {
87                if (is.null(x)) NA
88	        else sum(x^2)
89	    }
90	    trVi2 <- sum(sapply(Vi$weights, sumsq.Vi), na.rm=TRUE)
91	    t2a <- sum(diag(t(ViX) %*% ViX %*% XtXinv))
92	    t2b <- sum(diag(MViM %*% MViM))
93	    t2 <- trVi2 - 2*t2a + t2b
94	    e1 <- 0.5 * (t1 + sqrt(2*t2 - t1^2))
95	    en <- 0.5 * (t1 - sqrt(2*t2 - t1^2))
96            gamma <- c(c(e1), c(en))
97            obj <- sadLocalMoran(Ii, gamma, m, ii=select[i],
98                alternative=alternative)
99            sad.p <- obj$sad.p
100            sad.r <- obj$sad.r
101            sad.u <- obj$sad.u
102            omega <- obj$omega
103            p.sad <- obj$p.sad
104            gamma <- obj$gamma
105	}
106        statistic <- sad.p
107        attr(statistic, "names") <- "Saddlepoint approximation"
108        p.value <- p.sad
109        estimate <- c(Ii)
110        attr(estimate, "names") <- "Observed Moran Ii"
111        internal1 <- c(omega, sad.r, sad.u)
112        attr(internal1, "names") <- c("omega", "sad.r", "sad.u")
113        method <- paste("Saddlepoint approximation for local Moran I",
114            "(Barndorff-Nielsen formula)")
115        data.name <- paste("region:", select[i],
116	    attr(nb, "region.id")[select[i]],
117	    "\n", paste(strwrap(paste("model: ", gsub("[ ]+", " ",
118	    paste(dmc, sep="", collapse="")))),
119	    collapse="\n"),
120            "\nneighbours:", deparse(substitute(nb)),
121	    "style:", style, "\n")
122        obj <- list(statistic = statistic, p.value = p.value,
123            estimate = estimate, method = method,
124	    alternative = alternative, data.name = data.name,
125	    internal1 = internal1, df = (n-p), tau = gamma,
126	    i = paste(select[i], attr(nb, "region.id")[select[i]]),
127#	    if (save.Vi) {Vi = Vi}
128	    Vi = if(save.Vi) Vi else NULL)
129        class(obj) <- "moransad"
130	res[[i]] <- obj
131    }
132    class(res) <- "localmoransad"
133    if (save.M && cond.sad) attr(res, "M") <- list(M1=M1, M2=M2, type="cond")
134    if (save.M && !cond.sad) attr(res, "M") <- list(X=X, XtXinv=XtXinv,
135        type="null")
136    res
137}
138
139sadLocalMoranAlt <- function(Ii, Vi, M1, M2, n, tol.bounds=0.0001,
140    tol = .Machine$double.eps^0.5, maxiter = 1000, ii, alternative="greater") {
141    ViI <- listw2mat(Vi) - Ii * diag(n)
142    innerTerm <- M1 %*% ViI %*% M2
143    evalue <- eigen(innerTerm, only.values=TRUE)$values
144    tau <- c(evalue)
145    e1 <- tau[1]
146    en <- tau[length(tau)]
147    low <- (1 / (2*tau[length(tau)])) + tol.bounds #+ 0.01
148    high <- (1 / (2*tau[1])) - tol.bounds #- 0.01
149    f <- function(omega, tau) {sum(tau/(1 - (2*omega*tau)))}
150    root <- uniroot(f, lower=low, upper=high, tol=tol, maxiter=maxiter,
151      	tau=tau)
152    omega <- root$root
153# 0 should be expectation - maybe use try()
154    if (omega < 0 ) sad.r <- try(-sqrt(sum(log(1 - 2*omega*tau))))
155    else sad.r <- try(sqrt(sum(log(1 - 2*omega*tau))))
156    if (inherits(sad.r, "try.error")) {
157    	warning (paste("In zone:", ii, "sad.r not a number"))
158        sad.r <- sad.u <- sad.p <- NaN
159    } else {
160	sad.u <- omega * sqrt(2*sum(tau^2 / (1 - (2*omega*tau))^2))
161    	sad.p <- sad.r - ((1/sad.r)*log(sad.r/sad.u))
162    }
163        if (alternative == "two.sided") p.sad <- 2 * pnorm(abs(sad.p),
164	    lower.tail=FALSE)
165        else if (alternative == "greater")
166            p.sad <- pnorm(sad.p, lower.tail=FALSE)
167        else p.sad <- pnorm(sad.p)
168    obj <- list(p.sad=p.sad, sad.p=sad.p, sad.r=sad.r, sad.u=sad.u,
169        omega=omega, root=root, gamma=tau)
170    obj
171}
172
173sadLocalMoran <- function(Ii, gamma, m, ii, alternative="greater") {
174	e1 <- gamma[1]
175	en <- gamma[2]
176	l <- en
177	h <- e1
178	mi <- Ii
179	aroot= m*mi*(l+h-2*mi)+mi*(3*l+3*h-4*mi)-2*l*h
180        broot= (m+2)*mi*(l-mi)*(h-mi)
181        c1root= l**2 * mi**2 * (m+1)**2 + h**2 * mi**2 * (m+1)**2
182        c2root= 2*l*h * (2*l*h - 2*l*mi - 2*h*mi - 2*m*mi**2 -
183	    m**2 * mi**2 + mi**2)
184        omega= 0.25*((aroot-sqrt(c1root+c2root))/broot)
185	if (is.nan(omega)) {
186	    warning (paste("In zone:", ii, "omega not a number"))
187	    sad.r <- sad.u <- sad.p <- NaN
188	} else {
189            tau <- c(c(e1), rep(0, m), c(en))
190	    taumi <- tau - Ii
191            if (omega < 0 ) sad.r <- -sqrt(sum(log(1 - 2*omega*taumi)))
192            else sad.r <- sqrt(sum(log(1 - 2*omega*taumi)))
193            sad.u <- omega * sqrt(2*sum(taumi^2 / (1 - (2*omega*taumi))^2))
194            sad.p <- sad.r - ((1/sad.r)*log(sad.r/sad.u))
195	}
196        if (alternative == "two.sided") p.sad <- 2 * pnorm(abs(sad.p),
197	    lower.tail=FALSE)
198        else if (alternative == "greater")
199            p.sad <- pnorm(sad.p, lower.tail=FALSE)
200        else p.sad <- pnorm(sad.p)
201	obj <- list(p.sad=p.sad, sad.p=sad.p, sad.r=sad.r, sad.u=sad.u,
202            omega=omega, gamma=gamma)
203        obj
204}
205
206print.localmoransad <- function(x, ...) {
207    extract <- function(x, i) {x[[i]]}
208    regnames <- sapply(x, extract, 10)
209    est <- sapply(x, extract, 3)
210    sad <- sapply(x, extract, 1)
211    pval <- sapply(x, extract, 2)
212    res <- as.matrix(cbind(est, sad, pval))
213    rownames(res) <- regnames
214    colnames(res) <- c("Local Morans I", "Saddlepoint", "Pr. (Sad)")
215    print(res, ...)
216    invisible(res)
217}
218
219as.data.frame.localmoransad <- function(x, row.names=NULL, optional=FALSE, ...) {
220    n <- length(x)
221    if (n < 1) stop("x too short")
222    res <- matrix(0, nrow=n, ncol=14)
223    regnames <- NULL
224    if (!is.null(row.names))
225	if (length(row.names) == n) regnames <- row.names
226    if (is.null(regnames))for (i in 1:n) regnames <- c(regnames, x[[i]]$i)
227    for (i in 1:n) {
228        tau <- x[[i]]$tau
229	df <- x[[i]]$df
230        tau <- c(tau[1], rep(0, df-2), tau[2])
231        max.I <- tau[1]
232        min.I <- tau[length(tau)]
233        E.I <- sum(tau)/df
234        tau <- tau - E.I
235        V.I <- (2*sum(tau^2)) / (df*(df+2))
236        Z.I <- (x[[i]]$estimate - E.I) / sqrt(V.I)
237	if (x[[i]]$alternative == "two.sided")
238	    P.I <- 2 * (1 - pnorm(Z.I))
239        else if (x[[i]]$alternative == "greater")
240            P.I <- pnorm(Z.I, lower.tail=FALSE)
241        else P.I <- pnorm(Z.I)
242        Sk.I <- ((8*sum(tau^3))/(df*(df+2)*(df+4))) / (V.I^(3/2))
243        Kur.I <- ((48*sum(tau^4) + 12*(sum(tau^2))^2) /
244            (df*(df+2)*(df+4)*(df+6))) / (V.I^2)
245	res[i,] <- c(x[[i]]$estimate, Z.I, P.I, x[[i]]$statistic,
246	    x[[i]]$p.value, E.I, V.I, Sk.I, Kur.I, min.I, max.I,
247	    x[[i]]$internal1)
248    }
249    colnames(res) <- c("Local Morans I", "Stand. dev. (N)", "Pr. (N)",
250        "Saddlepoint", "Pr. (Sad)", "Expectation", "Variance",
251        "Skewness", "Kurtosis", "Minimum", "Maximum",
252        "omega", "sad.r", "sad.u")
253    rownames(res) <- regnames
254    res <- as.data.frame(res)
255    res
256}
257
258summary.localmoransad <- function(object, ...) {
259    res <- as.data.frame(object)
260    class(res) <- c("summary.localmoransad", class(res))
261    res
262}
263
264print.summary.localmoransad <- function(x, ...) {
265	print(as.data.frame(x), ...)
266	invisible(x)
267}
268
269
270
271