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