One of the relaxations of TRIM with respect to the Poission
regression assumptions is that variance does not need to equal variance,
but instead is allowed to be proportionally larger, i.e. \(\operatorname{var}[x] = \sigma^2
\operatorname{E}[x]\), where \(\sigma^2\) is an overdispersion parameter,
which is a scalar parameter estimated by TRIM.
This alternative constraint on variance is easily expressed using a
negative binomial distribution \(\operatorname{NB}(r,p)\) where \(r\) is a parameter and \(p\) is a . In ecological applications an
alternative parameterization is often used, using \(r\) and mean \(\mu\). In this case, the first two moments
are given by \[ \operatorname{E}[x] = \mu
\] and \[ \operatorname{var}[x] = \mu
+ \frac{\mu^2}{r} \] Using the overdispersion constraint \(\operatorname{var}[x] \equiv \sigma^2
\operatorname{E}[x]\) we get \[ \mu +
\frac{\mu^2}{r} \equiv \sigma^2 \mu \] \[ \Rightarrow 1 + \frac{\mu}{r} = \sigma^2
\] \[ \Rightarrow r =
\frac{\mu}{\sigma^2-1} \] Here is an example showing the CDF’s
for \(\sigma^2=1\) (Poisson) and \(\sigma^2=2,5,10\), using \(\mu=\lambda=100\).
mu <- 100L
sig2 <- c(1, 2, 5, 10)
n <- length(sig2)
colors <- c("#E41A1C","#377EB8","#4DAF4A","#984EA3") # ColorBrewer "Set1" colors 1..4
x <- seq.int(0, 2*mu, by=2)
plot(NULL, NULL, type="n", xlim=range(x), ylim=c(0,1),
xlab="Value", ylab="Cumulative distribution function", las=1)
for (i in 1:n) {
y <- if (sig2[i]==1) ppois(x, mu)
else pnbinom(x, mu=mu, size=mu/(sig2[i]-1))
points(x, y, col=colors[i], pch=16, cex=0.5)
}
# mark expected value
abline(v=mu, col=mgray, lty="dashed")
text(mu, 0.1, expression(lambda), pos=4, col=mgray)
# add legend
leg_msg <- sprintf("Overdispersion = %d", sig2)
leg_msg[1] <- "Poisson"
legend("topleft", legend=leg_msg, col=colors[1:n], pch=16)
As with the Poisson distribution, the Negative Binomial distribution
can well be approximated by a Gamma distribution. Again, we use the
methods of equal moments to find the Gamma parameters \(k,\theta\) corresponding to the relaxed
assumptions \(\operatorname{var}[x] = \sigma^2
\operatorname{E}[x] = \sigma^2\mu = \sigma^2\lambda\): \[ \operatorname{E}[x] = k\theta = \mu \Rightarrow
\theta = \frac{\mu}{k} \] \[\operatorname{var}[x] = \sigma^2\mu \Rightarrow k
\theta^2 = \sigma^2\mu\] \[
\Rightarrow k \left(\frac{\mu}{k}\right)^2 = \sigma^2 \mu \]
\[ \Rightarrow k = \frac{\mu}{\sigma^2}
\] Going back to our expression for \(\theta\): \[
\theta = \frac{\mu}{k} \Rightarrow \theta=\sigma^2 \]
Here is the previous example of four different overdispersion
parameters, for both the Poisson / Negative Binomial distibutions
enhanced by plotting the corresponding Gamma distributions.
plot(NULL, NULL, type="n", xlim=range(x), ylim=c(0,1),
xlab="", ylab="Cumulative distribution function", las=1)
# Plot discrete distributions: Poisson or Negative Binonmial
for (i in 1:n) {
y <- if (sig2[i]==1) ppois(x, mu)
else pnbinom(x, mu=mu, size=mu/(sig2[i]-1))
points(x, y, col=colors[i], pch=16, cex=0.5)
}
# Plot continuous distributions: Gamma
for (i in 1:n) {
y <- pgamma(x, shape=mu/sig2[i], scale=sig2[i])
lines(x, y, col=colors[i])
}
# mark expected value
abline(v=mu, col=mgray, lty="dashed")
text(mu, 0.1, expression(lambda), pos=4)
# add legend
leg_msg <- sprintf("Overdispersion = %d", sig2)
leg_msg[1] <- "Poisson"
legend("topleft", legend=leg_msg, col=colors[1:n], pch=16)
Multipliers
Multipliers relating standard errors with confidence interval bounds
can be computed using the same approach as demonstrated above for the
Poisson distribution. Here is an example for \(\sigma^2 = 1 \ldots 10\)
lambda <- exp(seq(log(2.0), log(1000.0), len=100)) # uniform spacing in log-space
alpha <- 0.05 # 95% CI
sig2 <- c(1, 2, 5, 10)
n <- length(sig2)
colors <- c("#E41A1C","#377EB8","#4DAF4A","#984EA3") # ColorBrewer "Set1" colors 1..4
plot(NULL, NULL, xlim=range(lambda), ylim=c(0,3), log="x",
xlab="Expected value", ylab="Multiplier", las=1)
for (i in 1:n) {
qhi <- qgamma(p=1-alpha/2, shape=lambda/sig2[i], scale=sig2[i])
qlo <- qgamma(p=alpha/2, shape=lambda/sig2[i], scale=sig2[i])
sd <- sqrt(sig2[i] * lambda)
umul <- (qhi-lambda) / sd
lmul <- (lambda-qlo) / sd
lines(lambda, umul, col=colors[i], lwd=2)
lines(lambda, lmul, col=colors[i], lwd=2)
# idx <- which(lambda>=10)[1]
# text(lambda[idx],umul[idx],"M_hi", pos=3, col=red)
# text(lambda[idx],lmul[idx],"M_lo", pos=1, col=red)\
}
m0 <- qnorm(0.975) # Multiplier for normal distibutions, i.e. 1.96
abline(h=m0, col=mgray, lty="dashed", lwd=2)
text(2.5,m0,sprintf("%.2f",m0), col=mgray, pos=1)
leg_msg <- sprintf("Overdispersion = %d", sig2)
leg_msg[1] <- "Poisson"
legend("bottomright", legend=leg_msg, col=colors[1:n], pch=16)
Note that while upper-bound multipliers increase with overdispersion,
lower-bound multipliers decrease. This is easily understood from the
skewness of these distributions and the constraints that the lower bound
is constrained to be \(>0\). The
following plot illustrates the principle for \(\lambda=20\): and \(\sigma^2=8\)
lambda <- 20
sig2 <- 8
alpha <- 0.05 # 95% CI
n <- length(sig2)
# Discrete: Poisson & Neg. Binomial
x <- 0 : (3*lambda)
y <- dnbinom(x, mu=lambda, size=lambda/(sig2-1))
plot(x,y, xlab="Value", ylab="Frequency", type='s', col=red, las=1)
points(x, y, col=red, pch=16, cex=0.5)
# Expected value as solid vertical line
abline(v=lambda, col=mgray, lwd=2)
# Standard errors as dashed lines
sd <- sqrt(sig2 * lambda)
abline(v=c(lambda-sd,lambda+sd), lty="dashed", col=blue, lwd=2)
# plot CI intervals as dotted lines
lo <- qgamma(p= alpha/2, shape=lambda/sig2, scale=sig2)
hi <- qgamma(p=1-alpha/2, shape=lambda/sig2, scale=sig2)
abline(v=c(lo,hi), lty="dotted", col=green, lwd=2)
legend("topright", c("mean","std.err","C.I."), lty=c("solid","dashed","dotted"),
lwd=2, col=c(mgray,blue,green))
Things change when \(\sigma^2 >
\lambda\), in this case the standard errors are larger than the
mean, such that the multipliers must be \(<1\) to ensure the confidence interval
be positive. Here is an example for \(\lambda=10\) and \(\sigma^2=12\)
lambda <- 10
sig2 <- 12
alpha <- 0.05 # 95% CI
n <- length(sig2)
# Discrete: Poisson & Neg. Binomial
x <- 0 : (3*lambda)
y <- dnbinom(x, mu=lambda, size=lambda/(sig2-1))
plot(x,y, xlab="Value", ylab="Frequency", type='s', col=red, las=1)
points(x, y, col=red, pch=16, cex=0.5)
# Expected value as solid vertical line
abline(v=lambda, col=mgray, lwd=2)
# Standard errors as dashed lines
sd <- sqrt(sig2 * lambda)
abline(v=c(lambda-sd,lambda+sd), lty="dashed", col=blue, lwd=2)
# plot CI intervals as dotted lines
lo <- qgamma(p= alpha/2, shape=lambda/sig2, scale=sig2)
hi <- qgamma(p=1-alpha/2, shape=lambda/sig2, scale=sig2)
abline(v=c(lo,hi), lty="dotted", col=green, lwd=2)
legend("topright", c("mean","std.err","C.I."), lty=c("solid","dashed","dotted"),
lwd=2, col=c(mgray,blue,green))
It is recommended that these cases (i.e. \(\sigma^2 > \lambda\)) be avoided when
fitting rtrim models.
The folowing diagram shows the multipliers for a range of
overdispersion parameters, indicating the domain of applicability
(i.e. \(\sigma^2 < \lambda\))
lambda <- exp(seq(log(2.0), log(1000.0), len=100)) # uniform spacing in log-space
alpha <- 0.05 # 95% CI
sig2 <- c(1, 3, 10, 30, 100)
n <- length(sig2)
colors <- c("#E41A1C","#377EB8","#4DAF4A","#984EA3","#FF7F00") # ColorBrewer "Set1" colors 1..5
plot(NULL, NULL, xlim=range(lambda), ylim=c(0,3), log="x",
xlab="Expected value", ylab="Multiplier")
for (i in 1:n) {
qhi <- qgamma(p=1-alpha/2, shape=lambda/sig2[i], scale=sig2[i])
qlo <- qgamma(p=alpha/2, shape=lambda/sig2[i], scale=sig2[i])
sd <- sqrt(sig2[i] * lambda)
umul <- (qhi-lambda) / sd
lmul <- (lambda-qlo) / sd
lines(lambda, umul, col=colors[i], lwd=1, lty="dashed")
lines(lambda, lmul, col=colors[i], lwd=1, lty="dashed")
# plot feasible range
ok <- lambda > sig2[i]
lines(lambda[ok], umul[ok], col=colors[i], lwd=2)
lines(lambda[ok], lmul[ok], col=colors[i], lwd=2)
}
abline(h=qnorm(0.975), col=mgray, lty="dashed", lwd=2)
leg_msg <- sprintf("Overdispersion = %d", sig2)
leg_msg[1] <- "Poisson"
legend("bottomright", legend=leg_msg, col=colors[1:n], lty="solid")
Note that because of our feasibility constraints all mulitplier–expected
value relations look similar now, expcept for the obvious shift to the
right.
Again, these general-case multipliers (i.e. both Poisson and
Overdispersion assumptions) can be applied to the Skylark example to
create and plot confidence intervals.
# Run TRIM with overdispersion
m1 <- trim(count ~ site + year, data=skylark2, model=3, overdisp=TRUE)
tt1 <- totals(m1, level=0.95)
plot(tt1, main=sprintf("Skylark; overdispersion=%.2f", overdispersion(m1)))
In this case, the effect of overdispersion is very limited. For an
example where overdispersion is large, see the vignettes Taming overdispersion and [rtrim
2.0 extensions](rtrim_2_extensions.html
The problem with the Poisson distribution in this context is it’s
discrete nature: the CDF is discontinuous, so a given confidence level
of, say, 95% is not uniquely related to corresponding standard errors of
counts.
This is perfectly acceptable.
showing again the acceptable degree of approximation.