fir analysis

Introduction

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.

Preliminaries

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)

plot of chunk plot1

NB2

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])

plot of chunk unnamed-chunk-4

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, ...)
    })

plot of chunk unnamed-chunk-5

NB1

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")

plot of chunk unnamed-chunk-9

qplot(sample = resid2, data = firdat)

plot of chunk unnamed-chunk-9

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

Scaling plots

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))

plot of chunk unnamed-chunk-10

(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()

plot of chunk unnamed-chunk-11

Predictions

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)

plot of chunk coefplot

Not much …