glmm.admb: Generalized Linear Mixed Models Using AD Model Builder

In bbolker/glmmadmb: Generalized Linear Mixed Models using 'AD Model Builder'

Description

Fits (extended) generalized linear mixed-effects models to data using a variety of distributions and link functions, including zero-inflated models.

Usage

glmmadmb(formula, data, family="poisson", link,start,
          random, corStruct="diag", easyFlag=TRUE,
          zeroInflation=FALSE, 
          admb.opts=admbControl(),
          mcmc=FALSE,mcmc.opts=mcmcControl(),save.dir,
          verbose=FALSE,extra.args="",
          bin_loc=NULL,debug=FALSE)
glmm.admb(...)  ## deprecated

Arguments

formula	a two-sided linear formula describing the part of the model. The response, which may be a numeric vector (for most models) or a matrix (for binomial models) is on the left of a ~ operator and the terms, separated by + operators, are on the right. Random effects terms should be enclosed in parentheses (x|g), where x represents the formula for the parameter(s) that vary across grouping variables and g represents the grouping variables; for example, (1|g) indicates a model where only intercepts vary across groups. Grouping variables must be factors. The model may include an offset term, which must be specified on the scale of the linear predictor and included in data: see offset. Alternatively, the random effects terms may be specified separately using the random argument.
data	a data frame containing the variables named in formula.
family	a string determining the response distribution: currently implemented options (and default link functions) are exponential family distributions: "poisson" (log), "binomial" (logit), "gamma" (log), "gaussian" (identity) negative binomial (log): "nbinom", "nbinom1" (negative binomial “type 1”: see Details), "nbinom2" (a synonym for "nbinom") extended distributions: "beta" (logit), "logistic" (identity), "betabinomial" (logit) truncated distributions (log): "truncpoiss", "truncnbinom", "truncnbinom1" (these can be used to fit hurdle models)
link	a string specifying the link function: currently implemented options are "log", "logit", "probit", "cloglog", and "identity".
start	an (optional) named list specifying starting values, with (optional) components fixed (fixed-effect parameters: default 0); pz (zero-inflation probability: default 0.02); log_alpha (log of dispersion parameter: default 1); RE_sd (log-standard deviations of random effects: default 0.25); RE_cor (off-diagonal elements of Cholesky decomposition of random-effect correlation matrix: default 0.0001); u (random effects [unscaled, i.e. not multiplied by corresponding RE standard deviation]: default 0)
random	a formula specifying the random effects. A single random effect may be specified as ~x|g, but parentheses must be used if more than one random effect is specified (e.g., ~(x1|g1)+(x2|g2)).
corStruct	a string specifying the covariance structure of the random effects vector. Two types of covariance matrices are currently implemented: "diag" (diagonal matrix) and "full" (positive definite matrix with all elements being estimated). At present the same covariance structure must be specified for all grouping variables in the model; at some point in the future this requirement may be relaxed.
easyFlag	(currently inactive: whether a faster but less robust optimization algorithm should be employed (only "poisson" and "nbinom" families).
zeroInflation	whether a zero-inflated model should be fitted (only "poisson" and "nbinom" families).
admb.opts	options for AD Model Builder: see admbControl
mcmc	(logical) run MCMC after fitting, based on estimated mode and parameter variance-covariance matrix
mcmc.opts	list of options for MCMC run: see mcmcControl
save.dir	Optionally, a directory name where all ADMB output files are saved. If save.dir is not specified, a temporary directory is created (with tempfile) and all ADMB input and output files are deleted after the ADMB run finishes and the results are read into R.
verbose	(logical) Print interim output from AD Model Builder run?
extra.args	(character) extra arguments to AD Model Builder (see ADMB manual)
bin_loc	location of glmmadmb binary file (normally determined automatically by glmmADMB:::get_bin_loc)
debug	(logical) print debugging statements?
...	arguments to pass through to glmmadmb

Details

• For advice on troubleshooting, see admbControl.

• Parameterization of the negative binomial distribution:

  • for family=="nbinom":

    Var(Y) = E(Y) * (1 + E(Y)/alpha)

    (i.e. Hardin and Hilbe NB2; alpha corresponds to the size parameter in dnbinom)

  • for family=="nbinom1":

    Var(Y) = E(Y) * alpha

    (i.e. Hardin and Hilbe NB1: alpha corresponds to the scale parameter in a quasi-Poisson model)

• Parameterization of the beta-binomial distribution: as in Morris (1997), in terms of N (number of trials), p (per-trial probability, determined by the inverse-link function applied to the linear predictor), alpha (overdispersion parameter: as alpha becomes large, the distribution converges to binomial). The density function is

  p(x) = (C(N,x)*Beta(N-x+theta*(1-p),x+theta*p))/Beta(theta*(1-p),theta*p)

  .

• Zero-inflation: With probability 1-pz, Y comes from a Poisson (or negative binomial) distribution, and with probability pz, Y is zero (Bohning et al. 1999).

• Parameters are estimated by maximum likelihood, using the Laplace approximation to evaluate the marginal likelihood. When impSamp >0, importance sampling is used to improve the Laplace approximation (Skaug and Fournier 2006).

• Due to technical limitations, offset expressions can contain at most one internal set of parentheses; i.e. offset(log(x)) works, but offset((log(x))) will not. If necessary, define a new variable in the data frame, e.g. mydata$myoffset <- with(mydata,log(3*exp(x))).

• Internally, glmmADMB works with an orthogonalized version of the fixed-effect design matrix. The phi matrix returned in the results can be used to convert from orthogonalized to 'real' fixed parameters, i.e. b_real == b_orth %*% phi or b_orth = b_real %*% solve(phi). Fixed effects specified in the starting values are automatically orthogonalized.

• If you are specifying starting values from a different GLMM fit for more than the fixed-effect parameters, you may want to use extra.args="-phase 5" to skip over the estimation phase where ADMB fits the model without random effects.

Value

An object of class "glmmadmb" representing the model fit, including (among others) components:

b	vector of fixed effects
S	covariance matrix of random effects
alpha	scale/overdispersion parameter (negative binomial, Gamma, beta)
pz	Zero-inflation parameter (only when zeroInflation=TRUE)
phi	Matrix for converting from 'orthogonalized' to 'real' parameters: see Details

Note

• glmm.admb is defined for backward compatibility; it will be removed in a future version.

• Tools for working with MCMC output are preliminary; see mcmcControl or the vignette for more information.

References

Bohning, D., E. Dietz, and P. Schlattmann. 1999. The zero-inflated Poisson model and the decayed, missing and filled teeth index in dental epidemiology. Journal of the Royal Statistical Society. Series A (Statistics in Society) 162:195–209.

Skaug, H.J. and D.A. Fournier. 2006. Automatic approximation of the marginal likelihood in non-Gaussian hierarchical models. Computational Statistics & Data Analysis 51:699–709.

Morris, W. 1997. American Naturalist 150:299-327

See Also

glmmADMB-package.

Examples

## Australian H. influenzae data from MASS
data(bacteria,package="MASS")
bacteria$present <- as.numeric(bacteria$y)-1
if (!check_rforge())
    (bfit <-  glmmadmb(present ~ trt + I(week > 2), random = ~ 1 | ID,
                       family = "binomial", data = bacteria))
## simple simulated zero-inflated Poisson example
### simulate values
set.seed(101)
d <- data.frame(f=factor(rep(LETTERS[1:10],each=10)),x=runif(100))
u <- rnorm(10,sd=2)
d$eta <- with(d,u[f]+1+4*x)
pz <- 0.3
zi <- rbinom(100,size=1,prob=pz)
d$y <- ifelse(zi,0,rpois(100,lambda=exp(d$eta)))
## fit
## Not run: 
    zipmodel <- glmmadmb(y~x+(1|f),data=d,family="poisson",zeroInflation=TRUE)

## End(Not run)

## Not run: 
## Epilepsy data
## all grouping variables must be factors:
  epil2$subject <- factor(epil2$subject)
  (fm <- glmmadmb(y~Base*trt+Age+Visit+(Visit|subject),
                  data=epil2, family="nbinom"))
## Owls data
 om <- glmmadmb(SiblingNegotiation~FoodTreatment*SexParent+
           (1|Nest)+offset(log(BroodSize)),
          zeroInflation=TRUE,family="nbinom",data=Owls)

## End(Not run)

bbolker/glmmadmb documentation built on May 11, 2019, 9:29 p.m.