ptmixed (version 0.4.1)

nbmixed: Negative binomial generalized linear mixed model

Description

Estimates the negative binomial generalized linear mixed model with random intercept (here, the NB distribution is obtained as special case of the Poisson-Tweedie distribution when a = 0). Likelihood approximation for the model is based on the adaptive Gauss-Hermite quadrature rule.

Usage

nbmixed(fixef.formula, id, offset = NULL, data, npoints = 10,
  hessian = T, trace = T, theta.start = NULL, reltol = 1e-08,
  maxit = c(10000, 100), freq.updates = 200, min.var.init = 0.001)

Arguments

fixef.formula

A formula for the fixed effects part of the model. It should be in the form y ~ x1 + x2

id

A variable to distinguish observations from the same subject.

offset

An offset to be added to the linear predictor. Default is NULL.

data

A data frame containing the variables declared in fixef.formula.

npoints

Number of quadrature points employed in the adaptive quadrature. Default is 10.

hessian

Logical value. If TRUE, the hessian matrix is evaluated at the MLE to derive the observed Fisher information matrix. Default is TRUE.

trace

Logical value. If TRUE, additional information is printed during the optimization. Default is TRUE.

theta.start

Numeric vector comprising initial parameter values for the vector of regression coefficients, the dispersion parameter (using the same parametrization of ptmixed) and the variance of the random intercept. Default is NULL: initial parameter estimates are computed automatically by the function.

reltol

Relative tolerance to be used in optim. Default to 1e-8

maxit

Vector containing the maximum number of iterations used in optim by the Nelder-Mead method and, if this fails, by the BFGS method

freq.updates

Number of iterations after which the quadrature points are updated when the Nelder-Mead algorithm is used for the optimization. Default value is 200. To update the quadrature points at every iteration (note that this typically makes the computation 2.5x slower), set freq.updates = 1 or freq.updates = NA. The function first tries to optimize the loglikelihood using the Nelder-Mead algorithm, updating the quadrature points every freq.updates iterations. If this fails to converge, a second attempt is made using the BFGS algorithm, for which the quadrature points are updated at every iteration.

min.var.init

If the initial estimate of the variance of the random intercept is smaller than this value, estimation is stopped and the user is advided to use the simpler Poisson-Tweedie GLM is used. Default is 1e-3.

Value

A list containing the following elements: function's call (call); maximum likelihood estimate (mle); value of the loglikelihood at the mle (logl); convergence value (if 0, the optimization converged); the observed Fisher information (fisher.info), if hessian = T; the number of quadrature points used (quad.points) and the starting value used in the optimization (theta.init); relevant warnings (warnings).

See Also

summary.ptglmm, ranef

Examples

Run this code
# NOT RUN {
# generate data
set.seed(123)
n = 6; t = 3
id = rep(1:n, each = t)
rand.int = rep(rnorm(n, sd = 0.7), each = t)
group = rep(c(0,1), each = n*t/2)
time = rep(0:(t-1), n)
offset = rnorm(n*t, sd = 0.3)

beta = c(3, 0.3, 0.1)
X = model.matrix(~group + time)
mu = exp(X %*% beta + rand.int + offset)
y = rep(NA, n*t)
library(tweeDEseq)
for (i in 1:(n*t)) y[i] = rPT(1, mu = mu[i], D = 2, a = 0, max = 1000)

data.long = data.frame(y, group, time, id, offset)
rm(list = setdiff(ls(), 'data.long'))

# 1) Quick example (5 quadrature points, hessian and SEs not computed)
# estimate the model
fit1 = nbmixed(fixef.formula = y ~ group + time, id = data.long$id,
              offset = data.long$offset, data = data.long, npoints = 5, 
              freq.updates = 200, hessian = FALSE, trace = TRUE)
# print summary:
summary(fit1, wald = FALSE)

# }
# NOT RUN {
# 2) Full computation, including hessian evaluation and using more quadrature points
# estimate the model
fit2 = nbmixed(fixef.formula = y ~ group + time, id = data.long$id,
              offset = data.long$offset, data = data.long, npoints = 10, 
              freq.updates = 200, hessian = TRUE, trace = TRUE)
# print and get summary:
results = summary(fit2, wald = TRUE)
ls(results)
# view table with estimates of regression coefficients, standard errors and Wald test:
results$coefficients
# }

Run the code above in your browser using DataCamp Workspace