Spatially correlated negative binomial models

Define a spatially correlated negative binomial model. For the moment, leave the spatial correlation model
relatively unspecified as $\mathbf V = \sigma^2 \mathbf{C}(r)$; $\mathbf{C}$ will probably end up being something relatively simple like an exponential model ($C_{ij} = e^{-r_{ij}/s}$, where $s$ is a correlation scale. The statistical model for the observed number of stems of species $k$ in location $i$ at time $t+1$ is:

(1)
\begin{split} L_{ik}(t+1) & = \beta_0+\sum_j \beta_{kj}N_{ij}(t) \\ M_{ik}(t+1) & \sim \mbox{MVN} ({\mathbf L}(t+1), \sigma^2 \mathbf C) \\ N_{ik}(t+1) & \sim \mbox{Poisson}(N_{ik}(t) \exp(M_{ik}(t+1))) \end{split}

In words: the expected densities follow a Ricker-like equation, with intra- and interspecific interaction coefficients $\beta_{kj}$ (not sure I have this the right/conventional way around). $L$ is the linear predictor, the expected log population growth rate. $\mathbf M$ is a "latent variable" that adds normally distributed variability to these expected densities, but also adds (through the correlation matrix $\mathbf C$) a potential for spatial correlation in that variation. The observed stem numbers are then a Poisson draw with a log-mean equal to the latent variable. (The reason to go through all this nonsense with the latent variable is that the multivariate normal distribution is by far the easiest way to define spatial correlation structures.) One weakness of this model (which I won't try to deal with at the moment) is that since there is no immigration term, local extinction is permanent —- if $N_{ik}(t)=0$ then $N_{ik}(t+1)=0$. There is also no observational error.

Simulating from this system:

> set.seed(1001)  ## set random number seed
> locs <- expand.grid(1:10,1:50) ## 10 x 50 grid of locations
> ## year-1 values for two species
> ## (seed with negative binomial even though we 
> ##  will be using LN-Poisson later on)
> init1 <- rnbinom(500,mu=2,size=1)
> init2 <- rnbinom(500,mu=1,size=2)
> N0 <- cbind(init1,init2)
> beta0 <- 1  ## dens-ind. growth rate
> ## interaction matrix -- inter > intra 
> betamat <- matrix(c(-0.2,-0.5,-0.5,-0.2),ncol=2)
> linpred <- beta0+ N0 %*% betamat ## linear predictor
> require(MASS)
> dmat <- as.matrix(dist(locs))
> Cor <- exp(-dmat/2)
> sigmasq <- c(0.5,0.25)
> ## assume same correlation structure for now
> M1 <- mvrnorm(1,mu=linpred[,1],Sigma=sigmasq[1]*Cor)
> M2 <- mvrnorm(1,mu=linpred[,2],Sigma=sigmasq[2]*Cor)
> ## could construct one giant block-diagonal var-cov
> ## matrix and pick a single MVN deviate, but probably
> ## not worth the trouble
> N1.1 <- rpois(500,N0[,1]*exp(M1))
> N1.2 <- rpois(500,N0[,2]*exp(M2))

Fitting without regard to spatial structure (negative-binomial rather than LN-Poisson):
> library(bbmle)
> mod1 <- mle2(N1.1 ~ dnbinom(mu = N0[, 1] * exp(beta0 + beta11 * 
+     N0[, 1] + beta12 * N0[, 2]), size = exp(logk)), start = list(beta0 = 0, 
+     beta11 = 0, beta12 = 0, logk = 0.5))

Or:
> library(MASS)
> mod2 <- glm.nb(N1.1 ~ N0[, 1] + N0[, 2] + offset(log(N0[, 1] + 
+     0.001)))

Have to cheat a bit on the offset here to avoid $\log(0)$. Could tweak this further (probably) and fit the whole thing in one shot, but leave that for later. Compare fits:
> (rbind(coef(mod1), c(coef(mod2), log(mod2$theta))))
        beta0     beta11     beta12      logk
[1,] 1.124450 -0.2033449 -0.5252317 0.2813700
[2,] 1.120297 -0.2026524 -0.5251838 0.2817174

That's good. Let's plot a map of the residuals and see if we can see anything.
> rmat <- t(matrix(residuals(mod2), nrow = 10))
> image(rmat, col = gray((1:50)/50))

spatnbinom-005.png
It does look as though there's some autocorrelation here (there should be!). Next step is to try Moran's $I$; this is not entirely trivial (hint: install the spdep package and figure out how to use moran.test …) WinBUGS? Or is there any easier way to fit this?
## construct linear predictor
for (k in 1:nspecies) {
  for (i in 1:nlocs) {
     linpred[i,k] <- beta0+inprod(betamat[k,],N1[i,])
  }
}
## construct var-cov matrix
## do we have to construct a PRECISION matrix in WinBUGS?
##  ugh ...
for (j in 1:nlocs) {
  for (i in 1:nlocs) {
    for (k in 1:nspecies) {
        vmat[i,j,k] <- sigmasq ....
    }
  }
}
##
 do something to get the MVN variates ...
for 
for (k in 1:nspecies) {
  for (i in 1:nlocs) {
     expN[i,k] <- N1[i,k]*exp(latent[i,k])
     N2[i,k] ~ dpois(expN[i,k]
  }
}
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License