gam modelling function is designed to be able to use
negbin family (a modification of MASS library
by Venables and Ripley), or the
nb function designed for integrated estimation of
theta. θ is the parameter such that var(y) = μ + μ^2/θ, where μ = E(y).
Two approaches to estimating
theta are available (with
negbin then if ‘performance iteration’ is used for smoothing parameter estimation
gam), then smoothing parameters are chosen by GCV and
theta is chosen in order to ensure that the Pearson estimate of the scale
parameter is as close as possible to 1, the value that the scale parameter should have.
If ‘outer iteration’ is used for smoothing parameter selection with the
nb family then
theta is estimated alongside the smoothing parameters by ML or REML.
To use the first option, set the
optimizer argument of
"perf" (it can sometimes fail to converge).
negbin(theta = stop("'theta' must be specified"), link = "log") nb(theta = NULL, link = "log")
Either i) a single value known value of theta or ii) two values of theta specifying the
endpoints of an interval over which to search for theta (this is an option only for
The link function: one of
nb allows estimation of the
theta parameter alongside the model smoothing parameters, but is only usable with
negbin, if a single value of
theta is supplied then it is always taken as the known fixed value and this is useable with
theta is two
theta>theta) then they are taken as specifying the range of values over which to search for
the optimal theta. This option is deprecated and should only be used with performance iteration estimation (see
optimizer), in which case the method
of estimation is to choose theta so that the GCV (Pearson) estimate
of the scale parameter is one (since the scale parameter
is one for the negative binomial). In this case theta estimation is nested within the IRLS loop
used for GAM fitting. After each call to fit an iteratively weighted additive model to the IRLS pseudodata,
the theta estimate is updated. This is done by conditioning on all components of the current GCV/Pearson
estimator of the scale parameter except theta and then searching for the
theta which equates this conditional estimator to one. The search is
a simple bisection search after an initial crude line search to bracket one. The search will
terminate at the upper boundary of the search region is a Poisson fit would have yielded an estimated
scale parameter <1.
negbin an object inheriting from class
family, with additional elements
the function giving the first derivative of the variance function w.r.t.
the function giving the second derivative of the variance function w.r.t.
A function for retrieving the value(s) of theta. This also useful for retriving the
nb an object inheriting from class
gamm does not support
The negative binomial functions from the MASS library are no longer supported.
Simon N. Wood firstname.lastname@example.org
modified from Venables and Ripley's
Venables, B. and B.R. Ripley (2002) Modern Applied Statistics in S, Springer.
Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter and model selection for general smooth models. Journal of the American Statistical Association 111, 1548-1575 doi: 10.1080/01621459.2016.1180986
library(mgcv) set.seed(3) n<-400 dat <- gamSim(1,n=n) g <- exp(dat$f/5) ## negative binomial data... dat$y <- rnbinom(g,size=3,mu=g) ## known theta fit ... b0 <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=negbin(3),data=dat) plot(b0,pages=1) print(b0) ## same with theta estimation... b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=nb(),data=dat) plot(b,pages=1) print(b) b$family$getTheta(TRUE) ## extract final theta estimate ## another example... set.seed(1) f <- dat$f f <- f - min(f)+5;g <- f^2/10 dat$y <- rnbinom(g,size=3,mu=g) b2 <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=nb(link="sqrt"), data=dat,method="REML") plot(b2,pages=1) print(b2) rm(dat)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.