Thanks Roger for the ideas! Finally the problem was with the input for the
glm.nb function. I copy bellow the modified function ME() for negative
binomial models. It is possible that gives some warnings because with some
eigenvectors the nb model might not converge.

ME.nb <- function (formula, data, weights,listw, 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")
      X <- model.matrix(mt, mf)
      Y <- model.extract(mf, "response")
      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_fit0 <- glm.nb(formula,data=data)
    glm_res <- glm_fit0$y - glm_fit0$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
    tres <- c(NA, mRES$statistic, pIZ)
    if (pIZ > alpha)
        stop("base correlation larger than alpha")
    Cent <- diag(n) - (matrix(1, n, n)/n)
    eV <- eigen(Cent %*% Wmat %*% Cent, EISPACK = TRUE)$vectors
    rm(Cent, Wmat)
    iZ <- numeric(n)
# Select eigenvectors that reduce Moran's Index
    for (i in 1:n) {
        iX <- cbind(X, eV[, i])
        i_glm <- glm.nb(Y~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 <- cbind(X, eV[, min_iZ])
    glm_fit <- glm.nb(Y~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 <- cbind(X, eV[, i])
              i_glm <- glm.nb(Y~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 <- cbind(X, eV[, min_iZ])
        glm_fit <- glm.nb(Y~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,eV=eV)
    class(res) <- "ME_res"
    res
}


2012/1/31 Roger Bivand <Roger.Bivand@nhh.no>

> 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@nhh.no
>
>


-- 

Pablo González Moreno
Department of Integrative Ecology
Estación Biológica de Doñana (EBD-CSIC)
Avda. Américo Vespucio s/n 41092 Sevilla
Tel: 954 466700 (Ext. 1463) Fax: 954 621125

	[[alternative HTML version deleted]]

