Dave Fournier suggests that “NB1” (i.e., a negative binomial model parameterized so that the variance is *proportional* to the mean, rather than the more traditional “NB2” with \( V= \mu(1+\mu/k) \)) is a better fit to the fir data set from Dodd and Silvertown 1999 (reproduced in Bolker 2008 *Ecological Models and Data*). He points out that NB2 is overwhelmingly more popular in the R community (and to be fair the statistics literature in general, although I see that SAS's PROC COUNTREG allows fitting of the NB1 model …), probably because NB2 with a known overdispersion parameter fits into the exponential family that can be handled by the generalized linear model framework.

Load various useful packages/preliminary data manipulation:

```
library(emdbook)
firdat <- na.omit(FirDBHFec[, c("TOTCONES", "DBH", "WAVE_NON")])
firdat <- transform(firdat, TOTCONES = round(TOTCONES))
library(ggplot2)
theme_set(theme_bw())
library(R2admb)
setup_admb()
library(bbmle)
library(reshape2)
```

A basic plot (although in this example we will ignore the distinction between “wave” and “non-wave” populations, at least to begin with)

```
qplot(DBH, TOTCONES, colour = WAVE_NON, shape = WAVE_NON, data = firdat)
```

We'll start by fitting the more traditional NB2 model.

Fitting with ADMB, using the TPL file found in Mollie Brooks' example on the ADMB project web site:

```
compile_admb("fir0")
write_dat("fir0", with(firdat, list(nobs = nrow(firdat), totcones = TOTCONES,
dbh = DBH)))
write_pin("fir0", list(a = 1, b = 1, k = 1)) ## fails with k=5 start
run_admb("fir0")
admb1 <- read_admb("fir0")
```

*FIXME: check order vs. TPL in write_dat !!*

Now try it in R, both with `mle2`

and with `glm.nb`

from the `MASS`

package:

```
r1 <- mle2(TOTCONES ~ dnbinom(mu = a * DBH^b, size = k), start = list(a = 1,
b = 1, k = 1), data = firdat, method = "Nelder-Mead")
r2 <- MASS::glm.nb(TOTCONES ~ log(DBH), data = firdat)
```

It takes slightly more work to extract the `glm.nb`

coefficients in the same format as the other models:

```
r2coef2 <- with(as.list(coef(r2)), c(a = exp(`(Intercept)`), b = `log(DBH)`,
k = r2$theta))
```

One of the nice things about `glm.nb`

(despite its limitations)
is that it makes diagnostic plots easy:

```
par(mfrow = c(2, 2), mar = c(4, 4, 1.5, 1), las = 1, bty = "l", mgp = c(2, 1,
0))
colvec <- adjustcolor(c("blue", "purple"), alpha.f = 0.6)
plot(r2, lwd = 2, col = colvec[firdat$WAVE_NON])
```

The most obvious problem is in the scale-location plot (lower left), where the decreasing line shows that our model variance increases faster than it should with the mean (so that the scaled residual variance decreases with the mean). This is a strong suggestion that the linear mean-variance relationship of NB1 might be more appropriate than the quadratic relationship of NB2 …

The Q-Q plot also has a wonky little wiggle in it. I'm not sure where this comes from – it may be related to the variance scaling problem. I worried that it might be related to the population grouping we were ignoring, so I tried looking at the Q-Q plot split by population (code taking from the examples in `?qqmath`

):

```
firdat$resid <- residuals(r2)
lattice::qqmath(~resid | WAVE_NON, data = firdat, prepanel = prepanel.qqmathline,
panel = function(x, ...) {
panel.qqmathline(x, ...)
panel.qqmath(x, ...)
})
```

Let's try NB1. We want \( V=\tau \mu \) rather than \( V=\mu(1+\mu/k) \), hence to convert between types we should take \( (1+\mu/k)=\tau \), or \( k=\mu/(\tau-1) \).

To do this most conveniently in `mle2`

we define a new density function for NB1:

```
dnbinom1 <- function(x, mu, tau, log = FALSE) {
k <- mu/(tau - 1)
dnbinom(x, mu = mu, size = k, log = log)
}
r3 <- mle2(TOTCONES ~ dnbinom1(mu = a * DBH^b, tau), start = list(a = 1, b = 1,
tau = 5), data = firdat, method = "Nelder-Mead")
```

DF actually used a *robust* NB1 model, where we assume a mixture of two different NB1 models with different \( \tau \) values.

```
compile_admb("firxr")
write_dat("firxr", with(firdat, list(nobs = nrow(firdat), totcones = TOTCONES,
dbh = DBH)))
write_pin("firxr", list(a = 0.3, b = 2, tau = 27, p = 0.05, w = 9))
run_admb("firxr")
admb2 <- read_admb("firxr")
```

*FIXME: make ADMB log-likelihoods get proper dfs*

```
coef(r3)
```

```
## a b tau
## 0.3646 2.2438 27.5132
```

```
AICtab(r1, r3)
```

```
## dAIC df
## r3 0.0 3
## r1 15.2 3
```

```
firdat$predict2 <- with(as.list(coef(r3)), a * firdat$DBH^b)
firdat <- within(firdat, {
predict2 <- with(as.list(coef(r3)), a * DBH^b)
resid2_unsc <- predict2 - TOTCONES
p2_var <- coef(r3)["tau"] * predict2
resid2 <- resid2_unsc/sqrt(p2_var)
})
qplot(predict2, sqrt(abs(resid2)), data = firdat) + geom_smooth(colour = "red",
se = FALSE, method = "loess")
```

```
qplot(sample = resid2, data = firdat)
```

Q-Q plot is still not perfect, but better …

Look at residuals from an `nls`

fit (assuming
constant variance), see how

```
n1 <- nls(TOTCONES ~ a * DBH^b, data = firdat, list(a = 1, b = 1))
ndat <- data.frame(v = residuals(n1)^2, m = fitted(n1))
f1 <- lm(v ~ m - 1, data = ndat)
f2 <- lm(v ~ I(m^2) + offset(m) - 1, data = ndat)
pp <- data.frame(m = seq(min(ndat$m), max(ndat$m), length = 101))
pp <- data.frame(pp, NB1 = predict(f1, newdata = pp), NB2 = predict(f2, newdata = pp))
mpp <- melt(pp, id.var = 1, value.name = "v")
g1 <- qplot(m, v, data = ndat) + geom_line(data = mpp, aes(colour = variable))
g1 + geom_smooth(method = "loess", colour = "purple") + coord_cartesian(ylim = c(0,
16000))
```

(Axis limits exclude one point at approx. m=50, v=30000)

On a log-log plot, it *looks* as though NB2 fits better …

```
g1 + scale_y_log10() + scale_x_log10()
```

How much do the coefficient estimates of all three (NB2, NB1, robust NB1) differ … ?

```
library(coefplot2)
```

```
## Loading required package: coda
```

```
modList <- list(NB2.R = r1, NB2.ADMB = admb1, NB1.R = r3, NB1rob.ADMB = admb2)
cList <- lapply(lapply(modList, coeftab), function(x) x[c("a", "b"), ])
colvec <- c(1, 2, 4, 5)
coefplot2(cList, intercept = TRUE, col = colvec, legend = TRUE)
```

Not much …