Compute Modelaveraged Effect Sizes (Multimodel Inference on Group Differences)
Description
This function modelaverages the effect size between two groups defined by a categorical variable based on the entire model set and computes the unconditional standard error and unconditional confidence intervals as described in Buckland et al. (1997) and Burnham and Anderson (2002). This can be particularly useful when dealing with data from an experiment (e.g., ANOVA) and when the focus is to determine the effect of a given factor. This is an informationtheoretic alternative to multiple comparisons (e.g., Burnham et al. 2011).
Usage
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120  modavgEffect(cand.set, modnames = NULL, newdata, second.ord = TRUE,
nobs = NULL, uncond.se = "revised", conf.level = 0.95,
...)
## S3 method for class 'AICaov.lm'
modavgEffect(cand.set, modnames = NULL, newdata,
second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, ...)
## S3 method for class 'AICglm.lm'
modavgEffect(cand.set, modnames = NULL, newdata,
second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, type = "response", c.hat = 1, gamdisp = NULL,
...)
## S3 method for class 'AICgls'
modavgEffect(cand.set, modnames = NULL, newdata,
second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, ...)
## S3 method for class 'AIClm'
modavgEffect(cand.set, modnames = NULL, newdata,
second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, ...)
## S3 method for class 'AIClme'
modavgEffect(cand.set, modnames = NULL, newdata,
second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, ...)
## S3 method for class 'AICmer'
modavgEffect(cand.set, modnames = NULL, newdata,
second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, type = "response", ...)
## S3 method for class 'AICglmerMod'
modavgEffect(cand.set, modnames = NULL,
newdata, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, type = "response", ...)
## S3 method for class 'AIClmerMod'
modavgEffect(cand.set, modnames = NULL,
newdata, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, ...)
## S3 method for class 'AICrlm.lm'
modavgEffect(cand.set, modnames = NULL, newdata,
second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, ...)
## S3 method for class 'AICsurvreg'
modavgEffect(cand.set, modnames = NULL, newdata,
second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, type = "response", ...)
## S3 method for class 'AICunmarkedFitOccu'
modavgEffect(cand.set, modnames = NULL,
newdata, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, type = "response", c.hat = 1,
parm.type = NULL, ...)
## S3 method for class 'AICunmarkedFitColExt'
modavgEffect(cand.set, modnames =
NULL, newdata, second.ord = TRUE, nobs = NULL, uncond.se =
"revised", conf.level = 0.95, type = "response",
c.hat = 1, parm.type = NULL, ...)
## S3 method for class 'AICunmarkedFitOccuRN'
modavgEffect(cand.set, modnames =
NULL, newdata, second.ord = TRUE, nobs = NULL, uncond.se =
"revised", conf.level = 0.95, type = "response",
c.hat = 1, parm.type = NULL, ...)
## S3 method for class 'AICunmarkedFitPCount'
modavgEffect(cand.set, modnames =
NULL, newdata, second.ord = TRUE, nobs = NULL, uncond.se =
"revised", conf.level = 0.95, type = "response",
c.hat = 1, parm.type = NULL, ...)
## S3 method for class 'AICunmarkedFitPCO'
modavgEffect(cand.set, modnames = NULL,
newdata, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, type = "response", c.hat = 1,
parm.type = NULL, ...)
## S3 method for class 'AICunmarkedFitDS'
modavgEffect(cand.set, modnames = NULL,
newdata, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, type = "response", c.hat = 1,
parm.type = NULL, ...)
## S3 method for class 'AICunmarkedFitGDS'
modavgEffect(cand.set, modnames = NULL,
newdata, second.ord = TRUE, nobs = NULL, uncond.se = "revised",
conf.level = 0.95, type = "response", c.hat = 1,
parm.type = NULL, ...)
## S3 method for class 'AICunmarkedFitOccuFP'
modavgEffect(cand.set, modnames =
NULL, newdata, second.ord = TRUE, nobs = NULL, uncond.se =
"revised", conf.level = 0.95, type = "response",
c.hat = 1, parm.type = NULL, ...)
## S3 method for class 'AICunmarkedFitMPois'
modavgEffect(cand.set, modnames =
NULL, newdata, second.ord = TRUE, nobs = NULL, uncond.se =
"revised", conf.level = 0.95, type = "response",
c.hat = 1, parm.type = NULL, ...)
## S3 method for class 'AICunmarkedFitGMM'
modavgEffect(cand.set, modnames =
NULL, newdata, second.ord = TRUE, nobs = NULL, uncond.se =
"revised", conf.level = 0.95, type = "response",
c.hat = 1, parm.type = NULL, ...)
## S3 method for class 'AICunmarkedFitGPC'
modavgEffect(cand.set, modnames =
NULL, newdata, second.ord = TRUE, nobs = NULL, uncond.se =
"revised", conf.level = 0.95, type = "response",
c.hat = 1, parm.type = NULL, ...)

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 
newdata 
a data frame with two rows and where the columns correspond to the
explanatory variables specified in the candidate models. Note that this
data set must have the same structure as that of the original data frame
for which we want to make predictions, specifically, the same variable
type and names that appear in the original data set. Each row of the
data set defines one of the two groups compared. The first row in

second.ord 
logical. If 
nobs 
this argument allows the specification of a numeric value other than
total sample size to compute the AICc (i.e., 
uncond.se 
either, 
conf.level 
the confidence level (1  α) requested for the computation of unconditional confidence intervals. To obtain confidence intervals corrected for multiple comparisons between pairs of treatments, it is possible to adjust the α level according to various strategies such as the Bonferroni correction (Dunn 1961). 
type 
the scale of prediction requested, one of 
c.hat 
value of overdispersion parameter (i.e., variance inflation factor) such
as that obtained from 
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

... 
additional arguments passed to the function. 
Details
The strategy used here to compute effect sizes is to work from the
newdata
object to create two predictions from a given model and
compute the differences and standard errors between both values. This
step is executed for each model in the candidate model set, to obtain a
modelaveraged estimate of the effect size and unconditional standard
error. As a result, the newdata
argument is restricted to two
rows, each for a given prediction. To specify each group, the values
entered in the column for each explanatory variable can be identical,
except for the grouping variable. In such a case, the function will
identify the variable and the assign group names based on the values of
the variable. If more than a single variable has different values in
its respective column, the function will print generic names in the
output to identify the two groups. A sensible choice of value for the
explanatory variables to be held constant is the average of the
variable.
Modelaveraging effect sizes is most useful in true experiments (e.g., ANOVAtype designs), where one wants to obtain the best estimate of effect size given the support of each candidate model. This can be considered as a informationtheoretic analog of traditional multiple comparisons, except that the information contained in the entire model set is used instead of being restricted to a single model. See 'Examples' below for applications.
modavgEffect
calls the appropriate method depending on the class
of objects in the list. The current classes supported include
aov
, glm
, gls
, lm
, lme
, mer
,
glmerMod
, lmerMod
, rlm
, survreg
, as well as
unmarkedFitOccu
, unmarkedFitColExt
,
unmarkedFitOccuFP
, unmarkedFitOccuRN
,
unmarkedFitPCount
, unmarkedFitPCO
, unmarkedFitDS
,
unmarkedFitGDS
, unmarkedFitMPois
, unmarkedFitGMM
,
and unmarkedFitGPC
classes.
Value
The result is an object of class modavgEffect
with the following
components:
Group.variable 
the grouping variable defining the two groups compared 
Group1 
the first group considered in the comparison 
Group2 
the second group considered in the comparison 
Type 
the scale on which the modelaveraged effect size was computed (e.g., response or link) 
Mod.avg.table 
the full model selection table including the entire set of candidate models 
Mod.avg.eff 
the modelaveraged effect size based on the entire candidate model set 
Uncond.SE 
the unconditional standard error for the modelaveraged effect size 
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) Modelbased 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 informationtheoretic 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.
Burnham, K. P., Anderson, D. R., Huyvaert, K. P. (2011) AIC model selection and multimodel inference in behaviorial ecology: some background, observations and comparisons. Behavioral Ecology and Sociobiology 65, 23–25.
Dail, D., Madsen, L. (2011) Models for estimating abundance from repeated counts of an open population. Biometrics 67, 577–587.
Dunn, O. J. (1961) Multiple comparisons among means. Journal of the American Statistical Association 56, 52–64.
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. AmphibiaReptilia 27, 169–180.
Royle, J. A. (2004) Nmixture models for estimating population size from spatially replicated counts. Biometrics 60, 108–115.
See Also
AICc
, aictab
, c_hat
,
confset
, evidence
, importance
,
modavgShrink
, modavgPred
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264  ##heights (cm) of plants grown under two fertilizers, Ex. 9.5 from
##Zar (1984): Biostatistical Analysis. Prentice Hall: New Jersey.
heights < data.frame(Height = c(48.2, 54.6, 58.3, 47.8, 51.4, 52.0,
55.2, 49.1, 49.9, 52.6, 52.3, 57.4, 55.6, 53.2,
61.3, 58.0, 59.8, 54.8),
Fertilizer = c(rep("old", 10), rep("new", 8)))
##run linear model hypothesizing an effect of fertilizer
m1 < lm(Height ~ Fertilizer, data = heights)
##run null model (no effect of fertilizer)
m0 < lm(Height ~ 1, data = heights)
##assemble models in list
Cands < list(m1, m0)
Modnames < c("Fert", "null")
##compute model selection table to compare
##both hypotheses
aictab(cand.set = Cands, modnames = Modnames)
##note that model with fertilizer effect is much better supported
##than the null
##compute modelaveraged effect sizes: one model hypothesizes a
##difference of 0, whereas the other assumes a difference
##prepare newdata object from which differences between groups
##will be computed
##the first row of the newdata data.frame relates to the first group,
##whereas the second row corresponds to the second group
pred.data < data.frame(Fertilizer = c("new", "old"))
##compute best estimate of effect size accounting for model selection
##uncertainty
modavgEffect(cand.set = Cands, modnames = Modnames,
newdata = pred.data)
##classical oneway ANOVA typedesign
## Not run:
##generate data for two groups and control
set.seed(seed = 15)
y < round(c(rnorm(n = 15, mean = 10, sd = 5),
rnorm(n = 15, mean = 15, sd = 5),
rnorm(n = 15, mean = 12, sd = 5)), digits = 2)
##groups
group < c(rep("cont", 15), rep("trt1", 15), rep("trt2", 15))
##combine in data set
aov.data < data.frame(Y = y, Group = group)
rm(y, group)
##run model with group effect
lm.eff < lm(Y ~ Group, data = aov.data)
##null model
lm.0 < lm(Y ~ 1, data = aov.data)
##compare both models
Cands < list(lm.eff, lm.0)
Mods < c("group effect", "no group effect")
aictab(cand.set = Cands, modnames = Mods)
##model with group effect has most of the weight
##compute modelaveraged effect sizes
##trt1  control
modavgEffect(cand.set = Cands, modnames = Modnames,
newdata = data.frame(Group = c("trt1", "cont")))
##trt1 differs from cont
##trt2  control
modavgEffect(cand.set = Cands, modnames = Modnames,
newdata = data.frame(Group = c("trt2", "cont")))
##trt2 does not differ from cont
## End(Not run)
##twoway ANOVA type design, Ex. 13.1 (Zar 1984) of plasma calcium
##concentration (mg/100 ml) in birds as a function of sex and hormone
##treatment
## Not run:
birds < data.frame(Ca = c(16.87, 16.18, 17.12, 16.83, 17.19, 15.86,
14.92, 15.63, 15.24, 14.8, 19.07, 18.77, 17.63,
16.99, 18.04, 17.2, 17.64, 17.89, 16.78, 16.92,
32.45, 28.71, 34.65, 28.79, 24.46, 30.54, 32.41,
28.97, 28.46, 29.65),
Sex = c("M", "M", "M", "M", "M", "F", "F", "F", "F",
"F", "M", "M", "M", "M", "M", "F", "F", "F", "F",
"F", "M", "M", "M", "M", "M", "F", "F", "F", "F",
"F"),
Hormone = as.factor(c(1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3)))
##candidate models
##interactive effects
m.inter < lm(Ca ~ Sex + Hormone + Sex:Hormone, data = birds)
##additive effects
m.add < lm(Ca ~ Sex + Hormone, data = birds)
##Sex only
m.sex < lm(Ca ~ Sex, data = birds)
##Hormone only
m.horm < lm(Ca ~ Hormone, data = birds)
##null
m.0 < lm(Ca ~ 1, data = birds)
##model selection
Cands < list(m.inter, m.add, m.sex, m.horm, m.0)
Mods < c("interaction", "additive", "sex only", "horm only", "null")
aictab(Cands, Mods)
##there is some support for a hormone only treatment, but also for
##additive effects
##compute modelaveraged effects of sex, and set the other variable
##to a constant value
##M  F
sex.data < data.frame(Sex = c("M", "F"), Hormone = c("1", "1"))
modavgEffect(Cands, Mods, newdata = sex.data)
##no support for a sex main effect
##hormone 1  3, but set Sex to a constant value
horm1.data < data.frame(Sex = c("M", "M"), Hormone = c("1", "3"))
modavgEffect(Cands, Mods, newdata = horm1.data)
##hormone 2  3, but set Sex to a constant value
horm2.data < data.frame(Sex = c("M", "M"), Hormone = c("2", "3"))
modavgEffect(Cands, Mods, newdata = horm2.data)
## End(Not run)
##Poisson regression with anuran larvae example from Mazerolle (2006)
## Not run:
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,
family = poisson, offset = log(Effort),
data = min.trap)
Cand.mod[[2]] < glm(Num_anura ~ log.Perimeter, family = poisson,
offset = log(Effort), data = min.trap)
Cand.mod[[3]] < glm(Num_anura ~ Type, family = poisson,
offset = log(Effort), data = min.trap)
Cand.mod[[4]] < glm(Num_anura ~ 1, family = poisson,
offset = log(Effort), data = min.trap)
##check chat for global model
vif.hat < c_hat(Cand.mod[[1]]) #uses Pearson's chisquare/df
##assign names to each model
Modnames < c("type + logperim", "type", "logperim", "intercept only")
##compute modelaveraged estimate of difference between abundance at bog
##pond and upland pond
##create newdata object to make predictions
pred.data < data.frame(Type = c("BOG", "UPLAND"),
log.Perimeter = mean(min.trap$log.Perimeter),
Effort = mean(min.trap$Effort))
modavgEffect(Cand.mod, Modnames, newdata = pred.data, c.hat = vif.hat,
type = "response")
##little suport for a pond type effect
## End(Not run)
##mixed linear model example from ?nlme
## Not run:
library(nlme)
Cand.models < list( )
Cand.models[[1]] < lme(distance ~ age, data = Orthodont, method="ML")
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")
Modnames < c("age", "age + sex", "null", "sex")
data.other < data.frame(age = mean(Orthodont$age),
Sex = factor(c("Male", "Female")))
modavgEffect(cand.set = Cand.models, modnames = Modnames,
newdata = data.other, conf.level = 0.95, second.ord = TRUE,
nobs = NULL, uncond.se = "revised")
detach(package:nlme)
## End(Not run)
##site occupancy analysis example
## Not run:
library(unmarked)
##single season model
data(frogs)
pferUMF < unmarkedFrameOccu(pfer.bin)
##create a bogus site group
site.group < c(rep(1, times = nrow(pfer.bin)/2), rep(0, nrow(pfer.bin)/2))
## add some fake covariates for illustration
siteCovs(pferUMF) < data.frame(site.group, sitevar1 =
rnorm(numSites(pferUMF)),
sitevar2 = runif(numSites(pferUMF)))
## observation covariates are in sitemajor, observationminor order
obsCovs(pferUMF) < data.frame(obsvar1 =
rnorm(numSites(pferUMF) * obsNum(pferUMF)))
fm1 < occu(~ obsvar1 ~ site.group, pferUMF)
fm2 < occu(~ obsvar1 ~ 1, pferUMF)
Cand.mods < list(fm1, fm2)
Modnames < c("fm1", "fm2")
##model selection table
aictab(cand.set = Cand.mods, modnames = Modnames, second.ord = TRUE)
##modelaveraged effect sizes comparing site.group 1  site.group 0
newer.dat < data.frame(site.group = c(0, 1))
modavgEffect(cand.set = Cand.mods, modnames = Modnames, type = "response",
second.ord = TRUE, newdata = newer.dat, parm.type = "psi")
##no support for an effect of site group
## End(Not run)
##single season Nmixture models
## Not run:
data(mallard)
##this variable was created to illustrate the use of modavgEffect
##with detection variables
mallard.site$site.group < c(rep(1, 119), rep(0, 120))
mallardUMF < unmarkedFramePCount(mallard.y, siteCovs = mallard.site,
obsCovs = mallard.obs)
siteCovs(mallardUMF)
tmp.covs < obsCovs(mallardUMF)
obsCovs(mallardUMF)$date2 < tmp.covs$date^2
(fm.mall < pcount(~ site.group ~ length + elev + forest, mallardUMF, K=30))
(fm.mallb < pcount(~ 1 ~ length + elev + forest, mallardUMF, K=30))
Cands < list(fm.mall, fm.mallb)
Modnames < c("one", "null")
##model averaged effect size of site.group 1  site.group 0 on response
##scale (point estimate)
modavgEffect(Cands, Modnames, newdata = data.frame(site.group = c(0, 1)),
parm.type = "detect", type = "response")
##model averaged effect size of site.group 1  site.group 0 on link
##scale (here, logit link)
modavgEffect(Cands, Modnames, newdata = data.frame(site.group = c(0, 1)),
parm.type = "detect", type = "link")
detach(package:unmarked)
## End(Not run)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker. Vote for new features on Trello.