Using Spatial Eigenvector Mapping for Negative Binomial GLM

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.

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
}