[R-sig-Geo] Using Spatial Eigenvector Mapping for Negative Binomial GLM
Roger Bivand
Roger.Bivand at nhh.no
Tue Jan 31 21:59:20 CET 2012
On Tue, 31 Jan 2012, Pablo Gonzalez wrote:
> Dear all,
>
> I am trying to modify the function ME() of the package spdep to reduce the
> spatial autocorrelation problem of a negative binomial model by Spatial
> Eigenvector Mapping (Dray et al 2006, Dormann et al 2007). I am having
> problems to use the eigenvectors created because they include complex
> numbers and the negative binomial function, glm.nb(), doesn't allow such
> numbers.
The usual reason for complex eigenvectors from weights matrices is that
they are asymmetric. In your editing of ME(), did you remove the line
making sure that W is symmetric:
listw <- listw2U(listw) # make weights symmetric
Try doing bits manually outside your draft function. If:
Wmat <- listw2mat(listw) # convert to full matrix form
Cent <- diag(n) - (matrix(1, n, n)/n)
eV <- eigen(Cent %*% Wmat %*% Cent, EISPACK=TRUE)$vectors
eV are complex, something odd is going on, and maybe the EISPACK=TRUE is
an issue. Are all the Im(eV) zero?
Roger
>
> Honestly, I don't know why the function eigen() create eigenvectors with
> complex numbers. It might be a problem of the spatial weight matrix which
> includes points without neighbours or the complexity itself of the matrix.
> When I use the function ME() in its original form for a glm with poisson
> family it works because the function glm() seems to accept complex
> predictors.
>
> Could be possible to force only real eigenvectors or to allow complex
> predictors in NB glms?? I appreciate any help. Thank you very much in
> advance!
>
> If it helps, I copy bellow the modified function:
>
> ME.nb <- function (formula, data, listw, weights,alpha = 0.05, nsim = 99,
> verbose = TRUE, stdev = FALSE)
> {
> #Moran's Index
> MoraneI.boot <- function(var, i, ...) {
> var <- var[i]
> I <- (n/S0) * (crossprod(sW %*% var, var))/cpvar
> return(c(as(I, "matrix")))
> }
> #Test Moran's Index
> MIR_a <- function(resids, sW, n, cpvar, S0, nsim, stdev = TRUE) {
> boot1 <- boot(resids, statistic = MoraneI.boot, R = nsim,
> sim = "permutation", sW = sW, n = n, S0 = S0, cpvar = cpvar)
> mi <- boot1$t0
> if (stdev) {
> zi <- (boot1$t0 - mean(boot1$t))/sqrt(var(boot1$t))
> pri <- pnorm(abs(zi), lower.tail = FALSE)
> }
> else {
> zi <- NA
> pri <- (sum(boot1$t >= mi) + 1)/(nsim + 1)
> }
> res <- list(estimate = mi, statistic = zi, p.value = pri)
> res
> }
> if (is.null(verbose))
> verbose <- get("verbose", env = .spdepOptions)
> stopifnot(is.logical(verbose))
> listw <- listw2U(listw)
> sW <- as_dgRMatrix_listw(listw)
> sW <- as(sW, "CsparseMatrix")
> Wmat <- listw2mat(listw)
> n <- ncol(Wmat)
> S0 <- Szero(listw)
>
> if (missing(data))
> data <- environment(formula)
> mf <- match.call(expand.dots = FALSE)
> m <- match(c("formula", "data", "weights", "offset"), names(mf),0)
> mf <- mf[c(1, m)]
> mf$drop.unused.levels <- TRUE
> mf[[1]] <- as.name("model.frame")
> mf <- eval(mf, parent.frame())
> mt <- attr(mf, "terms")
> Y <- model.extract(mf, "response")
> X <- model.matrix(mt, mf)
> weights <- model.weights(mf)
> if (!is.null(weights) && any(weights < 0))
> stop("negative weights not allowed")
> offset <- model.offset(mf)
> if (!is.null(offset) && length(offset) != NROW(Y))
> stop("number of offsets should equal number of observations")
> #Fit model
> glm_fit <- glm.nb(formula,data=data)
> glm_res <- glm_fit$y - glm_fit$fitted.values
> cpvar <- crossprod(glm_res)
> mRES <- MIR_a(glm_res, sW = sW, n = n, cpvar = cpvar, S0 = S0,
> nsim = nsim, stdev = stdev)
> pIZ <- mRES$p.value
> print(pIZ)
> tres <- c(NA, mRES$statistic, pIZ)
> if (pIZ > alpha)
> stop("base correlation larger than alpha")
> Cent <- diag(n) - (matrix(1, n, n)/n)
> #Create eigenvectors of the Spatial Weights matrix
> eV <- eigen(Cent %*% Wmat %*% Cent, EISPACK = TRUE)$vectors
> str(eV)
> rm(Cent, Wmat)
> iZ <- numeric(n)
> # Select eigenvectors that reduce Moran's Index
> for (i in 1:n) {
> iX <- eV[,i]
> i_glm <- update(glm_fit, .~.+ iX)
> glm_res <- i_glm$y - i_glm$fitted.values
> cpvar <- crossprod(glm_res)
> iZ[i] <- c(as((n/S0) * (crossprod(sW %*% glm_res,
> glm_res))/cpvar,"matrix"))
> }
> min_iZ <- which.min(abs(iZ))
> X <- eV[, min_iZ]
> glm_fit <- update(glm_fit, .~.+ X)
> glm_res <- glm_fit$y - glm_fit$fitted.values
> cpvar <- crossprod(glm_res)
> mRES <- MIR_a(glm_res, sW = sW, n = n, cpvar = cpvar, S0 = S0,
> nsim = nsim, stdev = stdev)
> pIZ <- mRES$p.value
> used <- rep(FALSE, n)
> used[min_iZ] <- TRUE
> min_v <- min_iZ
> if (verbose)
> cat("eV[,", min_iZ, "], I: ", mRES$estimate, " ZI: ",
> mRES$statistic, ", pr(ZI): ", pIZ, "\n", sep = "")
> tres <- rbind(tres, c(min_iZ, mRES$statistic, pIZ))
> # Try rest eigenvectors till Moran's index is not significant
> while (pIZ <= alpha) {
> for (i in 1:n) {
> if (!used[i]) {
> iX <- eV[,i]
> i_glm <- update(glm_fit, .~.+ iX)
> glm_res <- i_glm$y - i_glm$fitted.values
> cpvar <- crossprod(glm_res)
> iZ[i] <- c(as((n/S0) * (crossprod(sW %*% glm_res,
> glm_res))/cpvar, "matrix"))
> }
> else iZ[i] <- NA
> }
> min_iZ <- which.min(abs(iZ))
> X <- eV[, min_iZ]
> glm_fit <- update(glm_fit, .~.+ X)
> glm_res <- glm_fit$y - glm_fit$fitted.values
> cpvar <- crossprod(glm_res)
> mRES <- MIR_a(glm_res, sW = sW, n = n, cpvar = cpvar,
> S0 = S0, nsim = nsim, stdev = stdev)
> pIZ <- mRES$p.value
> used[min_iZ] <- TRUE
> min_v <- c(min_v, min_iZ)
> if (verbose)
> cat("eV[,", min_iZ, "], I: ", mRES$estimate, " ZI: ",
> mRES$statistic, ", pr(ZI): ", pIZ, "\n", sep = "")
> tres <- rbind(tres, c(min_iZ, mRES$statistic, pIZ))
> }
> sv <- eV[, min_v, drop = FALSE]
> colnames(sv) <- paste("vec", min_v, sep = "")
> colnames(tres) <- c("Eigenvector", "ZI", "pr(ZI)")
> rownames(tres) <- 0:(nrow(tres) - 1)
> res <- list(selection = tres, vectors = sv)
> class(res) <- "ME_res"
> res
> }
>
>
>
--
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no
More information about the R-sig-Geo
mailing list