This function model-averages the estimate of a parameter of interest
among a set of candidate models, computes the unconditional standard
error and unconditional confidence intervals as described in Buckland et
al. (1997) and Burnham and Anderson (2002). This model-averaged estimate
is also referred to as a natural average of the estimate by Burnham and
Anderson (2002, p. 152).
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 - α) requested for the computation of
unconditional confidence intervals.
exclude
this argument excludes models based on the terms specified for the
computation of a model-averaged estimate of parm. The
exclude argument is set to NULL by default and does not
exclude any models other than those without the parm. When
parm is a main effect but is also involved in
interactions/polynomial terms in some models, one should specify the
interaction/polynomial terms as a list to exclude models with these
terms from the computation of model-averaged estimate of the main effect
(e.g., exclude = list("sex:mass", "mass2")). See 'Details'
and 'Examples' below.
warn
logical. If TRUE, modavg performs a check and isssues a
warning when the value in parm occurs more than once in any given
model. This is a check for potential interaction/polynomial terms in
the model when such terms are constructed with the usual operators
(e.g., I( ) for polynomial terms, : for interaction
terms).
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 effect size
will be computed and is only relevant for models of
unmarkedFitOccu, unmarkedFitColExt,
unmarkedFitOccuFP, unmarkedFitOccuRN,
unmarkedFitMPois, unmarkedFitPCount,
unmarkedFitPCO, unmarkedFitDS, unmarkedFitGDS,
unmarkedFitGMM, and unmarkedFitGPC 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, possible
values are psi, gamma, epsilon, and
detect, for parameters on occupancy in the inital year,
colonization, extinction, and detectability, respectively. For
unmarkedFitOccuFP objects, one can specify psi,
detect, or fp, for occupancy, detectability, and
probability of assigning false-positives, 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 objects, one can enter lambda,
gamma, omega, or detect, to specify parameters on
abundance, recruitment, apparent survival, and detectability,
respectively. For unmarkedFitDS objects, only lambda is
supported for the moment. For unmarkedFitGDS, lambda
and phi denote abundance and availability, 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.
modavg includes checks to find variations of interaction terms
specified in the parm and exclude arguments. However, to
avoid problems, one should specify interaction terms consistently for
all models: e.g., either a:b or b:a for all models, but
not a mixture of both.
You must exercise caution when some models include interaction or
polynomial terms, because main effect terms do not have the same
interpretation when they also appear in an interaction/polynomial term
in the same model. In such cases, one should exclude models containing
interaction terms where the main effect is involved with the
exclude argument of modavg. Note that modavg
checks for potential cases of multiple instances of a variable appearing
more than once in a given model (presumably in an interaction) and
issues a warning. To correctly compute the model-averaged estimate of a
main effect involved in interaction/polynomial terms, specify the
interaction terms(s) that should not appear in the same model with the
exclude argument. This will effectively exclude models from the
computation of the model-averaged estimate.
When warn = TRUE, modavg looks for matches among the
labels of the estimates with identical. It then compares the
results to partial matches with regexpr, and issues a warning
whenever they are different. As a result, modavg may issue a
warning when some variables or levels of categorical variables have
nested names (e.g., treat, treat10; L, TL).
When this warning is only due to the presence of similarly named
variables in the models (and NOT due to interaction terms), you can
suppress this warning by setting warn = FALSE.
The model-averaging estimator implemented in modavg is known to
be biased away from 0 when there is substantial model selection
uncertainty (Cade 2015). In such instances, it is recommended to use
the model-averaging shrinkage estimator (i.e., modavgShrink) for
inference on beta estimates or to focus on model-averaged effect sizes
(modavgEffect) and model-averaged predictions
(modavgPred).
modavg is implemented for a list containing objects of
aov, betareg, clm, clmm, clogit,
coxme, coxph, glm, gls, hurdle,
lm, lme, lmekin, maxlikeFit, mer,
glmerMod, lmerMod, multinom, polr,
rlm, survreg, vglm, zeroinfl classes as well
as various models of unmarkedFit classes.
Value
modavg creates an object of class modavg with the following
components:
Parameter
the parameter for which a model-averaged estimate was
obtained.
Mod.avg.table
the reduced model selection table based on models
including the parameter of interest.
Mod.avg.beta
the model-averaged estimate based on all models
including the parameter of interest (see 'Details' above regarding the
exclusion of models where parameter of interest is involved in an
interaction).
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. Biometrics53, 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 Research33, 261–304.
Cade, B. S. (2015) Model averaging and muddled multimodel
inferences. Ecology96, 2370–2382.
Dail, D., Madsen, L. (2011) Models for estimating abundance from
repeated counts of an open population. Biometrics67,
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
Monographs62, 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. Ecology83,
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. Ecology84, 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-Reptilia27, 169–180.
Royle, J. A. (2004) N-mixture models for estimating population
size from spatially replicated counts. Biometrics60,
108–115.
##anuran larvae example modified from Mazerolle (2006)
##these are different models than in the paper
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
Cand.mod <- list( )
##global model
Cand.mod[[1]] <- glm(Num_anura ~ Type + log.Perimeter +
Type:log.Perimeter + Num_ranatra,
family = poisson, offset = log(Effort),
data = min.trap)
##interactive model
Cand.mod[[2]] <- glm(Num_anura ~ Type + log.Perimeter +
Type:log.Perimeter, family = poisson,
offset = log(Effort), data = min.trap)
##additive model
Cand.mod[[3]] <- glm(Num_anura ~ Type + log.Perimeter, family = poisson,
offset = log(Effort), data = min.trap)
##Predator model
Cand.mod[[4]] <- glm(Num_anura ~ Type + Num_ranatra, family = poisson,
offset = log(Effort), data = min.trap)
##check c-hat for global model
c_hat(Cand.mod[[1]]) #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
##assign names to each model
Modnames <- c("global model", "interactive model",
"additive model", "invertpred model")
##model selection
aictab(Cand.mod, Modnames)
##compute model-averaged estimates for parameters appearing in top
##models
modavg(parm = "Num_ranatra", cand.set = Cand.mod, modnames = Modnames)
##round to 4 digits after decimal point
print(modavg(parm = "Num_ranatra", cand.set = Cand.mod,
modnames = Modnames), digits = 4)
##model-averaging a variable involved in an interaction
##the following produces an error - because the variable is involved
##in an interaction in some candidate models
## Not run: modavg(parm = "TypeBOG", cand.set = Cand.mod,
modnames = Modnames)
## End(Not run)
##exclude models where the variable is involved in an interaction
##to get model-averaged estimate of main effect
modavg(parm = "TypeBOG", cand.set = Cand.mod, modnames = Modnames,
exclude = list("Type:log.Perimeter"))
##to get model-averaged estimate of interaction
modavg(parm = "TypeBOG:log.Perimeter", cand.set = Cand.mod,
modnames = Modnames)
##beware of variables that have similar names
set.seed(seed = 4)
resp <- rnorm(n = 40, mean = 3, sd = 1)
size <- rep(c("small", "medsmall", "high", "medhigh"), times = 10)
set.seed(seed = 4)
mass <- rnorm(n = 40, mean = 2, sd = 0.1)
mass2 <- mass^2
age <- rpois(n = 40, lambda = 3.2)
agecorr <- rpois(n = 40, lambda = 2)
sizecat <- rep(c("a", "ab"), times = 20)
data1 <- data.frame(resp = resp, size = size, sizecat = sizecat,
mass = mass, mass2 = mass2, age = age,
agecorr = agecorr)
##set up models in list
Cand <- list( )
Cand[[1]] <- lm(resp ~ size + agecorr, data = data1)
Cand[[2]] <- lm(resp ~ size + mass + agecorr, data = data1)
Cand[[3]] <- lm(resp ~ age + mass, data = data1)
Cand[[4]] <- lm(resp ~ age + mass + mass2, data = data1)
Cand[[5]] <- lm(resp ~ mass + mass2 + size, data = data1)
Cand[[6]] <- lm(resp ~ mass + mass2 + sizecat, data = data1)
Cand[[7]] <- lm(resp ~ sizecat, data = data1)
Cand[[8]] <- lm(resp ~ sizecat + mass + sizecat:mass, data = data1)
Cand[[9]] <- lm(resp ~ agecorr + sizecat + mass + sizecat:mass,
data = data1)
##create vector of model names
Modnames <- paste("mod", 1:length(Cand), sep = "")
aictab(cand.set = Cand, modnames = Modnames, sort = TRUE) #correct
##as expected, issues warning as mass occurs sometimes with "mass2" or
##"sizecatab:mass" in some of the models
## Not run: modavg(cand.set = Cand, parm = "mass", modnames = Modnames)
##no warning issued, because "age" and "agecorr" never appear in same model
modavg(cand.set = Cand, parm = "age", modnames = Modnames)
##as expected, issues warning because warn=FALSE, but it is a very bad
##idea in this example since "mass" occurs with "mass2" and "sizecat:mass"
##in some of the models - results are INCORRECT
## Not run: modavg(cand.set = Cand, parm = "mass", modnames = Modnames,
warn = FALSE)
## End(Not run)
##correctly excludes models with quadratic term and interaction term
##results are CORRECT
modavg(cand.set = Cand, parm = "mass", modnames = Modnames,
exclude = list("mass2", "sizecat:mass"))
##correctly computes model-averaged estimate because no other parameter
##occurs simultaneously in any of the models
modavg(cand.set = Cand, parm = "sizesmall", modnames = Modnames) #correct
##as expected, issues a warning because "sizecatab" occurs sometimes in
##an interaction in some models
## Not run: modavg(cand.set = Cand, parm = "sizecatab",
modnames = Modnames)
## End(Not run)
##exclude models with "sizecat:mass" interaction - results are CORRECT
modavg(cand.set = Cand, parm = "sizecatab", modnames = Modnames,
exclude = list("sizecat:mass"))
##example with multiple-season occupancy model modified from ?colext
##this is a bit longer
## Not run:
require(unmarked)
data(frogs)
umf <- formatMult(masspcru)
obsCovs(umf) <- scale(obsCovs(umf))
siteCovs(umf) <- rnorm(numSites(umf))
yearlySiteCovs(umf) <- data.frame(year = factor(rep(1:7,
numSites(umf))))
##set up model with constant transition rates
fm <- colext(psiformula = ~ 1, gammaformula = ~ 1, epsilonformula = ~ 1,
pformula = ~ JulianDate + I(JulianDate^2), data = umf,
control = list(trace=1, maxit=1e4))
##model with with year-dependent transition rates
fm.yearly <- colext(psiformula = ~ 1, gammaformula = ~ year,
epsilonformula = ~ year,
pformula = ~ JulianDate + I(JulianDate^2),
data = umf)
##store in list and assign model names
Cand.mods <- list(fm, fm.yearly)
Modnames <- c("psi1(.)gam(.)eps(.)p(Date + Date2)",
"psi1(.)gam(Year)eps(Year)p(Date + Date2)")
##compute model-averaged estimate of occupancy in the first year
modavg(cand.set = Cand.mods, modnames = Modnames, parm = "(Intercept)",
parm.type = "psi")
##compute model-averaged estimate of Julian Day squared on detectability
modavg(cand.set = Cand.mods, modnames = Modnames,
parm = "I(JulianDate^2)", parm.type = "detect")
## End(Not run)
##example of model-averaged estimate of area from distance model
##this is a bit longer
## Not run:
data(linetran) #example modified from ?distsamp
ltUMF <- with(linetran, {
unmarkedFrameDS(y = cbind(dc1, dc2, dc3, dc4),
siteCovs = data.frame(Length, area, habitat),
dist.breaks = c(0, 5, 10, 15, 20),
tlength = linetran$Length * 1000, survey = "line", unitsIn = "m")
})
## Half-normal detection function. Density output (log scale). No covariates.
fm1 <- distsamp(~ 1 ~ 1, ltUMF)
## Halfnormal. Covariates affecting both density and detection.
fm2 <- distsamp(~ area + habitat ~ area + habitat, ltUMF)
## Hazard function. Covariates affecting both density and detection.
fm3 <- distsamp(~ habitat ~ area + habitat, ltUMF, keyfun="hazard")
##assemble model list
Cands <- list(fm1, fm2, fm3)
Modnames <- paste("mod", 1:length(Cands), sep = "")
##model-average estimate of area on abundance
modavg(cand.set = Cands, modnames = Modnames, parm = "area", parm.type = "lambda")
detach(package:unmarked)
## End(Not run)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(AICcmodavg)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/AICcmodavg/modavg.Rd_%03d_medium.png", width=480, height=480)
> ### Name: modavg
> ### Title: Compute Model-averaged Parameter Estimate (Multimodel Inference)
> ### Aliases: modavg modavg.default modavg.AICaov.lm modavg.AICbetareg
> ### modavg.AICsclm.clm modavg.AICclmm modavg.AICcoxme modavg.AICcoxph
> ### modavg.AICglm.lm modavg.AICgls modavg.AIChurdle modavg.AIClm
> ### modavg.AIClme modavg.AIClmekin modavg.AICmaxlikeFit.list
> ### modavg.AICmer modavg.AIClmerMod modavg.AICglmerMod
> ### modavg.AICmultinom.nnet modavg.AICpolr modavg.AICrlm.lm
> ### modavg.AICsurvreg modavg.AICvglm modavg.AICzeroinfl
> ### modavg.AICunmarkedFitOccu modavg.AICunmarkedFitColExt
> ### modavg.AICunmarkedFitOccuRN modavg.AICunmarkedFitPCount
> ### modavg.AICunmarkedFitPCO modavg.AICunmarkedFitDS
> ### modavg.AICunmarkedFitGDS modavg.AICunmarkedFitOccuFP
> ### modavg.AICunmarkedFitMPois modavg.AICunmarkedFitGMM
> ### modavg.AICunmarkedFitGPC print.modavg
> ### Keywords: models
>
> ### ** Examples
>
> ##anuran larvae example modified from Mazerolle (2006)
> ##these are different models than in the paper
> 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
> Cand.mod <- list( )
> ##global model
> Cand.mod[[1]] <- glm(Num_anura ~ Type + log.Perimeter +
+ Type:log.Perimeter + Num_ranatra,
+ family = poisson, offset = log(Effort),
+ data = min.trap)
> ##interactive model
> Cand.mod[[2]] <- glm(Num_anura ~ Type + log.Perimeter +
+ Type:log.Perimeter, family = poisson,
+ offset = log(Effort), data = min.trap)
> ##additive model
> Cand.mod[[3]] <- glm(Num_anura ~ Type + log.Perimeter, family = poisson,
+ offset = log(Effort), data = min.trap)
> ##Predator model
> Cand.mod[[4]] <- glm(Num_anura ~ Type + Num_ranatra, family = poisson,
+ offset = log(Effort), data = min.trap)
>
> ##check c-hat for global model
> c_hat(Cand.mod[[1]]) #uses Pearson's chi-square/df
'c-hat' 1.1 (method: pearson estimator)
> ##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
>
> ##assign names to each model
> Modnames <- c("global model", "interactive model",
+ "additive model", "invertpred model")
>
> ##model selection
> aictab(Cand.mod, Modnames)
Model selection based on AICc:
K AICc Delta_AICc AICcWt Cum.Wt LL
invertpred model 3 54.03 0.00 0.87 0.87 -23.42
additive model 3 59.38 5.35 0.06 0.93 -26.09
global model 5 59.53 5.50 0.06 0.98 -23.10
interactive model 4 61.77 7.74 0.02 1.00 -25.83
>
> ##compute model-averaged estimates for parameters appearing in top
> ##models
> modavg(parm = "Num_ranatra", cand.set = Cand.mod, modnames = Modnames)
Multimodel inference on "Num_ranatra" based on AICc
AICc table used to obtain model-averaged estimate:
K AICc Delta_AICc AICcWt Estimate SE
global model 5 59.53 5.5 0.06 0.76 0.35
invertpred model 3 54.03 0.0 0.94 0.62 0.25
Model-averaged estimate: 0.63
Unconditional SE: 0.26
95% Unconditional confidence interval: 0.13, 1.13
> ##round to 4 digits after decimal point
> print(modavg(parm = "Num_ranatra", cand.set = Cand.mod,
+ modnames = Modnames), digits = 4)
Multimodel inference on "Num_ranatra" based on AICc
AICc table used to obtain model-averaged estimate:
K AICc Delta_AICc AICcWt Estimate SE
global model 5 59.5331 5.5012 0.0601 0.7594 0.3504
invertpred model 3 54.0319 0.0000 0.9399 0.6231 0.2473
Model-averaged estimate: 0.6312
Unconditional SE: 0.2567
95% Unconditional confidence interval: 0.1281, 1.1344
>
> ##model-averaging a variable involved in an interaction
> ##the following produces an error - because the variable is involved
> ##in an interaction in some candidate models
> ## Not run:
> ##D modavg(parm = "TypeBOG", cand.set = Cand.mod,
> ##D modnames = Modnames)
> ## End(Not run)
>
>
> ##exclude models where the variable is involved in an interaction
> ##to get model-averaged estimate of main effect
> modavg(parm = "TypeBOG", cand.set = Cand.mod, modnames = Modnames,
+ exclude = list("Type:log.Perimeter"))
Multimodel inference on "TypeBOG" based on AICc
AICc table used to obtain model-averaged estimate:
K AICc Delta_AICc AICcWt Estimate SE
additive model 3 59.38 5.35 0.06 -1.70 0.59
invertpred model 3 54.03 0.00 0.94 -1.35 0.56
Model-averaged estimate: -1.37
Unconditional SE: 0.57
95% Unconditional confidence interval: -2.48, -0.26
>
> ##to get model-averaged estimate of interaction
> modavg(parm = "TypeBOG:log.Perimeter", cand.set = Cand.mod,
+ modnames = Modnames)
Multimodel inference on "TypeBOG:log.Perimeter" based on AICc
AICc table used to obtain model-averaged estimate:
K AICc Delta_AICc AICcWt Estimate SE
global model 5 59.53 0.00 0.75 -0.52 0.99
interactive model 4 61.77 2.24 0.25 -0.69 0.95
Model-averaged estimate: -0.56
Unconditional SE: 0.99
95% Unconditional confidence interval: -2.49, 1.37
>
>
>
> ##beware of variables that have similar names
> set.seed(seed = 4)
> resp <- rnorm(n = 40, mean = 3, sd = 1)
> size <- rep(c("small", "medsmall", "high", "medhigh"), times = 10)
> set.seed(seed = 4)
> mass <- rnorm(n = 40, mean = 2, sd = 0.1)
> mass2 <- mass^2
> age <- rpois(n = 40, lambda = 3.2)
> agecorr <- rpois(n = 40, lambda = 2)
> sizecat <- rep(c("a", "ab"), times = 20)
> data1 <- data.frame(resp = resp, size = size, sizecat = sizecat,
+ mass = mass, mass2 = mass2, age = age,
+ agecorr = agecorr)
>
> ##set up models in list
> Cand <- list( )
> Cand[[1]] <- lm(resp ~ size + agecorr, data = data1)
> Cand[[2]] <- lm(resp ~ size + mass + agecorr, data = data1)
> Cand[[3]] <- lm(resp ~ age + mass, data = data1)
> Cand[[4]] <- lm(resp ~ age + mass + mass2, data = data1)
> Cand[[5]] <- lm(resp ~ mass + mass2 + size, data = data1)
> Cand[[6]] <- lm(resp ~ mass + mass2 + sizecat, data = data1)
> Cand[[7]] <- lm(resp ~ sizecat, data = data1)
> Cand[[8]] <- lm(resp ~ sizecat + mass + sizecat:mass, data = data1)
> Cand[[9]] <- lm(resp ~ agecorr + sizecat + mass + sizecat:mass,
+ data = data1)
>
> ##create vector of model names
> Modnames <- paste("mod", 1:length(Cand), sep = "")
>
> aictab(cand.set = Cand, modnames = Modnames, sort = TRUE) #correct
Model selection based on AICc:
K AICc Delta_AICc AICcWt Cum.Wt LL
mod3 4 -2525.96 0.00 0.57 0.57 1267.55
mod4 5 -2523.90 2.06 0.20 0.77 1267.83
mod6 5 -2521.86 4.10 0.07 0.85 1266.81
mod2 7 -2521.14 4.82 0.05 0.90 1269.32
mod8 5 -2520.75 5.21 0.04 0.94 1266.26
mod5 7 -2520.40 5.55 0.04 0.98 1268.95
mod9 6 -2519.68 6.28 0.02 1.00 1267.11
mod7 3 105.00 2630.95 0.00 1.00 -49.16
mod1 6 110.32 2636.28 0.00 1.00 -47.89
>
> ##as expected, issues warning as mass occurs sometimes with "mass2" or
> ##"sizecatab:mass" in some of the models
> ## Not run: modavg(cand.set = Cand, parm = "mass", modnames = Modnames)
>
> ##no warning issued, because "age" and "agecorr" never appear in same model
> modavg(cand.set = Cand, parm = "age", modnames = Modnames)
Multimodel inference on "age" based on AICc
AICc table used to obtain model-averaged estimate:
K AICc Delta_AICc AICcWt Estimate SE
mod3 4 -2525.96 0.00 0.74 0 0
mod4 5 -2523.90 2.06 0.26 0 0
Model-averaged estimate: 0
Unconditional SE: 0
95% Unconditional confidence interval: 0, 0
>
> ##as expected, issues warning because warn=FALSE, but it is a very bad
> ##idea in this example since "mass" occurs with "mass2" and "sizecat:mass"
> ##in some of the models - results are INCORRECT
> ## Not run:
> ##D modavg(cand.set = Cand, parm = "mass", modnames = Modnames,
> ##D warn = FALSE)
> ## End(Not run)
>
> ##correctly excludes models with quadratic term and interaction term
> ##results are CORRECT
> modavg(cand.set = Cand, parm = "mass", modnames = Modnames,
+ exclude = list("mass2", "sizecat:mass"))
Multimodel inference on "mass" based on AICc
AICc table used to obtain model-averaged estimate:
K AICc Delta_AICc AICcWt Estimate SE
mod2 7 -2521.14 4.82 0.08 10 0
mod3 4 -2525.96 0.00 0.92 10 0
Model-averaged estimate: 10
Unconditional SE: 0
95% Unconditional confidence interval: 10, 10
>
> ##correctly computes model-averaged estimate because no other parameter
> ##occurs simultaneously in any of the models
> modavg(cand.set = Cand, parm = "sizesmall", modnames = Modnames) #correct
Multimodel inference on "sizesmall" based on AICc
AICc table used to obtain model-averaged estimate:
K AICc Delta_AICc AICcWt Estimate SE
mod1 6 110.32 2631.46 0.00 0.14 0.39
mod2 7 -2521.14 0.00 0.59 0.00 0.00
mod5 7 -2520.40 0.73 0.41 0.00 0.00
Model-averaged estimate: 0
Unconditional SE: 0
95% Unconditional confidence interval: 0, 0
>
> ##as expected, issues a warning because "sizecatab" occurs sometimes in
> ##an interaction in some models
> ## Not run:
> ##D modavg(cand.set = Cand, parm = "sizecatab",
> ##D modnames = Modnames)
> ## End(Not run)
>
> ##exclude models with "sizecat:mass" interaction - results are CORRECT
> modavg(cand.set = Cand, parm = "sizecatab", modnames = Modnames,
+ exclude = list("sizecat:mass"))
Multimodel inference on "sizecatab" based on AICc
AICc table used to obtain model-averaged estimate:
K AICc Delta_AICc AICcWt Estimate SE
mod6 5 -2521.86 0.00 1 0.00 0.00
mod7 3 105.00 2626.85 0 -0.19 0.27
Model-averaged estimate: 0
Unconditional SE: 0
95% Unconditional confidence interval: 0, 0
>
>
>
> ##example with multiple-season occupancy model modified from ?colext
> ##this is a bit longer
> ## Not run:
> ##D require(unmarked)
> ##D data(frogs)
> ##D umf <- formatMult(masspcru)
> ##D obsCovs(umf) <- scale(obsCovs(umf))
> ##D siteCovs(umf) <- rnorm(numSites(umf))
> ##D yearlySiteCovs(umf) <- data.frame(year = factor(rep(1:7,
> ##D numSites(umf))))
> ##D
> ##D ##set up model with constant transition rates
> ##D fm <- colext(psiformula = ~ 1, gammaformula = ~ 1, epsilonformula = ~ 1,
> ##D pformula = ~ JulianDate + I(JulianDate^2), data = umf,
> ##D control = list(trace=1, maxit=1e4))
> ##D
> ##D ##model with with year-dependent transition rates
> ##D fm.yearly <- colext(psiformula = ~ 1, gammaformula = ~ year,
> ##D epsilonformula = ~ year,
> ##D pformula = ~ JulianDate + I(JulianDate^2),
> ##D data = umf)
> ##D
> ##D ##store in list and assign model names
> ##D Cand.mods <- list(fm, fm.yearly)
> ##D Modnames <- c("psi1(.)gam(.)eps(.)p(Date + Date2)",
> ##D "psi1(.)gam(Year)eps(Year)p(Date + Date2)")
> ##D
> ##D ##compute model-averaged estimate of occupancy in the first year
> ##D modavg(cand.set = Cand.mods, modnames = Modnames, parm = "(Intercept)",
> ##D parm.type = "psi")
> ##D
> ##D ##compute model-averaged estimate of Julian Day squared on detectability
> ##D modavg(cand.set = Cand.mods, modnames = Modnames,
> ##D parm = "I(JulianDate^2)", parm.type = "detect")
> ## End(Not run)
>
>
> ##example of model-averaged estimate of area from distance model
> ##this is a bit longer
> ## Not run:
> ##D data(linetran) #example modified from ?distsamp
> ##D
> ##D ltUMF <- with(linetran, {
> ##D unmarkedFrameDS(y = cbind(dc1, dc2, dc3, dc4),
> ##D siteCovs = data.frame(Length, area, habitat),
> ##D dist.breaks = c(0, 5, 10, 15, 20),
> ##D tlength = linetran$Length * 1000, survey = "line", unitsIn = "m")
> ##D })
> ##D
> ##D ## Half-normal detection function. Density output (log scale). No covariates.
> ##D fm1 <- distsamp(~ 1 ~ 1, ltUMF)
> ##D
> ##D ## Halfnormal. Covariates affecting both density and detection.
> ##D fm2 <- distsamp(~ area + habitat ~ area + habitat, ltUMF)
> ##D
> ##D ## Hazard function. Covariates affecting both density and detection.
> ##D fm3 <- distsamp(~ habitat ~ area + habitat, ltUMF, keyfun="hazard")
> ##D
> ##D ##assemble model list
> ##D Cands <- list(fm1, fm2, fm3)
> ##D Modnames <- paste("mod", 1:length(Cands), sep = "")
> ##D
> ##D ##model-average estimate of area on abundance
> ##D modavg(cand.set = Cands, modnames = Modnames, parm = "area", parm.type = "lambda")
> ##D detach(package:unmarked)
> ## End(Not run)
>
>
>
>
>
> dev.off()
null device
1
>