1# Copyright (c) 2007-2008 Markus Reder and Roger Bivand
2
3localmoran.exact.alt <- function(model, select, nb, glist = NULL, style = "W",
4    zero.policy = NULL, alternative = "greater", spChk=NULL,
5    resfun=weighted.residuals, Omega=NULL, save.Vi = FALSE, save.M=FALSE,
6    useTP=FALSE, truncErr=1e-6, zeroTreat=0.1) {
7# need to impose check on weights TODO!!
8# class to inherits Jari Oksanen 080603
9    if (!inherits(nb, "nb"))
10        stop(paste(deparse(substitute(nb)), "not an nb object"))
11        if (is.null(zero.policy))
12            zero.policy <- get("zeroPolicy", envir = .spdepOptions)
13        stopifnot(is.logical(zero.policy))
14#    if (class(model) != "lm")
15#        stop(paste(deparse(substitute(model)), "not an lm object"))
16    dmc <- deparse(model$call)
17    n <- length(nb)
18    if (!inherits(model, "lm"))
19     	stop(paste(deparse(substitute(model)), "not an lm object"))
20    if (is.null(Omega)) Omega <- diag(n)
21    else {
22        if (dim(Omega)[1] != n) stop("Omega of different size than data")
23        Omega <- chol(Omega)
24    }
25    u <- resfun(model)
26    if (n != length(u))
27        stop("objects of different length")
28    if (is.null(spChk)) spChk <- get.spChkOption()
29    if (spChk && !chkIDs(u, nb2listw(nb, zero.policy=zero.policy)))
30	stop("Check of data and weights ID integrity failed")
31    if (!(alternative %in% c("greater", "less", "two.sided")))
32	stop("alternative must be one of: \"greater\", \"less\", or \"two.sided\"")
33    if (missing(select)) select <- 1:n
34    u <- as.vector(u)
35    select <- unique(as.integer(select))
36    if (length(select) < 1L) stop("select too short")
37    if (any(select < 1) || any(select > n))
38        stop("select out of range")
39    utu <- c(crossprod(u))
40    p <- model$rank
41    p1 <- 1:p
42    nacoefs <- which(is.na(coefficients(model)))
43    m <- n - p - 2
44    XtXinv <- chol2inv(model$qr$qr[p1, p1, drop = FALSE])
45    X <- model.matrix(terms(model), model.frame(model))
46# fixed after looking at TOWN dummy in Boston data
47    if (length(nacoefs) > 0L) X <- X[,-nacoefs]
48    if (!is.null(wts <- weights(model))) {
49	X <- sqrt(diag(wts)) %*% X
50    }
51    M <- diag(n) - X %*% tcrossprod(XtXinv, X)
52    M1 <- Omega %*% M
53    M2 <- M %*% t(Omega)
54    B <- listw2U(nb2listw(nb, glist=glist, style="B",
55	zero.policy=zero.policy))
56    D <- NULL
57    a <- NULL
58    if (style == "W") {
59        D <- 1/sapply(B$weights, sum)
60    } else if (style == "S") {
61        D <- 1 / sqrt(sapply(B$weights, function(x) sum(x^2)))
62#        a <- sum(unlist(B$weights))
63# correction by Danlin Yu, 25 March 2004
64	a <- sum(sapply(B$weights, function(x) sqrt(sum(x^2))))
65    } else if (style == "C") a <- sum(unlist(B$weights))
66    res <- vector(mode="list", length=length(select))
67    for (i in 1:length(select)) {
68        Vi <- listw2star(B, select[i], style=style, n, D, a,
69	    zero.policy=zero.policy)
70        Viu <- lag.listw(Vi, u, zero.policy=TRUE)
71	Ii <- c(crossprod(u, Viu) / utu)
72
73        obj <- exactLocalMoranAlt(Ii=Ii, Vi=Vi, M1=M1, M2=M2, n=n,
74            alternative=alternative, useTP=useTP, truncErr=truncErr,
75                zeroTreat=zeroTreat)
76        data.name <- paste("region:", select[i],
77	    attr(nb, "region.id")[select[i]],
78	    "\n", paste(strwrap(paste("model: ", gsub("[ ]+", " ",
79	    paste(dmc, sep="", collapse="")))),
80	    collapse="\n"),
81            "\nneighbours:", deparse(substitute(nb)),
82	    "style:", style, "\n")
83        obj$data.name <- data.name
84        obj$df <- (n-p)
85        obj$i <- paste(select[i], attr(nb, "region.id")[select[i]])
86        obj$Vi <- if(save.Vi) Vi else NULL
87	res[[i]] <- obj
88    }
89    class(res) <- "localmoranex"
90    if (save.M) attr(res, "M") <- list(M1=M1, M2=M2)
91    res
92}
93
94exactLocalMoranAlt <- function(Ii, Vi, M1, M2, n, alternative,
95    type="Alternative", useTP=FALSE, truncErr=1e-6, zeroTreat=0.1) {
96    ViI <- listw2mat(Vi) - Ii * diag(n)
97    innerTerm <- M1 %*% ViI %*% M2
98    evalue <- eigen(innerTerm, only.values=TRUE)$values
99    gamma <- c(evalue)
100    obj <- exactMoran(Ii, gamma, alternative=alternative,
101        type=type, useTP=useTP, truncErr=truncErr, zeroTreat=zeroTreat)
102    obj
103}
104
105