making inferences about overdispersion in nbinom glmmabmb

Hi Paul,

Standard errors on variance components do not provide a consistent
means to assess either significance or create confidence intervals.
In part because variances are fundamentally bounded, so you really
cannot expect symmetric intervals.  Of course if they are large
enough, that may not be an issue, but as a general principle,
encouraging the use of standard errors assuming a standard normal is
not going to work well.

If you are bound to a frequentist world, a likelihood ratio test is
probably a good choice with a free and constrained estimate.  To
generate confidence intervals, a good approach can be used based on
profile likelihoods.  I am sure these are presented elsewhere, but
Terry Therneau has a good presentation of the technique and theory in
his text on survival models.

Below is an alternative Bayesian approach.  The parameterization of
course is slightly different as this uses MCMCglmm and a non zero
residual variance rather than negative binomial model, but the
posterior allows valid examination and computation of confidence
intervals and presentation of the distribution.  Note for larger
values, it is not an issue, but in the last example (admittedly chosen
as a degenerate case) normal based confidence intervals would perform
quite poorly.

Best,

Josh

require(MASS)
require(MCMCglmm)

set.seed(10)
d <- data.frame(x = rnorm(1000))

d <- within(d, {
  y1 = rnegbin(1000, 2 + .5 * x, 1)
  y2 = rnegbin(1000, 2 + .5 * x, 4)
  y3 = rnegbin(1000, 2 + .5 * x, 16)
})

m <- list(
  t1 = MCMCglmm(y1 ~ x, data = d, family = "poisson", nitt=1e5, verbose=FALSE),
  t4 = MCMCglmm(y2 ~ x, data = d, family = "poisson", nitt=1e5, verbose=FALSE),
  t16 = MCMCglmm(y3 ~ x, data = d, family = "poisson", nitt=1e5, verbose=FALSE))

plot(m$t1$VCV)
plot(m$t4$VCV)
plot(m$t16$VCV)

HPDinterval(m$t1$VCV)
HPDinterval(m$t4$VCV)
HPDinterval(m$t16$VCV)

p <- function(x) {
  hist(x, freq=FALSE)
  lines(seq(min(x), max(x), length.out = 1000),
  dnorm(seq(min(x), max(x), length.out = 1000), mean= mean(x), sd = sd(x)),
  type = "l", col = "blue")
}

p(m$t1$VCV[, 1])
p(m$t4$VCV[, 1])
p(m$t16$VCV[, 1])

On Wed, Jun 19, 2013 at 5:17 AM, Paul Johnson

making inferences about overdispersion in nbinom glmmabmb

Thread (4 messages)