[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