Generalized Linear Models with Extra Parameters

Estimation of generalized linear models with extra parameters, e.g., parametric links, or families with additional parameters (such as negative binomial).

Usage

glmx(formula, data, subset, na.action, weights, offset,
  family = negative.binomial, xlink = "log", control = glmx.control(...),
  model = TRUE, y = TRUE, x = FALSE, ...)

glmx.fit(x, y, weights = NULL, offset = NULL,
  family = negative.binomial, xlink = "log", control = glmx.control())

Arguments

formula: symbolic description of the model.
data, subset, na.action: arguments controlling formula processing via model.frame.
weights: optional numeric vector of case weights.
offset: optional numeric vector(s) with an a priori known component to be included in the linear predictor.
family: function that returns a "family" object, i.e., family(x) needs to be a "family" object when x is the numeric vector of extra parameters (by default assumed to be 1-dimensional).
xlink: link object or a character that can be passed to make.link. It should link the extra parameters to real parameters.
control: a list of control arguments as returned by glmx.control.
model, y, x: logicals. If TRUE the corresponding components of the fit (model frame, response, model matrix) are returned. For glmx.fit, x should be a numeric regressor matrix and y should be the response vector.
...: control arguments.

Details

The function glmx is a convenience interface that estimates generalized linear models (GLMs) with extra parameters. Examples would be binary response models with parametric link functions or count regression using a negative binomial family (which has one additional parameter).

Hence, glmx needs a family argument which is a family-generating function depending on one numeric argument for the extra parameters. Then, either profile-likelihood methods can be used for optimizing the extra parameters or all parameters can be optimized jointly.

If the generated family contains a list element loglik.extra for the derivative of the log-likelihood with respect to the extra parameters (i.e., score/gradient contributions), then this is used in the optimization process. This should be a function(y, mu, extra) depending on the observed response y, the estimated mean mu, and the extra parameters.

Value

glmx returns an object of class "glmx", i.e., a list with components as follows. glmx.fit returns an unclassed list with components up to converged.

coefficients: a list with elements "glm" and "extra" containing the coefficients from the respective models,
residuals: a vector of deviance residuals,
fitted.values: a vector of fitted means,
optim: list of optim outputs for maximizing the "profile" and "full" log-likelihood, respectively,
weights: the weights used (if any),
offset: the list of offset vectors used (if any),
n: number of observations,
nobs: number of observations with non-zero weights,
df: number of estimated parameters,
loglik: log-likelihood of the fitted model,
dispersion: estimate of the dispersion parameter (if any),
vcov: covariance matrix of all parameters in the model,
family: a list with elements "glm" and "extra" where the former contains the "family" object at the optimal extra parameters and the latter the family-generating function,
xlink: the link object for the extra parameters,
control: control options used,
converged: logical indicating successful convergence of optim,
call: the original function call,
formula: the formula,
terms: the terms object for the model,
levels: the levels of the categorical regressors,
contrasts: the contrasts corresponding to levels,
model: the full model frame (if model = TRUE),
y: the response vector (if y = TRUE),
x: the model matrix (if x = TRUE).

Examples

## artificial data from geometric regression
set.seed(1)
d <- data.frame(x = runif(200, -1, 1))
d$y <- rnbinom(200, mu = exp(0 + 3 * d$x), size = 1)

### negative binomial regression ###

## negative binomial regression via glmx
if(require("MASS")) {
m_nb1 <- glmx(y ~ x, data = d,
  family = negative.binomial, xlink = "log", xstart = 0)
summary(m_nb1)

## negative binomial regression via MASS::glm.nb
m_nb2 <- glm.nb(y ~ x, data = d)
summary(m_nb2)

## comparison
if(require("lmtest")) {
logLik(m_nb1)
logLik(m_nb2)
coeftest(m_nb1)
coeftest(m_nb2)
exp(coef(m_nb1, model = "extra"))
m_nb2$theta
exp(coef(m_nb1, model = "extra")) * sqrt(vcov(m_nb1, model = "extra"))
m_nb2$SE.theta
}}
#> Loading required package: MASS
#> [1] 0.251833

## if the score (or gradient) contribution of the extra parameters
## is supplied, then estimation can be speeded up:
negbin <- function(theta) {
  fam <- negative.binomial(theta)
  fam$loglik.extra <- function(y, mu, theta) digamma(y + theta) - digamma(theta) +
    log(theta) + 1 - log(mu + theta) - (y + theta)/(mu + theta)
  fam
}
m_nb3 <- glmx(y ~ x, data = d,
  family = negbin, xlink = "log", xstart = 0, profile = FALSE)
all.equal(coef(m_nb1), coef(m_nb3), tolerance = 1e-7)
#> [1] TRUE


### censored negative binomial hurdle regression (0 vs. > 0) ###

## negative binomial zero hurdle part via glmx
nbbin <- function(theta) binomial(link = nblogit(theta))
m_hnb1 <- glmx(factor(y > 0) ~ x, data = d,
  family = nbbin, xlink = "log", xstart = 0)
summary(m_hnb1)
#> 
#> Call:
#> glmx(formula = factor(y > 0) ~ x, data = d, family = nbbin, xlink = "log", 
#>     xstart = 0)
#> 
#> Deviance residuals:
#>    Min     1Q Median     3Q    Max 
#> 0.0254 0.2368 0.5597 1.3687 7.2594 
#> 
#> Coefficients (binomial model with nblogit link):
#>             Estimate Std. Error z value Pr(>|z|)   
#> (Intercept)  -0.1827     0.3417  -0.535  0.59275   
#> x             2.2709     0.7165   3.170  0.00153 **
#> 
#> Extra parameters (with log link):
#>            Estimate Std. Error z value Pr(>|z|)
#> log(theta)    1.412      2.345   0.602    0.547
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#> 
#> Log-likelihood: -95.02 on 3 Df
#> Dispersion: 1
#> Number of iterations in BFGS optimization: 9 (profile) 1 (full)

## negative binomial hurdle regression via pscl::hurdle
## (see only zero hurdle part)
if(require("pscl")) {
m_hnb2 <- hurdle(y ~ x, data = d, dist = "negbin", zero.dist = "negbin")
summary(m_hnb2)
}
#> Loading required package: pscl
#> Classes and Methods for R originally developed in the
#> Political Science Computational Laboratory
#> Department of Political Science
#> Stanford University (2002-2015),
#> by and under the direction of Simon Jackman.
#> hurdle and zeroinfl functions by Achim Zeileis.
#> 
#> Call:
#> hurdle(formula = y ~ x, data = d, dist = "negbin", zero.dist = "negbin")
#> 
#> Pearson residuals:
#>     Min      1Q  Median      3Q     Max 
#> -1.0850 -0.6256 -0.3790  0.4535  4.4274 
#> 
#> Count model coefficients (truncated negbin with log link):
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  -0.1178     0.2188  -0.538    0.590    
#> x             3.2573     0.3271   9.959   <2e-16 ***
#> Log(theta)    0.2066     0.2970   0.696    0.487    
#> Zero hurdle model coefficients (censored negbin with log link):
#>             Estimate Std. Error z value Pr(>|z|)   
#> (Intercept)  -0.1827     0.3416  -0.535  0.59274   
#> x             2.2709     0.7165   3.170  0.00153 **
#> Log(theta)    1.4115     2.3451   0.602  0.54725   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#> 
#> Theta: count = 1.2295, zero = 4.1021
#> Number of iterations in BFGS optimization: 29 
#> Log-likelihood: -342.8 on 6 Df