addreg (version 3.0)

addreg: Additive Regression for Discrete Data

Description

addreg fits additive (identity-link) Poisson, negative binomial and binomial regression models using a stable combinatorial EM algorithm.

Usage

addreg(formula, mono = NULL, family, data, standard, subset, na.action,
       start = NULL, offset, control = list(...), model = TRUE, 
       method = c("cem", "em"),
       accelerate = c("em", "squarem", "pem", "qn"), 
       control.method = list(), warn = TRUE, ...)

Arguments

formula

an object of class "formula" (or one that can be coerced into that class): a symbolic description of the model to be fitted. The details of model specification are given under "Details". Note that the model must contain an intercept, and 2nd-order terms (such as interactions) or above are currently not supported --- see "Note".

mono

a vector indicating which terms in formula should be restricted to have a monotonically non-decreasing relationship with the outcome. May be specified as names or indices of the terms.

family

a description of the error distribution to be used in the model. This can be a character string naming a family function, a family function or the result of a call to a family function (see family for details of family functions), but here it is restricted to be poisson, negbin1 or binomial family with identity link.

data

an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which addreg is called.

standard

a numeric vector of length equal to the number of cases, where each element is a positive constant that (multiplicatively) standardises the fitted value of the corresponding element of the response vector. Ignored for binomial family (two-column specification of response should be used instead).

subset

an optional vector specifying a subset of observations to be used in the fitting process.

na.action

a function which indicates what should happen when the data contain NAs. The default is set be the na.action setting of options, and is na.fail if that is unset. The `factory-fresh' default is na.omit. Another possible value is NULL, no action. Value na.exclude can be useful.

start

starting values for the parameters in the linear predictor, also with the starting value for the scale as the last element when family = negbin1.

offset

this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a non-negative numeric vector of length equal to the number of cases. One or more offset terms can be included in the formula instead or as well, and if more than one is specified their sum is used. See model.offset.

Ignored for binomial family; not yet implemented for negative binomial models.

control

list of parameters for controlling the fitting process, passed to addreg.control.

model

a logical value indicating whether the model frame (and, for binomial models, the equivalent Poisson model) should be included as a component of the returned value.

method

a character string that determines which algorithm to use to find the MLE: "cem" for the combinatorial EM algorithm, which cycles through a sequence of constrained parameter spaces, or "em" for a single EM algorithm based on an overparameterised model.

accelerate

a character string that determines the acceleration algorithm to be used, (partially) matching one of "em" (no acceleration --- the default), "squarem", "pem" or "qn". See turboem for further details. Note that "decme" is not permitted.

control.method

a list of control parameters for the acceleration algorithm, which are passed to the control.method argument of turboem.

If any items are not specified, the defaults are used.

warn

a logical indicating whether or not warnings should be provided for non-convergence or boundary values.

…

arguments to be used to form the default control argument if it is not supplied directly.

Value

addreg returns an object of class "addreg", which inherits from classes "glm" and "lm". The function summary.addreg can be used to obtain or print a summary of the results.

The generic accessor functions coefficients, fitted.values and residuals can be used to extract various useful features of the value returned by addreg. Note that effects will not work.

An object of class "addreg" is a list containing the same components as an object of class "glm" (see the "Value" section of glm), but without contrasts, qr, R or effects components. It also includes:

loglik

the maximised log-likelihood.

aic.c

a small-sample corrected version of Akaike's An Information Criterion (Hurvich, Simonoff and Tsai, 1998). This is used by addreg.smooth to choose the optimal number of knots for smooth terms.

xminmax

the minimum and maximum observed values for each of the continuous covariates, to help define the covariate space of the model.

As well as, for Poisson and negative binomial models:

nn.coefficients

estimated coefficients associated with the non-negative parameterisation corresponding to the MLE.

nn.x

non-negative model matrix associated with nn.coefficients.

standard

the standard argument.

Or, for binomial models:

model.addpois

if requested, the addreg object for the associated identity-link Poisson model.

The scale component of the result is fixed at 1 for Poisson and binomial models, and is the constant overdispersion parameter for negative binomial models (that is, scale = 1+\phi) where Var(\mu) = (1+\phi)\mu).

Details

addreg fits a generalised linear model (GLM) with a Poisson or binomial error distribution and identity link function, as well as additive NegBin I models (which are not GLMs). Predictors are assumed to be continuous, unless they are of class factor, or are character or logical (in which case they are converted to factors). Specifying a predictor as monotonic using the mono argument means that for continuous terms, the associated coefficient will be restricted to be non-negative, and for categorical terms, the coefficients will be non-decreasing in the order of the factor levels. This allows semi-parametric monotonic regression functions, in the form of unsmoothed step-functions. For smooth regression functions see addreg.smooth.

As well as allowing monotonicity constraints, the function is useful when a standard GLM routine, such as glm, fails to converge with an identity-link Poisson or binomial model. If glm does achieve successful convergence, and addreg converges to an interior point, then the two results will be identical. However, glm may still experience convergence problems even when addreg converges to an interior point. Note that if addreg converges to a boundary point, then it may differ slightly from glm even if glm successfully converges, because of differences in the definition of the parameter space. addreg produces valid fitted values for covariate values within the Cartesian product of the observed range of covariate values, whereas glm produces valid fitted values just for the observed covariate combinations (assuming it successfully converges). This issue is only relevant when addreg converges to a boundary point.

The computational method is a combinatorial EM algorithm (Marschner, 2014), which accommodates the parameter contraints in the model and is more stable than iteratively reweighted least squares. A collection of restricted parameter spaces is defined which covers the full parameter space, and the EM algorithm is applied within each restricted parameter space in order to find a collection of restricted maxima of the log-likelihood function, from which can be obtained the global maximum over the full parameter space. See Marschner (2010), Donoghoe and Marschner (2014) and Donoghoe and Marschner (2016) for further details.

Acceleration of the EM algorithm can be achieved through the methods of the turboEM package, specified through the accelerate argument. However, note that these methods do not have the guaranteed convergence of the standard EM algorithm, particularly when the MLE is on the boundary of its (possibly constrained) parameter space.

References

Donoghoe, M. W. and I. C. Marschner (2014). Stable computational methods for additive binomial models with application to adjusted risk differences. Computational Statistics and Data Analysis 80: 184--196.

Donoghoe, M. W. and I. C. Marschner (2016). Estimation of adjusted rate differences using additive negative binomial regression. Statistics in Medicine 35(18): 3166--3178.

Hurvich, C. M., J. S. Simonoff and C.-L. Tsai (1998). Smoothing parameter selection in nonparametric regression using an improved Akaike information criterion. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 60(2): 271--293.

Marschner, I. C. (2010). Stable computation of maximum likelihood estimates in identity link Poisson regression. Journal of Computational and Graphical Statistics 19(3): 666--683.

Marschner, I. C. (2014). Combinatorial EM algorithms. Statistics and Computing 24(6): 921--940.

Examples

# NOT RUN {
require(glm2)
data(crabs)

#============================================================================
#  identity-link Poisson model with periodic non-convergence when glm is used
#============================================================================

crabs.boot <- crabs[crabs$Rep1,-c(5:6)]
crabs.boot$width.shifted <- crabs.boot$Width - min(crabs$Width)

fit.glm <- glm(Satellites ~ width.shifted + factor(Dark) + factor(GoodSpine),
  family = poisson(identity), data = crabs.boot, start = rep(1,4), 
  control = glm.control(trace = TRUE))
  
fit.addreg <- addreg(formula(fit.glm), family = poisson, data = crabs.boot,
  trace = 1)

# Speed up convergence by using single EM algorithm
fit.addreg.em <- update(fit.addreg, method = "em")

# Speed up convergence by using acceleration methods
fit.addreg.acc <- update(fit.addreg, accelerate = "squarem")
fit.addreg.em.acc <- update(fit.addreg.em, accelerate = "squarem")
  
# Usual S3 methods work on addreg objects
summary(fit.addreg)
vcov(fit.addreg)
confint(fit.addreg)
summary(predict(fit.addreg), type = "response")

fit.addreg2 <- addreg(update(formula(fit.glm), ~ . - factor(GoodSpine)),
  family = poisson, data = crabs.boot, trace = 1)
anova(fit.addreg2, fit.addreg, test = "LRT")

# Account for overdispersion (use start to speed it up a little)
fit.addreg.od <- addreg(Satellites ~ factor(Dark) + factor(GoodSpine),
  family = negbin1, data = crabs.boot, trace = 1,
  start = c(4.3423675,-2.4059273,-0.4531984,5.969648))
summary(fit.addreg.od)
# }