Package 'AICcmodavg'

Model Selection and Multimodel Inference Based on (Q)AIC(c)

Description

Description: This package includes functions to create model selection tables based on Akaike's information criterion (AIC) and the second-order AIC (AICc), as well as their quasi-likelihood counterparts (QAIC, QAICc). The package also features functions to conduct classic model averaging (multimodel inference) for a given parameter of interest or predicted values, as well as a shrinkage version of model averaging parameter estimates. Other handy functions enable the computation of relative variable importance, evidence ratios, and confidence sets for the best model. The present version supports Cox proportional hazards models and conditional logistic regression (coxph and coxme classes), linear models (lm class), generalized linear models (glm, glm.nb, vglm, hurdle, and zeroinfl classes), linear models fit by generalized least squares (gls class), linear mixed models (lme class), generalized linear mixed models (mer, merMod, and glmmTMB classes), multinomial and ordinal logistic regressions (multinom, polr, clm, and clmm classes), robust regression models (rlm class), beta regression models (betareg class), parametric survival models (survreg class), nonlinear models (nls and gnls classes), nonlinear mixed models (nlme and nlmerMod classes), univariate models (fitdist and fitdistr classes), and certain types of latent variable models (lavaan class). The package also supports various models of unmarkedFit and maxLikeFit classes estimating demographic parameters after accounting for imperfect detection probabilities. Some functions also allow the creation of model selection tables for Bayesian models of the bugs and rjags classes. Objects following model selection and multimodel inference can be formatted to LaTeX using xtable methods included in the package.

Details

Package:	AICcmodavg
Type:	Package
Version:	2.3-3
Date:	2023-11-16
License:	GPL (>=2 )
LazyLoad:	yes

Many functions of the package require a list of models as the input to conduct model selection and multimodel inference. Thus, you should start by organizing the output of the models in a list (See 'Examples' below).

This package contains several useful functions for model selection and multimodel inference for several model classes:

AICc

Computes AIC, AICc, and their quasi-likelihood counterparts (QAIC, QAICc).

aictab

Constructs model selection tables with number of parameters, AIC, delta AIC, Akaike weights or variants based on AICc, QAIC, and QAICc for a set of candidate models.

bictab

Constructs model selection tables with number of parameters, BIC, delta BIC, BIC weights for a set of candidate models.

boot.wt

Computes summary statistics from detection histories.

confset

Determines the confidence set for the best model based on one of three criteria.

DIC

Extracts DIC.

dictab

Constructs model selection tables with number of parameters, DIC, delta DIC, DIC weights for a set of candidate models.

evidence

Computes the evidence ratio between the highest-ranked model based on the information criteria selected and a lower-ranked model.

importance

Computes importance values (w+) for the support of a given parameter among set of candidate models.

modavg

Computes model-averaged estimate, unconditional standard error, and unconditional confidence interval of a parameter of interest among a set of candidate models.

modavgEffect

Computes model-averaged effect sizes between groups based on the entire candidate model set.

modavgShrink

Computes shrinkage version of model-averaged estimate, unconditional standard error, and unconditional confidence interval of a parameter of interest among entire set of candidate models.

modavgPred

Computes model-average predictions, unconditional SE's, and confidence intervals among entire set of candidate models.

multComp

Performs multiple comparisons across levels of a factor in a model selection framework.

useBIC

Computes BIC or a quasi-likelihood counterparts (QBIC).

For models not yet supported by the functions above, the following can be useful for model selection and multimodel inference conducted from input values supplied by the user:

AICcCustom

Computes AIC, AICc, QAIC, and QAICc from user-supplied input values of log-likelihood and number of parameters.

aictabCustom

Creates model selection tables based on (Q)AIC(c) from user-supplied input values of log-likelihood and number of parameters.

bictabCustom

Creates model selection tables based on (Q)BIC from user-supplied input values of log-likelihood and number of parameters.

ictab

Creates model selection tables from user-supplied values of an information criterion.

modavgCustom

Computes model-averaged parameter estimate based on (Q)AIC(c) from user-supplied input values of log-likelihood, number of parameters, parameter estimates, and standard errors.

modavgIC

Computes model-averaged parameter estimate from user-supplied values of information criterion, parameter estimates, and standard errors.

useBICCustom

Computes BIC and QBIC from user-supplied input values of log-likelihoods and number of parameters.

A number of functions for model diagnostics are available:

c_hat

Estimates variance inflation factor for binomial or Poisson GLM's based on various estimators.

checkConv

Checks the convergence information of the algorithm for the model.

checkParms

Checks the occurrence of parameter estimates with high standard errors in a model.

countDist

Computes summary statistics from distance sampling data.

countHist

Computes summary statistics from count history data.

covDiag

Computes covariance diagnostics for lambda in N-mixture models.

detHist

Computes summary statistics from detection histories.

detTime

Computes summary statistics from time-to-detection data.

extractCN

Extracts condition number from models of certain classes.

mb.gof.test

Computes the MacKenzie and Bailey goodness-of-fit test for single season and dynamic occupancy models using the Pearson chi-square statistic.

Nmix.gof.test

Computes goodness-of-fit test for N-mixture models based on the Pearson chi-square statistic.

Other utility functions include:

anovaOD

Computes likelihood-ratio test statistic corrected for overdispersion between two models.

extractLL

Extracts log-likelihood from models of certain classes.

extractSE

Extracts standard errors from models of certain classes and adds the labels.

extractX

Extracts the predictors and associated information on variables from a list of candidate models.

fam.link.mer

Extracts the distribution family and link function from a generalized linear mixed model of classes mer and merMod.

predictSE

Computes predictions and associated standard errors models of certain classes.

summaryOD

Displays summary of model output adjusted for overdispersion.

xtable

Formats various objects resulting from model selection and multimodel inference to LaTeX or HTML tables.

Author(s)

Marc J. Mazerolle <[email protected]>.

References

Anderson, D. R. (2008) Model-based inference in the life sciences: a primer on evidence. Springer: New York.

Burnham, K. P., and Anderson, D. R. (2002) Model selection and multimodel inference: a practical information-theoretic approach. Second edition. Springer: New York.

Burnham, K. P., Anderson, D. R. (2004) Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods and Research 33, 261–304.

Mazerolle, M. J. (2006) Improving data analysis in herpetology: using Akaike's Information Criterion (AIC) to assess the strength of biological hypotheses. Amphibia-Reptilia 27, 169–180.

Examples

##Example 1:  Poisson GLM with offset
##anuran larvae example from Mazerolle (2006) 
data(min.trap)
##assign "UPLAND" as the reference level as in Mazerolle (2006)          
min.trap$Type <- relevel(min.trap$Type, ref = "UPLAND") 

##set up candidate models in a list
Cand.mod <- list()
##global model          
Cand.mod[[1]] <- glm(Num_anura ~ Type + log.Perimeter + Num_ranatra,
                     family = poisson, offset = log(Effort),
                     data = min.trap) 
Cand.mod[[2]] <- glm(Num_anura ~ Type + log.Perimeter, family = poisson,
                     offset = log(Effort), data = min.trap) 
Cand.mod[[3]] <- glm(Num_anura ~ Type + Num_ranatra, family = poisson,
                     offset = log(Effort), data = min.trap) 
Cand.mod[[4]] <- glm(Num_anura ~ Type, family = poisson,
                     offset = log(Effort), data = min.trap) 
Cand.mod[[5]] <- glm(Num_anura ~ log.Perimeter + Num_ranatra,
                     family = poisson, offset = log(Effort),
                     data = min.trap) 
Cand.mod[[6]] <- glm(Num_anura ~ log.Perimeter, family = poisson,
                     offset = log(Effort), data = min.trap) 
Cand.mod[[7]] <- glm(Num_anura ~ Num_ranatra, family = poisson,
                     offset = log(Effort), data = min.trap) 
Cand.mod[[8]] <- glm(Num_anura ~ 1, family = poisson,
                     offset = log(Effort), data = min.trap) 
          
##check c-hat for global model
c_hat(Cand.mod[[1]], method = "pearson") #uses Pearson's chi-square/df
##note the very low overdispersion: in this case, the analysis could be
##conducted without correcting for c-hat as its value is reasonably close
##to 1  

##output of model corrected for overdispersion
summaryOD(Cand.mod[[1]], c.hat = 1.04)

##assign names to each model
Modnames <- c("type + logperim + invertpred", "type + logperim",
              "type + invertpred", "type", "logperim + invertpred",
              "logperim", "invertpred", "intercept only") 

##model selection table based on AICc
aictab(cand.set = Cand.mod, modnames = Modnames)

##compute evidence ratio
evidence(aictab(cand.set = Cand.mod, modnames = Modnames))

##compute confidence set based on 'raw' method
confset(cand.set = Cand.mod, modnames = Modnames, second.ord = TRUE,
        method = "raw")  

##compute importance value for "TypeBOG" - same number of models
##with vs without variable
importance(cand.set = Cand.mod, modnames = Modnames, parm = "TypeBOG") 

##compute model-averaged estimate of "TypeBOG" using the natural average
modavg(cand.set = Cand.mod, modnames = Modnames, parm = "TypeBOG")

##compute model-averaged estimate of "TypeBOG" using shrinkage estimator
##same number of models with vs without variable
modavgShrink(cand.set = Cand.mod, modnames = Modnames,
             parm = "TypeBOG")

##compute model-averaged predictions for two types of ponds
##create a data set for predictions
dat.pred <- data.frame(Type = factor(c("BOG", "UPLAND")),
                       log.Perimeter = mean(min.trap$log.Perimeter),
                       Num_ranatra = mean(min.trap$Num_ranatra),
                       Effort = mean(min.trap$Effort))

##model-averaged predictions across entire model set
modavgPred(cand.set = Cand.mod, modnames = Modnames,
           newdata = dat.pred, type = "response")

##compute model-averaged effect size between two groups
##'newdata' data frame must be limited to two rows
modavgEffect(cand.set = Cand.mod, modnames = Modnames,
             newdata = dat.pred, type = "link")


## Not run: 
##Example 2:  single-season occupancy model example modified from ?occu
require(unmarked)
##single season
data(frogs)
pferUMF <- unmarkedFrameOccu(pfer.bin)
## add some fake covariates for illustration
siteCovs(pferUMF) <- data.frame(sitevar1 = rnorm(numSites(pferUMF)),
                                sitevar2 = rnorm(numSites(pferUMF))) 
     
## observation covariates are in site-major, observation-minor order
obsCovs(pferUMF) <- data.frame(obsvar1 = rnorm(numSites(pferUMF) *
                                 obsNum(pferUMF))) 

##check detection history data from data object
detHist(pferUMF)

##set up candidate model set
fm1 <- occu(~ obsvar1 ~ sitevar1, pferUMF)
##check detection history data from model object
detHist(fm1)

fm2 <- occu(~ 1 ~ sitevar1, pferUMF)
fm3 <- occu(~ obsvar1 ~ sitevar2, pferUMF)
fm4 <- occu(~ 1 ~ sitevar2, pferUMF)
Cand.models <- list(fm1, fm2, fm3, fm4)

##assign names to elements in list
##alternative to using 'modnames' argument
names(Cand.models) <- c("fm1", "fm2", "fm3", "fm4")

##check GOF of global model and estimate c-hat
mb.gof.test(fm4, nsim = 100) #nsim should be > 1000

##check for high SE's in models
lapply(Cand.models, checkParms, simplify = FALSE)

##compute table
print(aictab(cand.set = Cand.models,
             second.ord = TRUE), digits = 4)

##export as LaTeX table
if(require(xtable)) {
xtable(aictab(cand.set = Cand.models,
              second.ord = TRUE))
}

##compute evidence ratio
evidence(aictab(cand.set = Cand.models))
##evidence ratio between top model vs lowest-ranked model
evidence(aictab(cand.set = Cand.models), model.high = "fm2", model.low = "fm3")

##compute confidence set based on 'raw' method
confset(cand.set = Cand.models, second.ord = TRUE,
        method = "raw")  

##compute importance value for "sitevar1" on occupancy
##same number of models with vs without variable
importance(cand.set = Cand.models, parm = "sitevar1",
           parm.type = "psi") 

##compute model-averaged estimate of "sitevar1" on occupancy
##(natural average)
modavg(cand.set = Cand.models, parm = "sitevar1",
       parm.type = "psi")

##compute model-averaged estimate of "sitevar1"
##(shrinkage estimator)
##same number of models with vs without variable
modavgShrink(cand.set = Cand.models,
             parm = "sitevar1", parm.type = "psi")

##compute model-average predictions
##check explanatory variables appearing in models
extractX(Cand.models, parm.type = "psi")

##create a data set for predictions
dat.pred <- data.frame(sitevar1 = seq(from = min(siteCovs(pferUMF)$sitevar1),
                         to = max(siteCovs(pferUMF)$sitevar1), by = 0.5),
                       sitevar2 = mean(siteCovs(pferUMF)$sitevar2))

##model-averaged predictions of psi across range of values
##of sitevar1 and entire model set
modavgPred(cand.set = Cand.models, newdata = dat.pred,
           parm.type = "psi")
detach(package:unmarked)

## End(Not run)

## Not run: 
##Example 3:  example with user-supplied values of log-likelihoods and
##number of parameters

##vector with model LL's
LL <- c(-38.8876, -35.1783, -64.8970)
     
##vector with number of parameters
Ks <- c(7, 9, 4)
     
##create a vector of names to trace back models in set
Modnames <- c("Cm1", "Cm2", "Cm3")

##generate AICc table
aictabCustom(logL = LL, K = Ks, modnames = Modnames, nobs = 121,
             sort = TRUE) 
##generate AIC table
aictabCustom(logL = LL, K = Ks, modnames = Modnames,
             second.ord = FALSE, nobs = 121, sort = TRUE)

##model averaging parameter estimate
##vector of beta estimates for a parameter of interest
model.ests <- c(0.0478, 0.0480, 0.0478)
     
##vector of SE's of beta estimates for a parameter of interest
model.se.ests <- c(0.0028, 0.0028, 0.0034)
     
##compute model-averaged estimate and unconditional SE based on AICc
modavgCustom(logL = LL, K = Ks, modnames = Modnames, 
             estimate = model.ests, se = model.se.ests, nobs = 121)
##compute model-averaged estimate and unconditional SE based on BIC
modavgCustom(logL = LL, K = Ks, modnames = Modnames, 
             estimate = model.ests, se = model.se.ests, nobs = 121,
             useBIC = TRUE)

## End(Not run)

## Not run: 
##Example 4:  example with user-supplied values of information criterion
##model selection based on WAIC

##WAIC values
waic <- c(105.74, 107.36, 108.24, 100.57)
##number of effective parameters
effK <- c(7.45, 5.61, 6.14, 6.05)

##create a vector of names to trace back models in set
Modnames <- c("global model", "interactive model",
              "additive model", "invertpred model")

##generate WAIC model selection table
ictab(ic = waic, K = effK, modnames = Modnames,
      sort = TRUE, ic.name = "WAIC")

##compute model-averaged estimate
##vector of predictions
Preds <- c(0.106, 0.137, 0.067, 0.050)
##vector of SE's for prediction
Ses <- c(0.128, 0.159, 0.054, 0.039)

##compute model-averaged estimate and unconditional SE based on WAIC
modavgIC(ic = waic, K = effK, modnames = Modnames, 
         estimate = Preds, se = Ses,
         ic.name = "WAIC")

##export as LaTeX table
if(require(xtable)) {
##model-averaged estimate and confidence interval
xtable(modavgIC(ic = waic, K = effK, modnames = Modnames, 
       estimate = Preds, se = Ses,
       ic.name = "WAIC"))
##model selection table with estimate and SE's from each model
xtable(modavgIC(ic = waic, K = effK, modnames = Modnames, 
       estimate = Preds, se = Ses,
       ic.name = "WAIC"), print.table = TRUE)
}

## End(Not run)

---

Computing AIC, AICc, QAIC, and QAICc

Description

Functions to compute Akaike's information criterion (AIC), the second-order AIC (AICc), as well as their quasi-likelihood counterparts (QAIC, QAICc).

Usage

AICc(mod, return.K = FALSE, second.ord = TRUE, nobs = NULL, ...) 

## S3 method for class 'aov'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'betareg'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'clm'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...) 

## S3 method for class 'clmm'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)  

## S3 method for class 'coxme'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...) 

## S3 method for class 'coxph'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...) 

## S3 method for class 'fitdist'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'fitdistr'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'glm'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, c.hat = 1, ...)

## S3 method for class 'glmmTMB'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, c.hat = 1, ...)

## S3 method for class 'gls'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'gnls'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'hurdle'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'lavaan'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'lm'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'lme'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'lmekin'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...) 

## S3 method for class 'maxlikeFit'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, c.hat = 1, ...)

## S3 method for class 'mer'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'merMod'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...) 

## S3 method for class 'lmerModLmerTest'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...) 

## S3 method for class 'multinom'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, c.hat = 1, ...)

## S3 method for class 'negbin'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'nlme'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'nls'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'polr'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'rlm'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'survreg'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

## S3 method for class 'unmarkedFit'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, c.hat = 1, ...)

## S3 method for class 'vglm'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, c.hat = 1, ...)

## S3 method for class 'zeroinfl'
AICc(mod, return.K = FALSE, second.ord = TRUE,
     nobs = NULL, ...)

Arguments

mod	an object of class aov, betareg, clm, clmm, clogit, coxme, coxph, fitdist, fitdistr, glm, glmmTMB, gls, gnls, hurdle, lavaan, lm, lme, lmekin, maxlikeFit, mer, merMod, lmerModLmerTest, multinom, negbin, nlme, nls, polr, rlm, survreg, vglm, zeroinfl, and various unmarkedFit classes containing the output of a model.
return.K	logical. If FALSE, the function returns the information criterion specified. If TRUE, the function returns K (number of estimated parameters) for a given model.
second.ord	logical. If TRUE, the function returns the second-order Akaike information criterion (i.e., AICc).
nobs	this argument allows to specify a numeric value other than total sample size to compute the AICc (i.e., nobs defaults to total number of observations). This is relevant only for mixed models or various models of unmarkedFit classes where sample size is not straightforward. In such cases, one might use total number of observations or number of independent clusters (e.g., sites) as the value of nobs.
c.hat	value of overdispersion parameter (i.e., variance inflation factor) such as that obtained from c_hat. Note that values of c.hat different from 1 are only appropriate for binomial GLM's with trials > 1 (i.e., success/trial or cbind(success, failure) syntax), with Poisson GLM's, single-season occupancy models (MacKenzie et al. 2002), dynamic occupancy models (MacKenzie et al. 2003), or N-mixture models (Royle 2004, Dail and Madsen 2011). If c.hat > 1, AICc will return the quasi-likelihood analogue of the information criteria requested and multiply the variance-covariance matrix of the estimates by this value (i.e., SE's are multiplied by sqrt(c.hat)). This option is not supported for generalized linear mixed models of the mer or merMod classes.
...	additional arguments passed to the function.

Details

AICc computes one of the following four information criteria:

Akaike's information criterion (AIC, Akaike 1973),

$-2 * log-likelihood + 2 * K,$

where the log-likelihood is the maximum log-likelihood of the model and K corresponds to the number of estimated parameters.

Second-order or small sample AIC (AICc, Sugiura 1978, Hurvich and Tsai 1989, 1991),

$-2 * log-likelihood + 2 * K * (n/(n - K - 1)),$

where n is the sample size of the data set.

Quasi-likelihood AIC (QAIC, Burnham and Anderson 2002),

$QAIC = \frac{-2 * log-likelihood}{c-hat} + 2 * K,$

where c-hat is the overdispersion parameter specified by the user with the argument c.hat.

Quasi-likelihood AICc (QAICc, Burnham and Anderson 2002),

$QAIC = \frac{-2 * log-likelihood}{c-hat} + 2 * K * (n/(n - K - 1))$

.

Note that AIC and AICc values are meaningful to select among gls or lme models fit by maximum likelihood. AIC and AICc based on REML are valid to select among different models that only differ in their random effects (Pinheiro and Bates 2000).

Value

AICc returns the AIC, AICc, QAIC, or QAICc, or the number of estimated parameters, depending on the values of the arguments.

Note

The actual (Q)AIC(c) values are not really interesting in themselves, as they depend directly on the data, parameters estimated, and likelihood function. Furthermore, a single value does not tell much about model fit. Information criteria become relevant when compared to one another for a given data set and set of candidate models.

Author(s)

Marc J. Mazerolle

References

Akaike, H. (1973) Information theory as an extension of the maximum likelihood principle. In: Second International Symposium on Information Theory, pp. 267–281. Petrov, B.N., Csaki, F., Eds, Akademiai Kiado, Budapest.

Anderson, D. R. (2008) Model-based Inference in the Life Sciences: a primer on evidence. Springer: New York.

Burnham, K. P., Anderson, D. R. (2002) Model Selection and Multimodel Inference: a practical information-theoretic approach. Second edition. Springer: New York.

Burnham, K. P., Anderson, D. R. (2004) Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods and Research 33, 261–304.

Dail, D., Madsen, L. (2011) Models for estimating abundance from repeated counts of an open population. Biometrics 67, 577–587.

Hurvich, C. M., Tsai, C.-L. (1989) Regression and time series model selection in small samples. Biometrika 76, 297–307.

Hurvich, C. M., Tsai, C.-L. (1991) Bias of the corrected AIC criterion for underfitted regression and time series models. Biometrika 78, 499–509.

MacKenzie, D. I., Nichols, J. D., Lachman, G. B., Droege, S., Royle, J. A., Langtimm, C. A. (2002) Estimating site occupancy rates when detection probabilities are less than one. Ecology 83, 2248–2255.

MacKenzie, D. I., Nichols, J. D., Hines, J. E., Knutson, M. G., Franklin, A. B. (2003) Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84, 2200–2207.

Pinheiro, J. C., Bates, D. M. (2000) Mixed-effect models in S and S-PLUS. Springer Verlag: New York.

Royle, J. A. (2004) N-mixture models for estimating population size from spatially replicated counts. Biometrics 60, 108–115.

Sugiura, N. (1978) Further analysis of the data by Akaike's information criterion and the finite corrections. Communications in Statistics: Theory and Methods A7, 13–26.

See Also

AICcCustom, aictab, confset, importance, evidence, c_hat, modavg, modavgShrink, modavgPred, useBIC,

Examples

##cement data from Burnham and Anderson (2002, p. 101)
data(cement)
##run multiple regression - the global model in Table 3.2
glob.mod <- lm(y ~ x1 + x2 + x3 + x4, data = cement)

##compute AICc with full likelihood
AICc(glob.mod, return.K = FALSE)

##compute AIC with full likelihood 
AICc(glob.mod, return.K = FALSE, second.ord = FALSE)
##note that Burnham and Anderson (2002) did not use full likelihood
##in Table 3.2 and that the MLE estimate of the variance was
##rounded to 2 digits after decimal point  


##compute AICc for mixed model on Orthodont data set in Pinheiro and
##Bates (2000)
## Not run: 
require(nlme)
m1 <- lme(distance ~ age, random = ~1 | Subject, data = Orthodont,
          method= "ML")
AICc(m1, return.K = FALSE)

## End(Not run)

---

Compute AIC, AICc, QAIC, and QAICc from User-supplied Input

Description

This function computes Akaike's information criterion (AIC), the second-order AIC (AICc), as well as their quasi-likelihood counterparts (QAIC, QAICc) from user-supplied input instead of extracting the values automatically from a model object. This function is particularly useful for output imported from other software or for model classes that are not currently supported by AICc.

Usage

AICcCustom(logL, K, return.K = FALSE, second.ord = TRUE, nobs = NULL,
           c.hat = 1)

Arguments

logL	the value of the model log-likelihood.
K	the number of estimated parameters in the model.
return.K	logical. If FALSE, the function returns the information criterion specified. If TRUE, the function returns K (number of estimated parameters) for a given model.
second.ord	logical. If TRUE, the function returns the second-order Akaike information criterion (i.e., AICc).
nobs	the sample size required to compute the AICc or QAICc.
c.hat	value of overdispersion parameter (i.e., variance inflation factor) such as that obtained from c_hat. Note that values of c.hat different from 1 are only appropriate for binomial GLM's with trials > 1 (i.e., success/trial or cbind(success, failure) syntax), with Poisson GLM's, single-season or dynamic occupancy models (MacKenzie et al. 2002, 2003), N-mixture models (Royle 2004, Dail and Madsen 2011), or capture-mark-recapture models (e.g., Lebreton et al. 1992). If c.hat > 1, AICcCustom will return the quasi-likelihood analogue of the information criterion requested.

Details

AICcCustom computes one of the following four information criteria:

Akaike's information criterion (AIC, Akaike 1973), the second-order or small sample AIC (AICc, Sugiura 1978, Hurvich and Tsai 1989, 1991), the quasi-likelihood AIC (QAIC, Burnham and Anderson 2002), and the quasi-likelihood AICc (QAICc, Burnham and Anderson 2002).

Value

AICcCustom returns the AIC, AICc, QAIC, or QAICc, or the number of estimated parameters, depending on the values of the arguments.

Note

The actual (Q)AIC(c) values are not really interesting in themselves, as they depend directly on the data, parameters estimated, and likelihood function. Furthermore, a single value does not tell much about model fit. Information criteria become relevant when compared to one another for a given data set and set of candidate models.

Author(s)

Marc J. Mazerolle

References

Akaike, H. (1973) Information theory as an extension of the maximum likelihood principle. In: Second International Symposium on Information Theory, pp. 267–281. Petrov, B.N., Csaki, F., Eds, Akademiai Kiado, Budapest.

Burnham, K. P., Anderson, D. R. (2002) Model Selection and Multimodel Inference: a practical information-theoretic approach. Second edition. Springer: New York.

Dail, D., Madsen, L. (2011) Models for estimating abundance from repeated counts of an open population. Biometrics 67, 577–587.

Hurvich, C. M., Tsai, C.-L. (1989) Regression and time series model selection in small samples. Biometrika 76, 297–307.

Hurvich, C. M., Tsai, C.-L. (1991) Bias of the corrected AIC criterion for underfitted regression and time series models. Biometrika 78, 499–509.

Lebreton, J.-D., Burnham, K. P., Clobert, J., Anderson, D. R. (1992) Modeling survival and testing biological hypotheses using marked animals: a unified approach with case-studies. Ecological Monographs 62, 67–118.

MacKenzie, D. I., Nichols, J. D., Lachman, G. B., Droege, S., Royle, J. A., Langtimm, C. A. (2002) Estimating site occupancy rates when detection probabilities are less than one. Ecology 83, 2248–2255.

MacKenzie, D. I., Nichols, J. D., Hines, J. E., Knutson, M. G., Franklin, A. B. (2003) Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84, 2200–2207.

Royle, J. A. (2004) N-mixture models for estimating population size from spatially replicated counts. Biometrics 60, 108–115.

Sugiura, N. (1978) Further analysis of the data by Akaike's information criterion and the finite corrections. Communications in Statistics: Theory and Methods A7, 13–26.

See Also

AICc, aictabCustom, confset, evidence, c_hat, modavgCustom

Examples

##cement data from Burnham and Anderson (2002, p. 101)
data(cement)
##run multiple regression - the global model in Table 3.2
glob.mod <- lm(y ~ x1 + x2 + x3 + x4, data = cement)

##extract log-likelihood
LL <- logLik(glob.mod)[1]

##extract number of parameters
K.mod <- coef(glob.mod) + 1

##compute AICc with full likelihood
AICcCustom(LL, K.mod, nobs = nrow(cement))

---

Defunct Functions in AICcmodavg Package

Description

The functions listed below have been removed from the AICcmodavg package.

Usage

AICc.mult(...)
AICc.unmarked(...)
extract.LL(...)
extract.LL.coxph(...)
extract.LL.unmarked(...)
aictab.clm(...)
aictab.clmm(...)
aictab.coxph(...)
aictab.glm(...)
aictab.gls(...)
aictab.lm(...)
aictab.lme(...)
aictab.mer(...)
aictab.merMod(...)
aictab.mult(...)
aictab.nlme(...)
aictab.nls(...)
aictab.polr(...)
aictab.rlm(...)
aictab.unmarked(...)
dictab.bugs(...)
dictab.rjags(...)
modavg.clm(...)
modavg.clmm(...)
modavg.coxph(...)
modavg.glm(...)
modavg.gls(...)
modavg.lme(...)
modavg.mer(...)
modavg.merMod(...)
modavg.mult(...)
modavg.polr(...)
modavg.rlm(...)
modavg.unmarked(...)
modavg.effect(...)
modavg.effect.glm(...)
modavg.effect.gls(...)
modavg.effect.lme(...)
modavg.effect.mer(...)
modavg.effect.merMod(...)
modavg.effect.rlm(...)
modavg.effect.unmarked(...)
modavg.shrink(...)
modavg.shrink.clm(...)
modavg.shrink.clmm(...)
modavg.shrink.coxph(...)
modavg.shrink.glm(...)
modavg.shrink.gls(...)
modavg.shrink.lme(...)
modavg.shrink.mer(...)
modavg.shrink.merMod(...)
modavg.shrink.mult(...)
modavg.shrink.polr(...)
modavg.shrink.rlm(...)
modavg.shrink.unmarked(...)
modavgpred(...)
modavgpred.glm(...)
modavgpred.gls(...)
modavgpred.lme(...)
modavgpred.mer(...)
modavgpred.merMod(...)
modavgpred.rlm(...)
modavgpred.unmarked(...)
mult.comp(...)
predictSE.zip(...)

Arguments

...	arguments passed to the function.

Details

AICc.mult has been replaced by AICc.multinom.

AICc.unmarked has been replaced by AICc.unmarkedFit.

extract.LL has been replaced by extractLL.

extract.LL.coxph has been replaced by extractLL.coxph.

extract.LL.unmarked has been replaced by extractLL.unmarkedFit.

aictab.clm has been replaced by aictab.AICsclm.clm.

aictab.clmm has been replaced by aictab.AICclmm.

aictab.coxph has been replaced by aictab.AICcoxph.

aictab.glm has been replaced by aictab.AICglm.lm.

aictab.gls has been replaced by aictab.AICgls.

aictab.lm has been replaced by aictab.AIClm.

aictab.lme has been replaced by aictab.AIClme.

aictab.mer has been replaced by aictab.AICmer.

aictab.merMod has been replaced by aictab.AIClmerMod, aictab.AICglmerMod, or aictab.AICnlmerMod, depending on the class of the objects.

aictab.mult has been replaced by aictab.AICmultinom.nnet.

aictab.nlme has been replaced by aictab.AICnlme.

aictab.nls has been replaced by aictab.AICnls.

aictab.polr has been replaced by aictab.AICpolr.

aictab.rlm has been replaced by aictab.AICrlm.lm.

aictab.unmarked has been replaced by aictab.AICunmarkedFitOccu, aictab.AICunmarkedFitColExt, aictab.AICunmarkedFitOccuRN, aictab.AICunmarkedFitPCount, aictab.AICunmarkedFitPCO, aictab.AICunmarkedFitDS, aictab.AICunmarkedFitGDS, aictab.AICunmarkedFitOccuFP, aictab.AICunmarkedFitMPois, aictab.AICunmarkedFitGMM, or aictab.AICunmarkedFitGPC, depending on the class of the objects.

dictab.bugs has been replaced by dictab.AICbugs.

dictab.jags has been replaced by dictab.AICjags.

modavg.clm has been replaced by modavg.AICsclm.clm.

modavg.clmm has been replaced by modavg.AICsclm.clm.

modavg.coxph has been replaced by modavg.AICcoxph.

modavg.glm has been replaced by modavg.AIClm or modavg.AICglm.lm, depending on the class of the objects.

modavg.gls has been replaced by modavg.AICgls.

modavg.lme has been replaced by modavg.AIClme.

modavg.mer has been replaced by modavg.AICmer.

modavg.merMod has been replaced by modavg.AIClmerMod or modavg.AICglmerMod, depending on the class of the objects.

modavg.mult has been replaced by modavg.AICmultinom.nnet.

modavg.polr has been replaced by modavg.AICpolr.

modavg.rlm has been replaced by modavg.AICrlm.lm.

modavg.unmarked has been replaced by modavg.AICunmarkedFitOccu, modavg.AICunmarkedFitColExt, modavg.AICunmarkedFitOccuRN, modavg.AICunmarkedFitPCount, modavg.AICunmarkedFitPCO, modavg.AICunmarkedFitDS, modavg.AICunmarkedFitGDS, modavg.AICunmarkedFitOccuFP, modavg.AICunmarkedFitMPois, modavg.AICunmarkedFitGMM, or modavg.AICunmarkedFitGPC, depending on the class of the objects.

modavg.effect has been replaced by modavgEffect.

modavg.effect.glm has been replaced by modavgEffect.AICglm.lm or modavgEffect.AIClm, depending on the class of the objects.

modavg.effect.gls has been replaced by modavgEffect.AICgls.

modavg.effect.lme has been replaced by modavgEffect.AIClme.

modavg.effect.mer has been replaced by modavgEffect.AICmer.

modavg.effect.merMod has been replaced by modavgEffect.AICglmerMod or modavgEffect.AIClmerMod, depending on the class of the objects.

modavg.effect.rlm has been replaced by modavgEffect.AICrlm.lm.

modavg.effect.unmarked has been replaced by modavgEffect.AICunmarkedFitOccu, modavgEffect.AICunmarkedFitColExt, modavgEffect.AICunmarkedFitOccuRN, modavgEffect.AICunmarkedFitPCount, modavgEffect.AICunmarkedFitPCO, modavgEffect.AICunmarkedFitDS, modavgEffect.AICunmarkedFitGDS, modavgEffect.AICunmarkedFitOccuFP, modavgEffect.AICunmarkedFitMPois, modavgEffect.AICunmarkedFitGMM, or modavgEffect.AICunmarkedFitGPC, depending on the class of the objects.

modavg.shrink has been replaced by modavgShrink.

modavg.shrink.clm has been replaced by modavgShrink.AICsclm.clm.

modavg.shrink.clmm has been replaced by modavgShrink.AICclmm.

modavg.shrink.coxph has been replaced by modavgShrink.AICcoxph.

modavg.shrink.glm has been replaced by modavgShrink.AICglm.lm or modavgShrink.AICglm.lm, depending on the class of the objects.

modavg.shrink.gls has been replaced by modavgShrink.AICgls.

modavg.shrink.lme has been replaced by modavgShrink.AIClme.

modavg.shrink.mer has been replaced by modavgShrink.AICmer.

modavg.shrink.merMod has been replaced by modavgShrink.AICglmerMod or modavgShrink.AIClmerMod, depending on the class of the objects.

modavg.shrink.mult has been replaced by modavgShrink.AICmultinom.nnet.

modavg.shrink.polr has been replaced by modavgShrink.AICpolr.

modavg.shrink.rlm has been replaced by modavgShrink.AICrlm.lm

modavg.shrink.unmarked has been replaced by modavgShrink.AICunmarkedFitOccu, modavgShrink.AICunmarkedFitColExt, modavgShrink.AICunmarkedFitOccuRN, modavgShrink.AICunmarkedFitPCount, modavgShrink.AICunmarkedFitPCO, modavgShrink.AICunmarkedFitDS, modavgShrink.AICunmarkedFitGDS, modavgShrink.AICunmarkedFitOccuFP, modavgShrink.AICunmarkedFitMPois, modavgShrink.AICunmarkedFitGMM, or modavgShrink.AICunmarkedFitGPC, depending on the class of the objects.

modavgpred has been replaced by modavgPred.

modavgpred.glm has been replaced by modavgpred.AICglm.lm or modavgPred.AIClm, depending on the class of the objects.

modavgpred.gls has been replaced by modavgPred.AICgls.

modavgpred.lme has been replaced by modavgPred.AIClme.

modavgpred.mer has been replaced by modavgPred.AICmer.

modavgpred.merMod has been replaced by modavgpred.AICglmerMod or modavgPred.AIClmerMod, depending on the class of the objects.

modavgpred.rlm has been replaced by modavgPred.AICrlm.lm.

modavgpred.unmarked has been replaced by modavgPred.AICunmarkedFitOccu, modavgPred.AICunmarkedFitColExt, modavgPred.AICunmarkedFitOccuRN, modavgPred.AICunmarkedFitPCount, modavgPred.AICunmarkedFitPCO, modavgPred.AICunmarkedFitDS, modavgPred.AICunmarkedFitGDS, modavgPred.AICunmarkedFitOccuFP, modavgPred.AICunmarkedFitMPois, modavgPred.AICunmarkedFitGMM, or modavgPred.AICunmarkedFitGPC, depending on the class of the objects.

mult.comp has been replaced by multComp.

predictSE.zip has been replaced by predictSE.

Author(s)

Marc J. Mazerolle

See Also

aictab, confset, dictab, importance, evidence, extractLL, c_hat, modavg, modavgEffect, modavgShrink, modavgPred, multComp, predictSE

---

Create Model Selection Tables

Description

This function creates a model selection table based on one of the following information criteria: AIC, AICc, QAIC, QAICc. The table ranks the models based on the selected information criteria and also provides delta AIC and Akaike weights. aictab selects the appropriate function to create the model selection table based on the object class. The current version works with lists containing objects of aov, betareg, clm, clmm, clogit, coxme, coxph, fitdist, fitdistr, glm, glmmTMB, gls, gnls, hurdle, lavaan, lm, lme, lmekin, maxlikeFit, mer, merMod, lmerModLmerTest, multinom, negbin, nlme, nls, polr, rlm, survreg, vglm, and zeroinfl classes as well as various models of unmarkedFit classes but does not yet allow mixing of different classes.

Usage

aictab(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, ...)

## S3 method for class 'AICaov.lm'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICbetareg'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICsclm.clm'
aictab(cand.set, modnames = NULL,
       second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICclmm'
aictab(cand.set, modnames = NULL,
       second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICclm'
aictab(cand.set, modnames = NULL,
       second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICcoxme'
aictab(cand.set, modnames = NULL,
       second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICcoxph'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICfitdist'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICfitdistr'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICglm.lm'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)  

## S3 method for class 'AICglmmTMB'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)  

## S3 method for class 'AICgls'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICgnls.gls'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIChurdle'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIClavaan'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIClm'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIClme'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIClmekin'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICmaxlikeFit.list'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...) 

## S3 method for class 'AICmer'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIClmerMod'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIClmerModLmerTest'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICglmerMod'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...) 

## S3 method for class 'AICnlmerMod'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICmultinom.nnet'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICnegbin.glm.lm'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...) 

## S3 method for class 'AICnlme.lme'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...) 

## S3 method for class 'AICnls'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICpolr'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICrlm.lm'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICsurvreg'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICunmarkedFitOccu'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitColExt'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitOccuRN'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitPCount'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitPCO'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitDS'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitGDS'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitOccuFP'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitMPois'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitGMM'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitGPC'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitOccuMulti'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitOccuMS'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitOccuTTD'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitMMO'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitDSO'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICvglm'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICzeroinfl'
aictab(cand.set, modnames = NULL,
        second.ord = TRUE, nobs = NULL, sort = TRUE, ...)

Arguments

cand.set	a list storing each of the models in the candidate model set.
modnames	a character vector of model names to facilitate the identification of each model in the model selection table. If NULL, the function uses the names in the cand.set list of candidate models (i.e., a named list). If no names appear in the list and no character vector is provided, generic names (e.g., Mod1, Mod2) are supplied in the table in the same order as in the list of candidate models.
second.ord	logical. If TRUE, the function returns the second-order Akaike information criterion (i.e., AICc).
nobs	this argument allows to specify a numeric value other than total sample size to compute the AICc (i.e., nobs defaults to total number of observations). This is relevant only for mixed models or various models of unmarkedFit classes where sample size is not straightforward. In such cases, one might use total number of observations or number of independent clusters (e.g., sites) as the value of nobs.
sort	logical. If TRUE, the model selection table is ranked according to the (Q)AIC(c) values.
c.hat	value of overdispersion parameter (i.e., variance inflation factor) such as that obtained from c_hat. Note that values of c.hat different from 1 are only appropriate for binomial GLM's with trials > 1 (i.e., success/trial or cbind(success, failure) syntax), with Poisson GLM's, single-season occupancy models (MacKenzie et al. 2002), dynamic occupancy models (MacKenzie et al. 2003), or N-mixture models (Royle 2004, Dail and Madsen 2011). If c.hat > 1, aictab will return the quasi-likelihood analogue of the information criteria requested and multiply the variance-covariance matrix of the estimates by this value (i.e., SE's are multiplied by sqrt(c.hat)). This option is not supported for generalized linear mixed models of the mer or merMod classes.
...	additional arguments passed to the function.

Details

aictab internally creates a new class for the cand.set list of candidate models, according to the contents of the list. The current function is implemented for clogit, coxme, coxph, fitdist, fitdistr, glm, glmmTMB, gls, gnls, hurdle, lavaan, lm, lme, lmekin, maxlikeFit, mer, merMod, lmerModLmerTest, multinom, negbin, nlme, nls, polr, rlm, survreg, vglm, and zeroinfl classes as well as various unmarkedFit classes.

The function constructs a model selection table based on one of the four information criteria: AIC, AICc, QAIC, and QAICc.

Ten guidelines for model selection:

1) Carefully construct your candidate model set. Each model should represent a specific (interesting) hypothesis to test.

2) Keep your candidate model set short. It is ill-advised to consider as many models as there are data.

3) Check model fit. Use your global model (most complex model) or subglobal models to determine if the assumptions are valid. If none of your models fit the data well, information criteria will only indicate the most parsimonious of the poor models.

4) Avoid data dredging (i.e., looking for patterns after an initial round of analysis).

5) Avoid overfitting models. You should not estimate too many parameters for the number of observations available in the sample.

6) Be careful of missing values. Remember that values that are missing only for certain variables change the data set and sample size, depending on which variable is included in any given model. I suggest to remove missing cases before starting model selection.

7) Use the same response variable for all models of the candidate model set. It is inappropriate to run some models with a transformed response variable and others with the untransformed variable. A workaround is to use a different link function for some models (e.g., identity vs log link).

8) When dealing with models with overdispersion, use the same value of c-hat for all models in the candidate model set. For binomial models with trials > 1 (i.e., success/trial or cbind(success, failure) syntax) or with Poisson GLM's, you should estimate the c-hat from the most complex model (global model). If c-hat > 1, you should use the same value for each model of the candidate model set (where appropriate) and include it in the count of parameters (K). Similarly, for negative binomial models, you should estimate the dispersion parameter from the global model and use the same value across all models.

9) Burnham and Anderson (2002) recommend to avoid mixing the information-theoretic approach and notions of significance (i.e., P values). It is best to provide estimates and a measure of their precision (standard error, confidence intervals).

10) Determining the ranking of the models is just the first step. Akaike weights sum to 1 for the entire model set and can be interpreted as the weight of evidence in favor of a given model being the best one given the candidate model set considered and the data at hand. Models with large Akaike weights have strong support. Evidence ratios, importance values, and confidence sets for the best model are all measures that assist in interpretation. In cases where the top ranking model has an Akaike weight > 0.9, one can base inference on this single most parsimonious model. When many models rank highly (i.e., delta (Q)AIC(c) < 4), one should model-average effect sizes for the parameters with most support across the entire set of models. Model averaging consists in making inference based on the whole set of candidate models, instead of basing conclusions on a single 'best' model. It is an elegant way of making inference based on the information contained in the entire model set.

Value

aictab creates an object of class aictab with the following components:

Modname	the name of each model of the candidate model set.
K	the number of estimated parameters for each model.
(Q)AIC(c)	the information criterion requested for each model (AIC, AICc, QAIC, QAICc).
Delta_(Q)AIC(c)	the appropriate delta AIC component depending on the information criteria selected.
ModelLik	the relative likelihood of the model given the data (exp(-0.5*delta[i])). This is not to be confused with the likelihood of the parameters given the data. The relative likelihood can then be normalized across all models to get the model probabilities.
(Q)AIC(c)Wt	the Akaike weights, also termed "model probabilities" sensu Burnham and Anderson (2002) and Anderson (2008). These measures indicate the level of support (i.e., weight of evidence) in favor of any given model being the most parsimonious among the candidate model set.
Cum.Wt	the cumulative Akaike weights. These are only meaningful if results in table are sorted in decreasing order of Akaike weights (i.e., sort = TRUE).
c.hat	if c.hat was specified as an argument, it is included in the table.
LL	if c.hat = 1 and parameters estimated by maximum likelihood, the log-likelihood of each model.
Quasi.LL	if c.hat > 1, the quasi log-likelihood of each model.
Res.LL	if parameters are estimated by restricted maximum-likelihood (REML), the restricted log-likelihood of each model.

Author(s)

Marc J. Mazerolle

References

Anderson, D. R. (2008) Model-based Inference in the Life Sciences: a primer on evidence. Springer: New York.

Burnham, K. P., Anderson, D. R. (2002) Model Selection and Multimodel Inference: a practical information-theoretic approach. Second edition. Springer: New York.

Burnham, K. P., Anderson, D. R. (2004) Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods and Research 33, 261–304.

Dail, D., Madsen, L. (2011) Models for estimating abundance from repeated counts of an open population. Biometrics 67, 577–587.

MacKenzie, D. I., Nichols, J. D., Lachman, G. B., Droege, S., Royle, J. A., Langtimm, C. A. (2002) Estimating site occupancy rates when detection probabilities are less than one. Ecology 83, 2248–2255.

MacKenzie, D. I., Nichols, J. D., Hines, J. E., Knutson, M. G., Franklin, A. B. (2003) Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84, 2200–2207.

Mazerolle, M. J. (2006) Improving data analysis in herpetology: using Akaike's Information Criterion (AIC) to assess the strength of biological hypotheses. Amphibia-Reptilia 27, 169–180.

Royle, J. A. (2004) N-mixture models for estimating population size from spatially replicated counts. Biometrics 60, 108–115.

See Also

AICc, aictabCustom, bictab, confset, c_hat, evidence, importance, modavg, modavgEffect, modavgShrink, modavgPred

Examples

##Mazerolle (2006) frog water loss example
data(dry.frog)

##setup a subset of models of Table 1
Cand.models <- list( )
Cand.models[[1]] <- lm(log_Mass_lost ~ Shade + Substrate +
                       cent_Initial_mass + Initial_mass2,
                       data = dry.frog)
Cand.models[[2]] <- lm(log_Mass_lost ~ Shade + Substrate +
                       cent_Initial_mass + Initial_mass2 +
                       Shade:Substrate, data = dry.frog)
Cand.models[[3]] <- lm(log_Mass_lost ~ cent_Initial_mass +
                       Initial_mass2, data = dry.frog)
Cand.models[[4]] <- lm(log_Mass_lost ~ Shade + cent_Initial_mass +
                       Initial_mass2, data = dry.frog)
Cand.models[[5]] <- lm(log_Mass_lost ~ Substrate + cent_Initial_mass +
                       Initial_mass2, data = dry.frog)

##create a vector of names to trace back models in set
Modnames <- paste("mod", 1:length(Cand.models), sep = " ")

##generate AICc table
aictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE)
##round to 4 digits after decimal point and give log-likelihood
print(aictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE),
      digits = 4, LL = TRUE)



## Not run: 
##Burnham and Anderson (2002) flour beetle data
data(beetle)
##models as suggested by Burnham and Anderson p. 198          
Cand.set <- list( )
Cand.set[[1]] <- glm(Mortality_rate ~ Dose, family =
                     binomial(link = "logit"), weights = Number_tested,
                     data = beetle)
Cand.set[[2]] <- glm(Mortality_rate ~ Dose, family =
                     binomial(link = "probit"), weights = Number_tested,
                     data = beetle)
Cand.set[[3]] <- glm(Mortality_rate ~ Dose, family =
                     binomial(link ="cloglog"), weights = Number_tested,
                     data = beetle)

##check c-hat
c_hat(Cand.set[[1]])
c_hat(Cand.set[[2]])
c_hat(Cand.set[[3]])
##lowest value of c-hat < 1 for these non-nested models, thus use
##c.hat = 1 
       
##set up named list
names(Cand.set) <- c("logit", "probit", "cloglog")

##compare models
##model names will be taken from the list if modnames is not specified
res.table <- aictab(cand.set = Cand.set, second.ord = FALSE)
##note that delta AIC and Akaike weights are identical to Table 4.7
print(res.table, digits = 2, LL = TRUE) #print table with 2 digits and
##print log-likelihood in table
print(res.table, digits = 4, LL = FALSE) #print table with 4 digits and
##do not print log-likelihood

## End(Not run)


##two-way ANOVA with interaction
data(iron)
##full model
m1 <- lm(Iron ~ Pot + Food + Pot:Food, data = iron)
##additive model
m2 <- lm(Iron ~ Pot + Food, data = iron)
##null model
m3 <- lm(Iron ~ 1, data = iron)

##candidate models
Cand.aov <- list(m1, m2, m3)
Cand.names <- c("full", "additive", "null")
aictab(Cand.aov, Cand.names)



##single-season occupancy model example modified from ?occu
## Not run: 
require(unmarked)
##single season example modified from ?occu
data(frogs)
pferUMF <- unmarkedFrameOccu(pfer.bin)
##add fake covariates
siteCovs(pferUMF) <- data.frame(sitevar1 = rnorm(numSites(pferUMF)),
                                sitevar2 = runif(numSites(pferUMF))) 
     
##observation covariates
obsCovs(pferUMF) <- data.frame(obsvar1 = rnorm(numSites(pferUMF) *
                                 obsNum(pferUMF))) 

##set up candidate model set
fm1 <- occu(~ obsvar1 ~ sitevar1, pferUMF)
fm2 <- occu(~ 1 ~ sitevar1, pferUMF)
fm3 <- occu(~ obsvar1 ~ sitevar2, pferUMF)
fm4 <- occu(~ 1 ~ sitevar2, pferUMF)

##assemble models in named list (alternative to using 'modnames' argument)
Cand.mods <- list("fm1" = fm1, "fm2" = fm2, "fm3" = fm3, "fm4" = fm4)

##compute table
aictab(cand.set = Cand.mods, second.ord = TRUE)

detach(package:unmarked)

## End(Not run)

---

Create Model Selection Tables from User-supplied Input Based on (Q)AIC(c)

Description

This function creates a model selection table from model input (log-likelihood, number of estimated parameters) supplied by the user instead of extracting the values automatically from a list of candidate models. The models are ranked based on one of the following information criteria: AIC, AICc, QAIC, QAICc. The table ranks the models based on the selected information criteria and also provides delta AIC and Akaike weights.

Usage

aictabCustom(logL, K, modnames = NULL, second.ord = TRUE, nobs = NULL,
             sort = TRUE, c.hat = 1)

Arguments

logL	a vector of log-likelihood values for the models in the candidate model set.
K	a vector containing the number of estimated parameters for each model in the candidate model set.
modnames	a character vector of model names to facilitate the identification of each model in the model selection table. If NULL, the function uses the names in the cand.set list of candidate models (i.e., a named list). If no names appear in the list and no character vector is provided, generic names (e.g., Mod1, Mod2) are supplied in the table in the same order as in the list of candidate models.
second.ord	logical. If TRUE, the function returns the second-order Akaike information criterion (i.e., AICc).
nobs	the sample size required to compute the AICc or QAICc.
sort	logical. If TRUE, the model selection table is ranked according to the (Q)AIC(c) values.
c.hat	value of overdispersion parameter (i.e., variance inflation factor) such as that obtained from c_hat. Note that values of c.hat different from 1 are only appropriate for binomial GLM's with trials > 1 (i.e., success/trial or cbind(success, failure) syntax), with Poisson GLM's, single-season or dynamic occupancy models (MacKenzie et al. 2002, 2003), N-mixture models (Royle 2004, Dail and Madsen 2011), or capture-mark-recapture models (e.g., Lebreton et al. 1992). If c.hat > 1, aictabCustom will return the quasi-likelihood analogue of the information criterion requested.

Details

aictabCustom constructs a model selection table based on one of the four information criteria: AIC, AICc, QAIC, and QAICc. This function is most useful when model input is imported into R from other software (e.g., Program MARK, PRESENCE) or for model classes that are not yet supported by aictab.

Value

aictabCustom creates an object of class aictab with the following components:

Modname	the name of each model of the candidate model set.
K	the number of estimated parameters for each model.
(Q)AIC(c)	the information criteria requested for each model (AICc, AICc, QAIC, QAICc).
Delta_(Q)AIC(c)	the appropriate delta AIC component depending on the information criteria selected.
ModelLik	the relative likelihood of the model given the data (exp(-0.5*delta[i])). This is not to be confused with the likelihood of the parameters given the data. The relative likelihood can then be normalized across all models to get the model probabilities.
(Q)AIC(c)Wt	the Akaike weights, also termed "model probabilities" sensu Burnham and Anderson (2002) and Anderson (2008). These measures indicate the level of support (i.e., weight of evidence) in favor of any given model being the most parsimonious among the candidate model set.
Cum.Wt	the cumulative Akaike weights. These are only meaningful if results in table are sorted in decreasing order of Akaike weights (i.e., sort = TRUE).
c.hat	if c.hat was specified as an argument, it is included in the table.
LL	if c.hat = 1 and parameters estimated by maximum likelihood, the log-likelihood of each model.
Quasi.LL	if c.hat > 1, the quasi log-likelihood of each model.

Author(s)

Marc J. Mazerolle

References

Anderson, D. R. (2008) Model-based Inference in the Life Sciences: a primer on evidence. Springer: New York.

Burnham, K. P., Anderson, D. R. (2002) Model Selection and Multimodel Inference: a practical information-theoretic approach. Second edition. Springer: New York.

Dail, D., Madsen, L. (2011) Models for estimating abundance from repeated counts of an open population. Biometrics 67, 577–587.

Lebreton, J.-D., Burnham, K. P., Clobert, J., Anderson, D. R. (1992) Modeling survival and testing biological hypotheses using marked animals: a unified approach with case-studies. Ecological Monographs 62, 67–118.

MacKenzie, D. I., Nichols, J. D., Lachman, G. B., Droege, S., Royle, J. A., Langtimm, C. A. (2002) Estimating site occupancy rates when detection probabilities are less than one. Ecology 83, 2248–2255.

MacKenzie, D. I., Nichols, J. D., Hines, J. E., Knutson, M. G., Franklin, A. B. (2003) Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84, 2200–2207.

Mazerolle, M. J. (2006) Improving data analysis in herpetology: using Akaike's Information Criterion (AIC) to assess the strength of biological hypotheses. Amphibia-Reptilia 27, 169–180.

Royle, J. A. (2004) N-mixture models for estimating population size from spatially replicated counts. Biometrics 60, 108–115.

See Also

AICcCustom, bictabCustom, confset, c_hat, evidence, ictab, modavgCustom

Examples

##vector with model LL's
LL <- c(-38.8876, -35.1783, -64.8970)

##vector with number of parameters
Ks <- c(7, 9, 4)

##create a vector of names to trace back models in set
Modnames <- c("Cm1", "Cm2", "Cm3")

##generate AICc table
aictabCustom(logL = LL, K = Ks, modnames = Modnames, nobs = 121,
             sort = TRUE)

---

Likelihood-Ratio Test Corrected for Overdispersion

Description

Compute likelihood-ratio test between a given model and a simpler model.

Usage

anovaOD(mod.simple, mod.complex, c.hat = 1, 
        nobs = NULL, ...)

## S3 method for class 'glm'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitOccu'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitColExt'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitOccuRN'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitPCount'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitPCO'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitDS'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitGDS'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitOccuFP'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitMPois'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitGMM'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitGPC'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitOccuMS'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitOccuTTD'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitMMO'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitDSO'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'glmerMod'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'maxlikeFit'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'multinom'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'vglm'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

Arguments

mod.simple	an object of class glm, glmmTMB, maxlikeFit, mer, merMod, multinom, vglm, and various unmarkedFit classes containing the output of a model. This model should be a simpler version of mod.complex resulting from a deletion of certain terms (i.e., nested model).
mod.complex	an object of the same class as mod.simple.
c.hat	value of overdispersion parameter (i.e., variance inflation factor) such as that obtained from c_hat, mb.gof.test, or Nmix.gof.test. Typically, this value should be computed for the most complex model and applied to simpler models.
nobs	the number of observations used in the analysis. If nobs = NULL, the total number of rows are used as the sample size to compute the residual degrees of freedom as $nobs - K$, where K is the number of estimated parameters. This is relevant only for mixed models or various models of unmarkedFit classes where sample size is not straightforward. In such cases, one might use total number of observations or number of independent clusters (e.g., sites) as the value of nobs.
...	additional arguments passed to the function.

Details

This function applies a correction for overdispersion on the likelihood-ratio test between a model and its simpler counterpart. The simpler model must be nested within the more complex model, typically as the result of deleting terms. You should supply the c.hat value of the most complex of the two models you are comparing.

When $1 < c.hat < 4$, the likelihood-ratio test is computed as:

$LR = \frac{-2 * (LL.simple - LL.complex)}{(K.complex - K.simple) * c.hat}$

where LL.simple and LL.complex are the log-likelihoods of the simple and complex models, respectively, and where K.complex and K.simple are the number of estimated parameters in each model. The test statistic is approximately distributed as $F_{K.complex - K.simple, n - K.complex}$, where n is the number of observations (i.e., nobs) used in the analysis (Venables and Ripley 2002).

When nobs = NULL, the number of observations is based on the number of rows of the data frame used in the analysis. For mixed models or various models of unmarkedFit, sample size is less straightforward, and nobs could be based on the total number of observations or on the number of independent clusters (e.g., sites), among other choices.

When c.hat = 1, the likelihood-ratio test simplifies to:

$LR = -2 * (LL.simple - LL.complex)$

where in this case the test statistic is distributed as a $\chi^2_{K.complex - K.simple}$ (McCullagh and Nelder 1989).

The function supports different model types such as Poisson GLM's and GLMM's, single-season and dynamic occupancy models (MacKenzie et al. 2002, 2003), and various N-mixture models (Royle 2004, Dail and Madsen 2011).

Value

anovaOD returns an object of class anovaOD as a list with the following components:

form.simple	a character string of the parameters estimated in mod.simple.
form.complex	a character string of the parameters estimated in mod.complex.
c.hat	the c.hat estimate used to adjust the likelihood-ratio test.
devMat	a matrix storing as columns the number of parameters estimated (K), the log-likelihood of each model logLik, the difference in estimated parameters between the two models (Kdiff), minus twice the difference in log-likelihoods between the models (-2LL), the test statistic, and the associated P-value.

Author(s)

Marc J. Mazerolle

References

Dail, D., Madsen, L. (2011) Models for estimating abundance from repeated counts of an open population. Biometrics 67, 577–587.

MacKenzie, D. I., Nichols, J. D., Lachman, G. B., Droege, S., Royle, J. A., Langtimm, C. A. (2002) Estimating site occupancy rates when detection probabilities are less than one. Ecology 83, 2248–2255.

MacKenzie, D. I., Nichols, J. D., Hines, J. E., Knutson, M. G., Franklin, A. B. (2003) Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84, 2200–2207.

Mazerolle, M. J. (2006) Improving data analysis in herpetology: using Akaike's Information Criterion (AIC) to assess the strength of biological hypotheses. Amphibia-Reptilia 27, 169–180.

McCullagh, P., Nelder, J. A. (1989) Generalized Linear Models. Second edition. Chapman and Hall: New York.

Royle, J. A. (2004) N-mixture models for estimating population size from spatially replicated counts. Biometrics 60, 108–115.

Venables, W. N., Ripley, B. D. (2002) Modern Applied Statistics with S. Second edition. Springer-Verlag: New York.

See Also

c_hat, mb.gof.test, Nmix.gof.test, summaryOD

Examples

##anuran larvae example from Mazerolle (2006)
data(min.trap)
##assign "UPLAND" as the reference level as in Mazerolle (2006)          
min.trap$Type <- relevel(min.trap$Type, ref = "UPLAND") 

##run model
m1 <- glm(Num_anura ~ Type + log.Perimeter + Num_ranatra,
          family = poisson, offset = log(Effort),
          data = min.trap) 
##null model
m0 <- glm(Num_anura ~ 1,
          family = poisson, offset = log(Effort),
          data = min.trap) 

##check c-hat for global model
c_hat(m1) #uses Pearson's chi-square/df

##likelihood ratio test corrected for overdispersion
anovaOD(mod.simple = m0, mod.complex = m1, c.hat = c_hat(m1))
##compare without overdispersion correction
anovaOD(mod.simple = m0, mod.complex = m1)


##example with occupancy model
## Not run: 
##load unmarked package
if(require(unmarked)){
   
   data(bullfrog)
     
   ##detection data
   detections <- bullfrog[, 3:9]

   ##assemble in unmarkedFrameOccu
   bfrog <- unmarkedFrameOccu(y = detections)
     
   ##run model
   fm <- occu(~ 1 ~ Reed.presence, data = bfrog)
   ##null model
   fm0 <- occu(~ 1 ~ 1, data = bfrog)

   ##check GOF
   ##GOF <- mb.gof.test(fm, nsim = 1000)
   ##estimate of c-hat:  1.89

   ##display results after overdispersion adjustment
   anovaOD(fm0, fm, c.hat = 1.89)

   detach(package:unmarked)
}

## End(Not run)

---

Flour Beetle Data

Description

This data set illustrates the acute mortality of flour beetles (Tribolium confusum) following 5 hour exposure to carbon disulfide gas.

Usage

data(beetle)

Format

A data frame with 8 rows and 4 variables.

Dose

dose of carbon disulfide in mg/L.

Number_tested

number of beetles exposed to given dose of carbon disulfide.

Number_killed

number of beetles dead after 5 hour exposure to given dose of carbon disulfide.

Mortality_rate

proportion of total beetles found dead after 5 hour exposure.

Details

Burnham and Anderson (2002, p. 195) use this data set originally from Young and Young (1998) to show model selection for binomial models with different link functions (logit, probit, cloglog).

Source

Burnham, K. P., Anderson, D. R. (2002) Model Selection and Multimodel Inference: a practical information-theoretic approach. Second edition. Springer: New York.

Young, L. J., Young, J. H. (1998) Statistical Ecology. Kluwer Academic Publishers: London.

Examples

data(beetle)
## maybe str(beetle) ; plot(beetle) ...

---

Create Model Selection Tables Based on BIC

Description

This function creates a model selection table based on the Bayesian information criterion (Schwarz 1978, Burnham and Anderson 2002). The table ranks the models based on the BIC and also provides delta BIC and BIC model weights. The function adjusts for overdispersion in model selection by using the QBIC when c.hat > 1. bictab selects the appropriate function to create the model selection table based on the object class. The current version works with lists containing objects of aov, betareg, clm, clmm, clogit, coxme, coxph, fitdist, fitdistr, glm, glmmTMB, gls, gnls, hurdle, lavaan, lm, lme, lmekin, maxlikeFit, mer, merMod, lmerModLmerTest, multinom, nlme, nls, polr, rlm, survreg, vglm, and zeroinfl classes as well as various models of unmarkedFit classes but does not yet allow mixing of different classes.

Usage

bictab(cand.set, modnames = NULL, nobs = NULL,
       sort = TRUE, ...)

## S3 method for class 'AICaov.lm'
bictab(cand.set, modnames = NULL,
        nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICbetareg'
bictab(cand.set, modnames = NULL,
        nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICsclm.clm'
bictab(cand.set, modnames = NULL,
        nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICclmm'
bictab(cand.set, modnames = NULL,
        nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICclm'
bictab(cand.set, modnames = NULL,
        nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICcoxme'
bictab(cand.set, modnames = NULL,
        nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICcoxph'
bictab(cand.set, modnames = NULL,
        nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICfitdist'
bictab(cand.set, modnames = NULL,
        nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICfitdistr'
bictab(cand.set, modnames = NULL,
        nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICglm.lm'
bictab(cand.set, modnames = NULL,
        nobs = NULL, sort = TRUE, c.hat = 1, ...)  

## S3 method for class 'AICglmmTMB'
bictab(cand.set, modnames = NULL,
        nobs = NULL, sort = TRUE, c.hat = 1, ...)  

## S3 method for class 'AICgls'
bictab(cand.set, modnames = NULL,
        nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICgnls.gls'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIChurdle'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIClavaan'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIClm'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIClme'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIClmekin'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICmaxlikeFit.list'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...) 

## S3 method for class 'AICmer'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIClmerMod'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AIClmerModLmerTest'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICglmerMod'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...) 

## S3 method for class 'AICnlmerMod'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICmultinom.nnet'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICnlme.lme'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...) 

## S3 method for class 'AICnls'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICpolr'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICrlm.lm'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICsurvreg'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

## S3 method for class 'AICunmarkedFitOccu'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitColExt'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitOccuRN'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitPCount'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitPCO'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitDS'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitGDS'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitOccuFP'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitMPois'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitGMM'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitGPC'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitOccuMulti'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitOccuMS'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitOccuTTD'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitMMO'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICunmarkedFitDSO'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICvglm'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, c.hat = 1, ...)

## S3 method for class 'AICzeroinfl'
bictab(cand.set, modnames = NULL,
         nobs = NULL, sort = TRUE, ...)

Arguments

cand.set	a list storing each of the models in the candidate model set.
modnames	a character vector of model names to facilitate the identification of each model in the model selection table. If NULL, the function uses the names in the cand.set list of candidate models (i.e., a named list). If no names appear in the list and no character vector is provided, generic names (e.g., Mod1, Mod2) are supplied in the table in the same order as in the list of candidate models.
nobs	this argument allows to specify a numeric value other than total sample size to compute the BIC (i.e., nobs defaults to total number of observations). This is relevant only for mixed models or various models of unmarkedFit classes where sample size is not straightforward. In such cases, one might use total number of observations or number of independent clusters (e.g., sites) as the value of nobs.
sort	logical. If TRUE, the model selection table is ranked according to the BIC values.
c.hat	value of overdispersion parameter (i.e., variance inflation factor) such as that obtained from c_hat. Note that values of c.hat different from 1 are only appropriate for binomial GLM's with trials > 1 (i.e., success/trial or cbind(success, failure) syntax), with Poisson GLM's, single-season occupancy models (MacKenzie et al. 2002), dynamic occupancy models (MacKenzie et al. 2003), or N-mixture models (Royle 2004, Dail and Madsen 2011). If c.hat > 1, bictab will return the quasi-likelihood analogue of the BIC (QBIC) and multiply the variance-covariance matrix of the estimates by this value (i.e., SE's are multiplied by sqrt(c.hat)). This option is not supported for generalized linear mixed models of the mer or merMod classes.
...	additional arguments passed to the function.

Details

BIC tends to favor simpler models than AIC whenever n > 8 (Schwarz 1978, Link and Barker 2006, Anderson 2008). BIC assigns uniform prior probabilities across all models (i.e., equal 1/R), whereas in AIC and AICc, prior probabilities increase with sample size (Burnham and Anderson 2004, Link and Barker 2010). Some authors argue that BIC requires the true model to be included in the model set, whereas AIC or AICc does not (Burnham and Anderson 2002). However, Link and Barker (2006, 2010) consider both as assuming that a model in the model set approximates truth.

bictab internally creates a new class for the cand.set list of candidate models, according to the contents of the list. The current function is implemented for clogit, coxme, coxph, fitdist, fitdistr, glm, glmmTMB, gls, gnls, hurdle, lavaan, lm, lme, lmekin, maxlikeFit, mer, merMod, lmerModLmerTest, multinom, nlme, nls, polr, rlm, survreg, vglm, and zeroinfl classes as well as various unmarkedFit classes. The function constructs a model selection table based on BIC.

Value

bictab creates an object of class bictab with the following components:

Modname	the name of each model of the candidate model set.
K	the number of estimated parameters for each model.
(Q)BIC	the Bayesian information criterion for each model.
Delta_(Q)BIC	the delta BIC component.
ModelLik	the relative likelihood of the model given the data (exp(-0.5*delta[i])). This is not to be confused with the likelihood of the parameters given the data. The relative likelihood can then be normalized across all models to get the model probabilities.
(Q)BICWt	the BIC model weights, also termed "model probabilities" (Burnham and Anderson 2002, Link and Barker 2006, Anderson 2008). These measures indicate the level of support (i.e., weight of evidence) in favor of any given model being the most parsimonious among the candidate model set.
Cum.Wt	the cumulative BIC weights. These are only meaningful if results in table are sorted in decreasing order of BIC weights (i.e., sort = TRUE).
c.hat	if c.hat was specified as an argument, it is included in the table.
LL	the log-likelihood of each model.
Res.LL	if parameters are estimated by restricted maximum-likelihood (REML), the restricted log-likelihood of each model.

Author(s)

Marc J. Mazerolle

References

Anderson, D. R. (2008) Model-based Inference in the Life Sciences: a primer on evidence. Springer: New York.

Burnham, K. P., Anderson, D. R. (2002) Model Selection and Multimodel Inference: a practical information-theoretic approach. Second edition. Springer: New York.

Burnham, K. P., Anderson, D. R. (2004) Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods and Research 33, 261–304.

Dail, D., Madsen, L. (2011) Models for estimating abundance from repeated counts of an open population. Biometrics 67, 577–587.

Link, W. A., Barker, R. J. (2006) Model weights and the foundations of multimodel inference. Ecology 87, 2626–2635.

Link, W. A., Barker, R. J. (2010) Bayesian Inference with Ecological Applications. Academic Press: Boston.

MacKenzie, D. I., Nichols, J. D., Lachman, G. B., Droege, S., Royle, J. A., Langtimm, C. A. (2002) Estimating site occupancy rates when detection probabilities are less than one. Ecology 83, 2248–2255.

MacKenzie, D. I., Nichols, J. D., Hines, J. E., Knutson, M. G., Franklin, A. B. (2003) Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84, 2200–2207.

Royle, J. A. (2004) N-mixture models for estimating population size from spatially replicated counts. Biometrics 60, 108–115.

Schwarz, G. (1978) Estimating the dimension of a model. Annals of Statistics 6, 461–464.

See Also

aictab, bictabCustom, confset, evidence, importance, useBIC,

Examples

##Mazerolle (2006) frog water loss example
data(dry.frog)

##setup a subset of models of Table 1
Cand.models <- list( )
Cand.models[[1]] <- lm(log_Mass_lost ~ Shade + Substrate +
                       cent_Initial_mass + Initial_mass2,
                       data = dry.frog)
Cand.models[[2]] <- lm(log_Mass_lost ~ Shade + Substrate +
                       cent_Initial_mass + Initial_mass2 +
                       Shade:Substrate, data = dry.frog)
Cand.models[[3]] <- lm(log_Mass_lost ~ cent_Initial_mass +
                       Initial_mass2, data = dry.frog)
Cand.models[[4]] <- lm(log_Mass_lost ~ Shade + cent_Initial_mass +
                       Initial_mass2, data = dry.frog)
Cand.models[[5]] <- lm(log_Mass_lost ~ Substrate + cent_Initial_mass +
                       Initial_mass2, data = dry.frog)

##create a vector of names to trace back models in set
Modnames <- paste("mod", 1:length(Cand.models), sep = " ")

##generate BIC table
bictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE)
##round to 4 digits after decimal point and give log-likelihood
print(bictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE),
      digits = 4, LL = TRUE)



## Not run: 
##Burnham and Anderson (2002) flour beetle data
data(beetle)
##models as suggested by Burnham and Anderson p. 198          
Cand.set <- list( )
Cand.set[[1]] <- glm(Mortality_rate ~ Dose, family =
                     binomial(link = "logit"), weights = Number_tested,
                     data = beetle)
Cand.set[[2]] <- glm(Mortality_rate ~ Dose, family =
                     binomial(link = "probit"), weights = Number_tested,
                     data = beetle)
Cand.set[[3]] <- glm(Mortality_rate ~ Dose, family =
                     binomial(link ="cloglog"), weights = Number_tested,
                     data = beetle)

##set up named list
names(Cand.set) <- c("logit", "probit", "cloglog")

##compare models
##model names will be taken from the list if modnames is not specified
bictab(cand.set = Cand.set)

## End(Not run)


##two-way ANOVA with interaction
data(iron)
##full model
m1 <- lm(Iron ~ Pot + Food + Pot:Food, data = iron)
##additive model
m2 <- lm(Iron ~ Pot + Food, data = iron)
##null model
m3 <- lm(Iron ~ 1, data = iron)

##candidate models
Cand.aov <- list(m1, m2, m3)
Cand.names <- c("full", "additive", "null")
bictab(Cand.aov, Cand.names)



##single-season occupancy model example modified from ?occu
## Not run: 
require(unmarked)
##single season example modified from ?occu
data(frogs)
pferUMF <- unmarkedFrameOccu(pfer.bin)
##add fake covariates
siteCovs(pferUMF) <- data.frame(sitevar1 = rnorm(numSites(pferUMF)),
                                sitevar2 = runif(numSites(pferUMF))) 
     
##observation covariates
obsCovs(pferUMF) <- data.frame(obsvar1 = rnorm(numSites(pferUMF) *
                                 obsNum(pferUMF))) 

##set up candidate model set
fm1 <- occu(~ obsvar1 ~ sitevar1, pferUMF)
fm2 <- occu(~ 1 ~ sitevar1, pferUMF)
fm3 <- occu(~ obsvar1 ~ sitevar2, pferUMF)
fm4 <- occu(~ 1 ~ sitevar2, pferUMF)

##assemble models in named list (alternative to using 'modnames' argument)
Cand.mods <- list("fm1" = fm1, "fm2" = fm2, "fm3" = fm3, "fm4" = fm4)

##compute table based on QBIC that accounts for c.hat
bictab(cand.set = Cand.mods, c.hat = 3.9)

detach(package:unmarked)

## End(Not run)

---

Create Model Selection Tables from User-supplied Input Based on (Q)BIC

Description

This function creates a model selection table from model input (log-likelihood, number of estimated parameters) supplied by the user instead of extracting the values automatically from a list of candidate models. The models are ranked based on the BIC (Schwarz 1978) or on a quasi-likelihood analogue (QBIC) corrected for overdispersion. The table ranks the models based on the selected information criteria and also provides delta BIC and BIC weights.

Usage

bictabCustom(logL, K, modnames = NULL, nobs = NULL, sort = TRUE,
             c.hat = 1)

Arguments

logL	a vector of log-likelihood values for the models in the candidate model set.
K	a vector containing the number of estimated parameters for each model in the candidate model set.
modnames	a character vector of model names to facilitate the identification of each model in the model selection table. If NULL, the function uses the names in the cand.set list of candidate models (i.e., a named list). If no names appear in the list and no character vector is provided, generic names (e.g., Mod1, Mod2) are supplied in the table in the same order as in the list of candidate models.
nobs	the sample size required to compute the AICc or QAICc.
sort	logical. If TRUE, the model selection table is ranked according to the (Q)BIC values.
c.hat	value of overdispersion parameter (i.e., variance inflation factor) such as that obtained from c_hat. Note that values of c.hat different from 1 are only appropriate for binomial GLM's with trials > 1 (i.e., success/trial or cbind(success, failure) syntax), with Poisson GLM's, single-season or dynamic occupancy models (MacKenzie et al. 2002, 2003), N-mixture models (Royle 2004, Dail and Madsen 2011), or capture-mark-recapture models (e.g., Lebreton et al. 1992). If c.hat > 1, bictabCustom will return the quasi-likelihood analogue of the information criterion requested.

Details

bictabCustom constructs a model selection table based on BIC or QBIC. This function is most useful when model input is imported into R from other software (e.g., Program MARK, PRESENCE) or for model classes that are not yet supported by bictab.

Value

bictabCustom creates an object of class bictab with the following components:

Modname	the name of each model of the candidate model set.
K	the number of estimated parameters for each model.
(Q)BIC	the information criteria requested for each model (BIC, QBIC).
Delta_(Q)BIC	the appropriate delta BIC component depending on the information criteria selected.
ModelLik	the relative likelihood of the model given the data (exp(-0.5*delta[i])). This is not to be confused with the likelihood of the parameters given the data. The relative likelihood can then be normalized across all models to get the model probabilities.
(Q)BICWt	the BIC weights, also termed "model probabilities" sensu Burnham and Anderson (2002) and Anderson (2008). These measures indicate the level of support (i.e., weight of evidence) in favor of any given model being the most parsimonious among the candidate model set.
Cum.Wt	the cumulative BIC weights. These are only meaningful if results in table are sorted in decreasing order of BIC weights (i.e., sort = TRUE).
c.hat	if c.hat was specified as an argument, it is included in the table.
LL	if c.hat = 1 and parameters estimated by maximum likelihood, the log-likelihood of each model.
Quasi.LL	if c.hat > 1, the quasi log-likelihood of each model.

Author(s)

Marc J. Mazerolle

References

Anderson, D. R. (2008) Model-based Inference in the Life Sciences: a primer on evidence. Springer: New York.

Burnham, K. P., Anderson, D. R. (2002) Model Selection and Multimodel Inference: a practical information-theoretic approach. Second edition. Springer: New York.

Dail, D., Madsen, L. (2011) Models for estimating abundance from repeated counts of an open population. Biometrics 67, 577–587.

Lebreton, J.-D., Burnham, K. P., Clobert, J., Anderson, D. R. (1992) Modeling survival and testing biological hypotheses using marked animals: a unified approach with case-studies. Ecological Monographs 62, 67–118.

MacKenzie, D. I., Nichols, J. D., Lachman, G. B., Droege, S., Royle, J. A., Langtimm, C. A. (2002) Estimating site occupancy rates when detection probabilities are less than one. Ecology 83, 2248–2255.

MacKenzie, D. I., Nichols, J. D., Hines, J. E., Knutson, M. G., Franklin, A. B. (2003) Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84, 2200–2207.

Mazerolle, M. J. (2006) Improving data analysis in herpetology: using Akaike's Information Criterion (AIC) to assess the strength of biological hypotheses. Amphibia-Reptilia 27, 169–180.

Royle, J. A. (2004) N-mixture models for estimating population size from spatially replicated counts. Biometrics 60, 108–115.

Schwarz, G. (1978) Estimating the dimension of a model. Annals of Statistics 6, 461–464.

See Also

AICcCustom, aictabCustom, confset, c_hat, evidence, ictab, modavgCustom

Examples

##vector with model LL's
LL <- c(-38.8876, -35.1783, -64.8970)

##vector with number of parameters
Ks <- c(7, 9, 4)

##create a vector of names to trace back models in set
Modnames <- c("Cm1", "Cm2", "Cm3")

##generate BIC table
bictabCustom(logL = LL, K = Ks, modnames = Modnames, nobs = 121,
             sort = TRUE)

---

Compute Model Selection Relative Frequencies

Description

This function computes the model selection relative frequencies based on the nonparametric bootstrap (Burnham and Anderson 2002). Models are ranked based on the AIC, AICc, QAIC, or QAICc. The function currently supports objects of aov, betareg, clm, glm, hurdle, lm, multinom, polr, rlm, survreg, vglm, and zeroinfl classes.

Usage

boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
        sort = TRUE, nsim = 100, ...)

## S3 method for class 'AICaov.lm'
boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, nsim = 100, ...)

## S3 method for class 'AICsurvreg'
boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, nsim = 100, ...)

## S3 method for class 'AICsclm.clm'
boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, nsim = 100, ...)

## S3 method for class 'AICglm.lm'
boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, nsim = 100, c.hat = 1, ...)

## S3 method for class 'AIChurdle'
boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, nsim = 100, ...)
       
## S3 method for class 'AIClm'
boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, nsim = 100, ...)

## S3 method for class 'AICmultinom.nnet'
boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, nsim = 100, c.hat = 1, ...)

## S3 method for class 'AICpolr'
boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, nsim = 100, ...)

## S3 method for class 'AICrlm.lm'
boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, nsim = 100, ...)

## S3 method for class 'AICsurvreg'
boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, nsim = 100, ...)

## S3 method for class 'AICvglm'
boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, nsim = 100, c.hat = 1, ...)

## S3 method for class 'AICzeroinfl'
boot.wt(cand.set, modnames = NULL, second.ord = TRUE, nobs = NULL,
       sort = TRUE, nsim = 100, ...)

Arguments

cand.set	a list storing each of the models in the candidate model set.
modnames	a character vector of model names to facilitate the identification of each model in the model selection table. If NULL, the function uses the names in the cand.set list of candidate models. If no names appear in the list, generic names (e.g., Mod1, Mod2) are supplied in the table in the same order as in the list of candidate models.
second.ord	logical. If TRUE, the function returns the second-order Akaike information criterion (i.e., AICc).
nobs	this argument allows to specify a numeric value other than total sample size to compute the AICc (i.e., nobs defaults to total number of observations). This is relevant only for certain types of models such as mixed models where sample size is not straightforward. In such cases, one might use total number of observations or number of independent clusters (e.g., sites) as the value of nobs.
sort	logical. If TRUE, the model selection table is ranked according to the (Q)AIC(c) values.
c.hat	value of overdispersion parameter (i.e., variance inflation factor) such as that obtained from c_hat. Note that values of c.hat different from 1 are only appropriate for binomial GLM's with trials > 1 (i.e., success/trial or cbind(success, failure) syntax) or with Poisson GLM's. If c.hat > 1, boot.wt will return the quasi-likelihood analogue of the information criterion requested.
nsim	the number of bootstrap iterations. Burnham and Anderson (2002) recommend at least 1000 and up to 10 000 iterations for certain problems.
...	additional arguments passed to the function.

Details

boot.wt is implemented for aov, betareg, glm, hurdle, lm, multinom, polr, rlm, survreg, vglm, and zeroinfl classes. During each bootstrap iteration, the data are resampled with replacement, all the models specified in cand.set are updated with the new data set, and the top-ranked model is saved. When all iterations are completed, the relative frequency of selection is computed for each model appearing in the candidate model set.

Relative frequencies of the models are often similar to Akaike weights, and the latter are often preferred due to their link with a Bayesian perspective (Burnham and Anderson 2002). boot.wt is most useful for teaching purposes of sampling-theory based relative frequencies of model selection. The current implementation is only appropriate with completely randomized designs. For more complex data structures (e.g., blocks or random effects), the bootstrap should be modified accordingly.

Value

boot.wt creates an object of class boot.wt with the following components:

Modname	the names of each model of the candidate model set.
K	the number of estimated parameters for each model.
(Q)AIC(c)	the information criteria requested for each model (AICc, AICc, QAIC, QAICc).
Delta_(Q)AIC(c)	the appropriate delta AIC component depending on the information criteria selected.
ModelLik	the relative likelihood of the model given the data (exp(-0.5*delta[i])). This is not to be confused with the likelihood of the parameters given the data. The relative likelihood can then be normalized across all models to get the model probabilities.
(Q)AIC(c)Wt	the Akaike weights, also termed "model probabilities" sensu Burnham and Anderson (2002) and Anderson (2008). These measures indicate the level of support (i.e., weight of evidence) in favor of any given model being the most parsimonious among the candidate model set.
PiWt	the relative frequencies of model selection from the bootstrap.
c.hat	if c.hat was specified as an argument, it is included in the table.

Author(s)

Marc J. Mazerolle

References

Anderson, D. R. (2008) Model-based Inference in the Life Sciences: a primer on evidence. Springer: New York.

Burnham, K. P., Anderson, D. R. (2002) Model Selection and Multimodel Inference: a practical information-theoretic approach. Second edition. Springer: New York.

Burnham, K. P., Anderson, D. R. (2004) Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods and Research 33, 261–304.

Mazerolle, M. J. (2006) Improving data analysis in herpetology: using Akaike's Information Criterion (AIC) to assess the strength of biological hypotheses. Amphibia-Reptilia 27, 169–180.

See Also

AICc, confset, c_hat, evidence, importance, modavg, modavgShrink, modavgPred

Examples

##Mazerolle (2006) frog water loss example
data(dry.frog)

##setup a subset of models of Table 1
Cand.models <- list( )
Cand.models[[1]] <- lm(log_Mass_lost ~ Shade + Substrate +
                       cent_Initial_mass + Initial_mass2,
                       data = dry.frog)
Cand.models[[2]] <- lm(log_Mass_lost ~ Shade + Substrate +
                       cent_Initial_mass + Initial_mass2 +
                       Shade:Substrate, data = dry.frog)
Cand.models[[3]] <- lm(log_Mass_lost ~ cent_Initial_mass +
                       Initial_mass2, data = dry.frog)
Cand.models[[4]] <- lm(log_Mass_lost ~ Shade + cent_Initial_mass +
                       Initial_mass2, data = dry.frog)
Cand.models[[5]] <- lm(log_Mass_lost ~ Substrate + cent_Initial_mass +
                       Initial_mass2, data = dry.frog)

##create a vector of names to trace back models in set
Modnames <- paste("mod", 1:length(Cand.models), sep = " ")

##generate AICc table with bootstrapped relative
##frequencies of model selection
boot.wt(cand.set = Cand.models, modnames = Modnames, sort = TRUE,
        nsim = 10) #number of iterations should be much higher


##Burnham and Anderson (2002) flour beetle data
## Not run: 
data(beetle)
##models as suggested by Burnham and Anderson p. 198          
Cand.set <- list( )
Cand.set[[1]] <- glm(Mortality_rate ~ Dose, family =
                     binomial(link = "logit"), weights = Number_tested,
                     data = beetle)
Cand.set[[2]] <- glm(Mortality_rate ~ Dose, family =
                     binomial(link = "probit"), weights = Number_tested,
                     data = beetle)
Cand.set[[3]] <- glm(Mortality_rate ~ Dose, family =
                     binomial(link ="cloglog"), weights = Number_tested,
                     data = beetle)

##create a vector of names to trace back models in set
Modnames <- paste("Mod", 1:length(Cand.set), sep = " ")

##model selection table with bootstrapped
##relative frequencies
boot.wt(cand.set = Cand.set, modnames = Modnames)

## End(Not run)

---

Bullfrog Occupancy and Common Reed Invasion

Description

This is a data set from Mazerolle et al. (2014) on the occupancy of Bullfrogs (Lithobates catesbeianus) in 50 wetlands sampled in 2009 in the area of Montreal, QC.

Usage

data(bullfrog)

Format

A data frame with 50 observations on the following 23 variables.

Location

a factor with a unique identifier for each wetland.

Reed.presence

a binary variable, either 1 (reed present) or 0 (reed absent).

V1

a binary variable for detection (1) or non detection (0) of bullfrogs during the first survey.

V2

a binary variable for detection (1) or non detection (0) of bullfrogs during the second survey.

V3

a binary variable for detection (1) or non detection (0) of bullfrogs during the third survey.

V4

a binary variable for detection (1) or non detection (0) of bullfrogs during the fourth survey.

V5

a binary variable for detection (1) or non detection (0) of bullfrogs during the fifth survey.

V6

a binary variable for detection (1) or non detection (0) of bullfrogs during the sixth survey.

V7

a binary variable for detection (1) or non detection (0) of bullfrogs during the seventh survey.

Effort1

a numeric variable for the centered number of sampling stations during the first survey.

Effort2

a numeric variable for the centered number of sampling stations during the second survey.

Effort3

a numeric variable for the centered number of sampling stations during the third survey.

Effort4

a numeric variable for the centered number of sampling stations during the fourth survey.

Effort5

a numeric variable for the centered number of sampling stations during the fifth survey.

Effort6

a numeric variable for the centered number of sampling stations during the sixth survey.

Effort7

a numeric variable for the centered number of sampling stations during the seventh survey.

Type1

a binary variable to identify the survey type, either minnow trap (1) or call survey (0) during the first sampling occasion.

Type2

a binary variable to identify the survey type, either minnow trap (1) or call survey (0) during the second sampling occasion.

Type3

a binary variable to identify the survey type, either minnow trap (1) or call survey (0) during the third sampling occasion.

Type4

a binary variable to identify the survey type, either minnow trap (1) or call survey (0) during the fourth sampling occasion.

Type5

a binary variable to identify the survey type, either minnow trap (1) or call survey (0) during the fifth sampling occasion.

Type6

a binary variable to identify the survey type, either minnow trap (1) or call survey (0) during the sixth sampling occasion.

Type7

a binary variable to identify the survey type, either minnow trap (1) or call survey (0) during the seventh sampling occasion.

Details

This data set is used to illustrate single-species single-season occupancy models (MacKenzie et al. 2002) in Mazerolle (2015). The average number of sampling stations on each visit was 8.665714, and was used to center Effort on each visit.

Source

MacKenzie, D. I., Nichols, J. D., Lachman, G. B., Droege, S., Royle, J. A., Langtimm, C. A. (2002). Estimating site occupancy rates when detection probabilities are less than one. Ecology 83, 2248–2255.

Mazerolle, M. J., Perez, A., Brisson, J. (2014) Common reed (Phragmites australis) invasion and amphibian distribution in freshwater wetlands. Wetlands Ecology and Management 22, 325–340.

Mazerolle, M. J. (2015) Estimating detectability and biological parameters of interest with the use of the R environment. Journal of Herpetology 49, 541–559.

Examples

data(bullfrog)
str(bullfrog)

---

Estimate Dispersion for Poisson and Binomial GLM's and GLMM's

Description

Functions to compute an estimate of c-hat for binomial or Poisson GLM's and GLMM's using different estimators of overdispersion.

Usage

c_hat(mod, method = "pearson", ...)

## S3 method for class 'glm'
c_hat(mod, method = "pearson", ...)

## S3 method for class 'glmmTMB'
c_hat(mod, method = "pearson", ...)

## S3 method for class 'merMod'
c_hat(mod, method = "pearson", ...)

## S3 method for class 'vglm'
c_hat(mod, method = "pearson", ...)

Arguments

mod	an object of class glm, glmmTMB, merMod, or vglm for which a c-hat estimate is required.
method	this argument defines the estimator used. The default "pearson" uses the Pearson chi-square divided by the residual degrees of freedom. Other methods include "deviance" consisting of the residual deviance divided by the residual degrees of freedom, "farrington" for the estimator suggested by Farrington (1996), and "fletcher" for the estimator suggested by Fletcher (2012).
...	additional arguments passed to the function.

Details

Poisson and binomial GLM's do not have a parameter for the variance and it is usually held fixed to 1 (i.e., mean = variance). However, one must check whether this assumption is appropriate by estimating the overdispersion parameter (c-hat). Though one can obtain an estimate of c-hat by dividing the residual deviance by the residual degrees of freedom (i.e., method = "deviance"), McCullagh and Nelder (1989) and Venables and Ripley (2002) recommend using Pearson's chi-square divided by the residual degrees of freedom (method = "pearson"). An estimator based on Farrington (1996) is also implemented by the function using the argument method = "farrington". Recent work by Fletcher (2012) suggests that an alternative estimator performs better than the above-mentioned methods in the presence of sparse data and is now implemented with method = "fletcher". For GLMM's, only the Pearson chi-square estimator of overdispersion is currently implemented.

Note that values of c-hat > 1 indicate overdispersion (variance > mean), but that values much higher than 1 (i.e., > 4) probably indicate lack-of-fit. In cases of moderate overdispersion, one usually multiplies the variance-covariance matrix of the estimates by c-hat. As a result, the SE's of the estimates are inflated (c-hat is also known as a variance inflation factor).

In model selection, c-hat should be estimated from the global model of the candidate model set and the same value of c-hat applied to the entire model set. Specifically, a global model is the most complex model which can be simplified to obtain all the other (nested) models of the set. When no single global model exists in the set of models considered, such as when sample size does not allow a complex model, one can estimate c-hat from 'subglobal' models. Here, 'subglobal' models denote models from which only a subset of the models of the candidate set can be derived. In such cases, one can use the smallest value of c-hat for model selection (Burnham and Anderson 2002).

Note that c-hat counts as an additional parameter estimated and should be added to K. All functions in package AICcmodavg automatically add 1 when the c.hat argument > 1 and apply the same value of c-hat for the entire model set. When c.hat > 1, functions compute quasi-likelihood information criteria (either QAICc or QAIC, depending on the value of the second.ord argument) by scaling the log-likelihood of the model by c.hat. The value of c.hat can influence the ranking of the models: as c-hat increases, QAIC or QAICc will favor models with fewer parameters. As an additional check against this potential problem, one can create several model selection tables by incrementing values of c-hat to assess the model selection uncertainty. If ranking changes little up to the c-hat value observed, one can be confident in making inference.

In cases of underdispersion (c-hat < 1), it is recommended to keep the value of c.hat to 1. However, note that values of c-hat << 1 can also indicate lack-of-fit and that an alternative model (and distribution) should be investigated.

Note that c_hat only supports the estimation of c-hat for binomial models with trials > 1 (i.e., success/trial or cbind(success, failure) syntax) or with Poisson GLM's or GLMM's.

Value

c_hat returns an object of class c_hat with the estimated c-hat value and an attribute for the type of estimator used.

Author(s)

Marc J. Mazerolle

References

Anderson, D. R. (2008) Model-based Inference in the Life Sciences: a primer on evidence. Springer: New York.

Burnham, K. P., Anderson, D. R. (2002) Model Selection and Multimodel Inference: a practical information-theoretic approach. Second edition. Springer: New York.

Burnham, K. P., Anderson, D. R. (2004) Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods and Research 33, 261–304.

Farrington, C. P. (1996) On assessing goodness of fit of generalized linear models to sparse data. Journal of the Royal Statistical Society B 58, 349–360.

Fletcher, D. J. (2012) Estimating overdispersion when fitting a generalized linear model to sparse data. Biometrika 99, 230–237.

Mazerolle, M. J. (2006) Improving data analysis in herpetology: using Akaike's Information Criterion (AIC) to assess the strength of biological hypotheses. Amphibia-Reptilia 27, 169–180.

McCullagh, P., Nelder, J. A. (1989) Generalized Linear Models. Second edition. Chapman and Hall: New York.

Venables, W. N., Ripley, B. D. (2002) Modern Applied Statistics with S. Second edition. Springer: New York.

See Also

AICc, confset, evidence, modavg, importance, modavgPred, mb.gof.test, Nmix.gof.test, anovaOD, summaryOD

Examples

#binomial glm example
set.seed(seed = 10)
resp <- rbinom(n = 60, size = 1, prob = 0.5)
set.seed(seed = 10)
treat <- as.factor(sample(c(rep(x = "m", times = 30), rep(x = "f",
                                           times = 30))))
age <- as.factor(c(rep("young", 20), rep("med", 20), rep("old", 20)))
#each invidual has its own response (n = 1)
mod1 <- glm(resp ~ treat + age, family = binomial)
## Not run: 
c_hat(mod1) #gives an error because model not appropriate for
##computation of c-hat

## End(Not run)

##computing table to summarize successes
table(resp, treat, age)
dat2 <- as.data.frame(table(resp, treat, age)) #not quite what we need
data2 <- data.frame(success = c(9, 4, 2, 3, 5, 2),
                    sex = c("f", "m", "f", "m", "f", "m"),
                    age = c("med", "med", "old", "old", "young",
                      "young"), total = c(13, 7, 10, 10, 7, 13))
data2$prop <- data2$success/data2$total
data2$fail <- data2$total - data2$success

##run model with success/total syntax using weights argument
mod2 <- glm(prop ~ sex + age, family = binomial, weights = total,
            data = data2)
c_hat(mod2)

##run model with other syntax cbind(success, fail)
mod3 <- glm(cbind(success, fail) ~ sex + age, family = binomial,
            data = data2) 
c_hat(mod3)

---

Blood Calcium Concentration in Birds

Description

This data set features calcium concentration in the plasma of birds of both sexes following a hormonal treatment.

Usage

data(calcium)

Format

A data frame with 20 rows and 3 variables.

Calcium

calcium concentration in mg/100 ml in the blood of birds.

Hormone

a factor with two levels indicating whether the bird received a hormonal treatment or not.

Sex

a factor with two levels coding for the sex of birds.

Details

Zar (1984, p. 206) illustrates a two-way ANOVA with interaction with this data set.

Source

Zar, J. H. (1984) Biostatistical analysis. Second edition. Prentice Hall: Englewood Cliffs, New Jersey.

Examples

data(calcium)
str(calcium)

---