modavgShrink {AICcmodavg}

Compute Model-averaged Parameter Estimate with Shrinkage (Multimodel Inference)

Description

This function computes an alternative version of model-averaging parameter estimates that consists in shrinking estimates toward 0 to reduce model selection bias as in Burnham and Anderson (2002, p. 152), Anderson (2008, pp. 130-132) and Lukacs et al. (2010). Specifically, models without the parameter of interest have an estimate and variance of 0. modavgShrink also returns unconditional standard errors and unconditional confidence intervals as described in Buckland et al. (1997) and Burnham and Anderson (2002).

Usage

modavgShrink(cand.set, parm, modnames = NULL, second.ord = TRUE,
              nobs = NULL, uncond.se = "revised", conf.level = 0.95,
              ...)
## S3 method for class 'AICaov.lm'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICbetareg'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICsclm.clm'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICclm'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICclmm'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICcoxme'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICcoxph'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICglm.lm'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, gamdisp = NULL, ...)

## S3 method for class 'AICgls'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICglmmTMB'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, ...)

## S3 method for class 'AIChurdle'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AIClm'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AIClme'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AIClmekin'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICmer'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICglmerMod'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AIClmerMod'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AIClmerModLmerTest'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICmaxlikeFit.list'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, ...)

## S3 method for class 'AICmultinom.nnet'
modavgShrink(cand.set, parm, modnames =
        NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, ...)

## S3 method for class 'AICnegbin.glm.lm'
modavgShrink(cand.set, parm, modnames = NULL, 
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICpolr'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICrlm.lm'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICsurvreg'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICvglm'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, ...)

## S3 method for class 'AICzeroinfl'
modavgShrink(cand.set, parm, modnames = NULL,
        second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, ...)

## S3 method for class 'AICunmarkedFitOccu'
modavgShrink(cand.set, parm, modnames =
        NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitColExt'
modavgShrink(cand.set, parm, modnames
           = NULL, second.ord = TRUE, nobs = NULL, uncond.se =
           "revised", conf.level = 0.95, c.hat = 1, parm.type = NULL,
           ...)

## S3 method for class 'AICunmarkedFitOccuRN'
modavgShrink(cand.set, parm, modnames
        = NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitPCount'
modavgShrink(cand.set, parm, modnames
        = NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitPCO'
modavgShrink(cand.set, parm, modnames =
        NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitDS'
modavgShrink(cand.set, parm, modnames =
        NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitGDS'
modavgShrink(cand.set, parm, modnames =
        NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitOccuFP'
modavgShrink(cand.set, parm, modnames
        = NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitMPois'
modavgShrink(cand.set, parm, modnames
        = NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitGMM'
modavgShrink(cand.set, parm, modnames
        = NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitGPC'
modavgShrink(cand.set, parm, modnames
        = NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitOccuMulti'
modavgShrink(cand.set, parm, modnames =
        NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitOccuMS'
modavgShrink(cand.set, parm, modnames =
        NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitOccuTTD'
modavgShrink(cand.set, parm, modnames =
        NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitMMO'
modavgShrink(cand.set, parm, modnames =
        NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

## S3 method for class 'AICunmarkedFitDSO'
modavgShrink(cand.set, parm, modnames =
        NULL, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
        conf.level = 0.95, c.hat = 1, parm.type = NULL, ...)

Arguments

cand.set

a list storing each of the models in the candidate model set.

parm

the parameter of interest, enclosed between quotes, for which a model-averaged estimate is required. For a categorical variable, the label of the estimate must be included as it appears in the output (see 'Details' below).

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 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.

uncond.se

either, "old", or "revised", specifying the equation used to compute the unconditional standard error of a model-averaged estimate. With uncond.se = "old", computations are based on equation 4.9 of Burnham and Anderson (2002), which was the former way to compute unconditional standard errors. With uncond.se = "revised", equation 6.12 of Burnham and Anderson (2002) is used. Anderson (2008, p. 111) recommends use of the revised version for the computation of unconditional standard errors and it is now the default. Note that versions of package AICcmodavg < 1.04 used the old method to compute unconditional standard errors.

conf.level

the confidence level (1 - \alpha) requested for the computation of unconditional confidence intervals.

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, modavgShrink 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.

gamdisp

if gamma GLM is used, the dispersion parameter should be specified here to apply the same value to each model.

parm.type

this argument specifies the parameter type on which the model-averaged estimate of a predictor will be computed and is only relevant for models of unmarkedFit classes. The character strings supported vary with the type of model fitted. For unmarkedFitOccu objects, either psi or detect can be supplied to indicate whether the parameter is on occupancy or detectability, respectively. For unmarkedFitColExt objects, possible values are psi, gamma, epsilon, and detect, for parameters on occupancy in the inital year, colonization, extinction, and detectability, respectively. For unmarkedFitOccuTTD objects, possible values are psi, gamma, epsilon, and detect, for parameters on occupancy in the inital year, colonization, extinction, and time-to-dection (lambda rate parameter), respectively. For unmarkedFitOccuFP objects, one can specify psi, detect, falsepos, and certain, for occupancy, detectability, probability of assigning false-positives, and probability detections are certain, respectively. For unmarkedFitOccuRN objects, either lambda or detect can be entered for abundance and detectability parameters, respectively. For unmarkedFitPCount and unmarkedFitMPois objects, lambda or detect denote parameters on abundance and detectability, respectively. For unmarkedFitPCO, unmarkedFitMMO, and unmarkedFitDSO objects, one can enter lambda, gamma, omega, iota, or detect, to specify parameters on abundance, recruitment, apparent survival, immigration, and detectability, respectively. For unmarkedFitDS objects, lambda and detect are supported. For unmarkedFitGDS, lambda, phi, and detect denote abundance, availability, and detection probability, respectively. For unmarkedFitGMM and unmarkedFitGPC objects, lambda, phi, and detect denote abundance, availability, and detectability, respectively.

...

additional arguments passed to the function.

Details

The parameter for which a model-averaged estimate is requested must be specified with the parm argument and must be identical to its label in the model output (e.g., from summary). For factors, one must specify the name of the variable and the level of interest. The shrinkage version of model averaging is only appropriate for cases where each parameter is given an equal weighting in the model (i.e., each parameter must appear the same number of times in the models) and has the same interpretation across all models. As a result, models with interaction terms or polynomial terms are not supported by modavgShrink.

modavgShrink is implemented for a list containing objects of aov, betareg, clm, clmm, clogit, coxme, coxph, glm, glmmTMB, gls, hurdle, lm, lme, lmekin, maxlikeFit, mer, glmerMod, lmerMod, lmerModLmerTest, multinom, polr, rlm, survreg, vglm, zeroinfl classes as well as various models of unmarkedFit classes.

Value

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

Parameter

the parameter for which a model-averaged estimate with shrinkage was obtained.

Mod.avg.table

the model selection table based on models including the parameter of interest.

Mod.avg.beta

the model-averaged estimate based on all models.

Uncond.SE

the unconditional standard error for the model-averaged estimate (as opposed to the conditional SE based on a single model).

Conf.level

the confidence level used to compute the confidence interval.

Lower.CL

the lower confidence limit.

Upper.CL

the upper confidence limit.

Author(s)

Marc J. Mazerolle

References

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

Buckland, S. T., Burnham, K. P., Augustin, N. H. (1997) Model selection: an integral part of inference. Biometrics 53, 603–618.

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.

Lukacs, P. M., Burnham, K. P., Anderson, D. R. (2010) Model selection bias and Freedman's paradox. Annals of the Institute of Statistical Mathematics 62, 117–125.

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, aictab, c_hat, importance, confset, evidence, modavg, modavgCustom, modavgPred

Examples

##cement example in Burnham and Anderson 2002
data(cement)
##setup same model set as in Table 3.2, p. 102         
Cand.models <- list( )
Cand.models[[1]] <- lm(y ~ x1 + x2, data = cement)
Cand.models[[2]] <- lm(y ~ x1 + x2 + x4, data = cement)          
Cand.models[[3]] <- lm(y ~ x1 + x2 + x3, data = cement)
Cand.models[[4]] <- lm(y ~ x1 + x4, data = cement)
Cand.models[[5]] <- lm(y ~ x1 + x3 + x4, data = cement)
Cand.models[[6]] <- lm(y ~ x2 + x3 + x4, data = cement)
Cand.models[[7]] <- lm(y ~ x1 + x2 + x3 + x4, data = cement)
Cand.models[[8]] <- lm(y ~ x3 + x4, data = cement)
Cand.models[[9]] <- lm(y ~ x2 + x3, data = cement)
Cand.models[[10]] <- lm(y ~ x4, data = cement)
Cand.models[[11]] <- lm(y ~ x2, data = cement)
Cand.models[[12]] <- lm(y ~ x2 + x4, data = cement)
Cand.models[[13]] <- lm(y ~ x1, data = cement)
Cand.models[[14]] <- lm(y ~ x1 + x3, data = cement)
Cand.models[[15]] <- lm(y ~ x3, data = cement)

##vector of model names
Modnames <- paste("mod", 1:15, sep="")

##AICc          
aictab(cand.set = Cand.models, modnames = Modnames)

##compute model-averaged estimate with shrinkage - each parameter
##appears 8 times in the models 
modavgShrink(cand.set = Cand.models, modnames = Modnames, parm = "x1")

##compare against classic model-averaging
modavg(cand.set = Cand.models, modnames = Modnames, parm = "x1")
##note that model-averaged estimate with shrinkage is closer to 0 than
##with the classic version

##remove a few models from the set and run again
Cand.unbalanced <- Cand.models[-c(3, 14, 15)]

##set up model names
Modnames <- paste("mod", 1:length(Cand.unbalanced), sep="")

##issues an error because some parameters appear more often than others
## Not run: modavgShrink(cand.set = Cand.unbalanced,
                       modnames = Modnames, parm = "x1")
## End(Not run)



##example on Orthodont data set in nlme
## Not run: 
require(nlme)

##set up candidate model list
##age and sex parameters appear in the same number of models
##same number of models with and without these parameters
Cand.models <- list( )
Cand.models[[1]] <- lme(distance ~ age, data = Orthodont, method = "ML") 
##random is ~ age | Subject as it is a grouped data frame
Cand.models[[2]] <- lme(distance ~ age + Sex, data = Orthodont,
                        random = ~ 1, method = "ML")
Cand.models[[3]] <- lme(distance ~ 1, data = Orthodont, random = ~ 1, 
                        method = "ML") 
Cand.models[[4]] <- lme(distance ~ Sex, data = Orthodont, random = ~ 1,
                        method = "ML")  

##create a vector of model names
Modnames <- paste("mod", 1:length(Cand.models), sep = "")

##compute importance values for age
imp.age <- importance(cand.set = Cand.models, parm = "age",
                      modnames = Modnames, second.ord = TRUE,
                      nobs = NULL)

##compute shrinkage version of model averaging on age
mod.avg.age.shrink <- modavgShrink(cand.set = Cand.models,
                                    parm = "age", modnames = Modnames,
                                    second.ord = TRUE, nobs = NULL)

##compute classic version of model averaging on age
mod.avg.age.classic <- modavg(cand.set = Cand.models, parm = "age",
                              modnames = Modnames, second.ord = TRUE,
                              nobs = NULL)

##correspondence between shrinkage version and classic version of
##model averaging 
mod.avg.age.shrink$Mod.avg.beta/imp.age$w.plus
mod.avg.age.classic$Mod.avg.beta
detach(package:nlme)

## End(Not run)


##example of N-mixture model modified from ?pcount
## Not run: 
require(unmarked)
data(mallard)
mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs = mallard.site,
                                  obsCovs = mallard.obs)
##set up models so that each variable on abundance appears twice
fm.mall.one <- pcount(~ ivel + date  ~ length + forest, mallardUMF,
                      K = 30)
fm.mall.two <- pcount(~ ivel + date  ~ elev + forest, mallardUMF,
                      K = 30)
fm.mall.three <- pcount(~ ivel + date  ~ length + elev, mallardUMF,
                        K = 30)

##model list and names
Cands <- list(fm.mall.one, fm.mall.two, fm.mall.three)
Modnames <- c("length + forest", "elev + forest", "length + elev")

##compute model-averaged estimate with shrinkage for elev on abundance
modavgShrink(cand.set = Cands, modnames = Modnames, parm = "elev",
              parm.type = "lambda")
detach(package:unmarked)

## End(Not run)