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