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.biv...@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.biv...@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]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
