Double Exponential Binomial Distribution Family Function
Fits a double exponential binomial distribution by maximum likelihood estimation. The two parameters here are the mean and dispersion parameter.
double.expbinomial(lmean = "logitlink", ldispersion = "logitlink",
idispersion = 0.25, zero = "dispersion")lmean, ldispersion |
Link functions applied to the two parameters, called
mu and theta respectively below.
See |
idispersion |
Initial value for the dispersion parameter. If given, it must be in range, and is recyled to the necessary length. Use this argument if convergence failure occurs. |
zero |
A vector specifying which
linear/additive predictor is to be modelled as intercept-only.
If assigned, the single value can be either |
This distribution provides a way for handling overdispersion in a binary response. The double exponential binomial distribution belongs the family of double exponential distributions proposed by Efron (1986). Below, equation numbers refer to that original article. Briefly, the idea is that an ordinary one-parameter exponential family allows the addition of a second parameter theta which varies the dispersion of the family without changing the mean. The extended family behaves like the original family with sample size changed from n to n*theta. The extended family is an exponential family in mu when n and theta are fixed, and an exponential family in theta when n and mu are fixed. Having 0 < theta < 1 corresponds to overdispersion with respect to the binomial distribution. See Efron (1986) for full details.
This VGAM family function implements an approximation
(2.10) to the exact density (2.4). It replaces the normalizing
constant by unity since the true value nearly equals 1.
The default model fitted is eta1 =logit(mu)
and eta2 = logit(theta).
This restricts both parameters to lie between 0 and 1, although
the dispersion parameter can be modelled over a larger parameter space by
assigning the arguments ldispersion and edispersion.
Approximately, the mean (of Y) is mu. The effective sample size is the dispersion parameter multiplied by the original sample size, i.e., n*theta. This family function uses Fisher scoring, and the two estimates are asymptotically independent because the expected information matrix is diagonal.
An object of class "vglmff" (see vglmff-class).
The object is used by modelling functions such as vglm.
Numerical difficulties can occur; if so, try using idispersion.
This function processes the input in the same way
as binomialff, however multiple responses are
not allowed (binomialff(multiple.responses = FALSE)).
T. W. Yee
Efron, B. (1986). Double exponential families and their use in generalized linear regression. Journal of the American Statistical Association, 81, 709–721.
# This example mimics the example in Efron (1986).
# The results here differ slightly.
# Scale the variables
toxop <- transform(toxop,
phat = positive / ssize,
srainfall = scale(rainfall), # (6.1)
sN = scale(ssize)) # (6.2)
# A fit similar (should be identical) to Section 6 of Efron (1986).
# But does not use poly(), and M = 1.25 here, as in (5.3)
cmlist <- list("(Intercept)" = diag(2),
"I(srainfall)" = rbind(1, 0),
"I(srainfall^2)" = rbind(1, 0),
"I(srainfall^3)" = rbind(1, 0),
"I(sN)" = rbind(0, 1),
"I(sN^2)" = rbind(0, 1))
fit <- vglm(cbind(phat, 1 - phat) * ssize ~
I(srainfall) + I(srainfall^2) + I(srainfall^3) +
I(sN) + I(sN^2),
double.expbinomial(ldisp = extlogitlink(min = 0, max = 1.25),
idisp = 0.2, zero = NULL),
toxop, trace = TRUE, constraints = cmlist)
# Now look at the results
coef(fit, matrix = TRUE)
head(fitted(fit))
summary(fit)
vcov(fit)
sqrt(diag(vcov(fit))) # Standard errors
# Effective sample size (not quite the last column of Table 1)
head(predict(fit))
Dispersion <- extlogitlink(predict(fit)[,2], min = 0, max = 1.25, inverse = TRUE)
c(round(weights(fit, type = "prior") * Dispersion, digits = 1))
# Ordinary logistic regression (gives same results as (6.5))
ofit <- vglm(cbind(phat, 1 - phat) * ssize ~
I(srainfall) + I(srainfall^2) + I(srainfall^3),
binomialff, toxop, trace = TRUE)
# Same as fit but it uses poly(), and can be plotted (cf. Figure 1)
cmlist2 <- list("(Intercept)" = diag(2),
"poly(srainfall, degree = 3)" = rbind(1, 0),
"poly(sN, degree = 2)" = rbind(0, 1))
fit2 <- vglm(cbind(phat, 1 - phat) * ssize ~
poly(srainfall, degree = 3) + poly(sN, degree = 2),
double.expbinomial(ldisp = extlogitlink(min = 0, max = 1.25),
idisp = 0.2, zero = NULL),
toxop, trace = TRUE, constraints = cmlist2)
## Not run: par(mfrow = c(1, 2))
plot(as(fit2, "vgam"), se = TRUE, lcol = "blue", scol = "orange") # Cf. Figure 1
# Cf. Figure 1(a)
par(mfrow = c(1,2))
ooo <- with(toxop, sort.list(rainfall))
with(toxop, plot(rainfall[ooo], fitted(fit2)[ooo], type = "l",
col = "blue", las = 1, ylim = c(0.3, 0.65)))
with(toxop, points(rainfall[ooo], fitted(ofit)[ooo], col = "orange",
type = "b", pch = 19))
# Cf. Figure 1(b)
ooo <- with(toxop, sort.list(ssize))
with(toxop, plot(ssize[ooo], Dispersion[ooo], type = "l", col = "blue",
las = 1, xlim = c(0, 100)))
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.