3
%- Also NEED an '\alias' for EACH other topic documented here.
4
\title{ Negative binomial canonical link function }
6
Computes the negative binomial canonical link transformation,
7
including its inverse and the first two derivatives.
11
nbcanlink(theta, earg = list(), inverse = FALSE, deriv = 0,
12
short = TRUE, tag = FALSE)
14
%- maybe also 'usage' for other objects documented here.
18
Typically the mean of a negative binomial (NB) distribution.
19
See below for further details.
25
Extra argument for passing in additional information.
26
Here, a \code{size} component contains the \eqn{k} matrix which
27
must be of a conformable dimension as \code{theta}.
28
Also, if \code{deriv > 0} then a \code{wrt.eta} component
29
which is either 1 or 2 (1 for with respect to the first
30
linear predictor, and 2 for with respect to the second
31
linear predictor (a function of \eqn{k})).
35
\item{inverse}{ Logical. If \code{TRUE} the inverse function is computed. }
36
\item{deriv}{ Order of the derivative. Integer with value 0, 1 or 2. }
38
Used for labelling the \code{blurb} slot of a
39
\code{\link{vglmff-class}} object.
43
Used for labelling the linear/additive predictor in the
44
\code{initialize} slot of a \code{\link{vglmff-class}} object.
45
Contains a little more information if \code{TRUE}.
51
The negative binomial (NB) canonical link is
52
\eqn{\log(\theta/ (\theta + k))}{log(theta/(theta + k))}
53
where \eqn{\theta}{theta} is the mean of a NB
54
distribution. The canonical link is used for theoretically
55
relating the NB to GLM class.
58
This link function was specifically written for
59
\code{\link{negbinomial}} and
60
\code{\link{negbinomial.size}},
61
and should not be used elsewhere
62
(these \pkg{VGAM} family functions have code that
63
specifically handles \code{nbcanlink()}.)
68
For \code{deriv = 0}, the above equation
69
when \code{inverse = FALSE}, and if \code{inverse = TRUE} then
70
\code{kmatrix / expm1(-theta)}.
71
For \code{deriv = 1}, then the function returns
72
\emph{d} \code{theta} / \emph{d} \code{eta} as a function of \code{theta}
73
if \code{inverse = FALSE},
74
else if \code{inverse = TRUE} then it returns the reciprocal.
81
Two-parameter reduced-rank vector generalized linear models.
82
\emph{In preparation}.
86
\emph{Negative Binomial Regression},
88
Cambridge: Cambridge University Press.
92
\author{ Thomas W. Yee }
95
This function currently does not work very well with \code{\link{negbinomial}}!
96
The NB-C model is sensitive to the initial values and may converge to a local solution.
97
Pages 210 and 309 of Hilbe (2011) notes convergence difficulties (of
98
Newton-Raphson type algorithms), and this applies here.
99
This function should work okay with \code{\link{negbinomial.size}}.
100
Currently trying something like \code{imethod = 3} or \code{imu},
101
\code{stepsize = 0.5}, \code{maxit = 100}, \code{zero = -2} should help;
102
see the example below.
109
While theoretically nice, this function is not recommended
110
in general since its value is always negative (linear predictors
111
ought to be unbounded in general). A \code{\link{loge}}
112
link for argument \code{lmu} is recommended instead.
115
Numerical instability may occur when \code{theta} is close to 0 or 1.
116
For the \code{earg} argument,
117
values of \code{theta} which are less than or equal to 0 can be
118
replaced by the \code{bvalue} component of the list \code{earg}
119
before computing the link function value.
120
The component name \code{bvalue} stands for ``boundary value''.
121
See \code{\link{Links}} for general information about \code{earg}.
128
\code{\link{negbinomial}},
129
\code{\link{negbinomial.size}}.
134
nbcanlink("mu", short = FALSE)
136
mymu = 1:10 # Test some basic operations:
137
kmatrix = matrix(runif(length(mymu)), length(mymu), 1)
138
eta1 = nbcanlink(mymu, earg = list(size = kmatrix))
139
ans2 = nbcanlink(eta1, earg = list(size = kmatrix), inverse = TRUE)
140
max(abs(ans2 - mymu)) # Should be 0
142
\dontrun{ mymu = c(seq(0.5, 10, length = 101))
143
kmatrix = matrix(10, length(mymu), 1)
144
plot(nbcanlink(mymu, earg = list(size = kmatrix)) ~ mymu, las = 1,
145
type = "l", col = "blue", lwd = 1.5, xlab = expression({mu})) }
147
# Estimate the parameters from some simulated data (see Warning section)
149
ndata <- data.frame(x2 = runif(nn <- 1000 ))
150
size1 = exp(1); size2 = exp(2)
151
ndata <- transform(ndata, eta1 = -1 - 2 * x2, # eta1 < 0
154
ndata <- transform(ndata,
155
mu1 = nbcanlink(eta1, earg = list(size = size1), inv = TRUE),
156
mu2 = nbcanlink(eta1, earg = list(size = size2), inv = TRUE))
157
ndata <- transform(ndata, y1 = rnbinom(nn, mu = mu1, size = size1),
158
y2 = rnbinom(nn, mu = mu2, size = size2))
162
fit <- vglm(cbind(y1, y2) ~ x2, negbinomial("nbcanlink", imethod = 3),
163
stepsize = 0.5, ndata, # Deliberately slow the convergence rate
164
maxit = 100, trace = TRUE) # Warning: may converge to a local soln
165
coef(fit, matrix = TRUE)
171
% abline(h = 0, col = "lightgray", lty = "dashed", lwd = 2.0)
172
% The variance-covariance matrix may be wrong when the
173
% canonical link is used.
174
% vcov(fit) # May be wrong