Standard error estimates for model averaged parameter estimates

[Terri]

HI Richard,

I am learning unmarked for the single season occupancy model case and
wanted to ask a couple of questions about model averaging. Bill
Deluca has been helping with code, but I wanted to double check to
make sure I understand the calculations fully.

We have run several models and used the fitList function and modSel
function to do model comparison:

fits <- fitList(Null = m1, m2 = m2, m3 = m3, m4 = m4)
(ms.wp <- modSel(fits, nullmod='Null'))

To get the model averaged parameter estimates, we extract the
information and coefficients from the models using

allInfo <- ms.wp@Full #extracts info from model selection object
coefs <- coef(ms.wp) #extracts coefficient estimates
coefs[is.na(coefs)] <- 0 #sets parameters to 0 that were not included
in a model.
wt <- allInfo$AICwt #extracts AIC wt from allInfo
colSums(coefs * wt) #multiplies coefs by wt and summarizes by
column

I just want to make sure I understand exactly what this computation is
doing. If a model does not contain a parameter (for example,
elevation), this code sets that estimate to 0. Then each model in the
model set contains an estimate for the covariate....those models that
did not contain the covariate were set to 0. Then each estimate is
multiplied by the model weight to provide the model averaged
estimate? This makes sense to me but I wanted to clarify because
there seems to be more than one way to model average. For instance,
some folks will include only those models where the covariate was
included in the model, and then rescale the weights so that they add
to 1 (so that models that do not include the covariate are not used in
the model averaging equation).

Also, can you help me with the code used to estimate the standard
errors for a model averaged estimate? Thank you for your
consideration.

[Matthew Giovanni]

Hi Terri.  I don't know if the functions are compatible with models in unmarked (I sort of doubt it) but I would highly recommend Marc Mazerolle's AICcmodavg package for R (http://cran.r-project.org/web/packages/AICcmodavg/index.html).  It is very handy for use with glm's and glmm's.

My impression is that there are no distinct rules regarding whether all models (inclusive of non-nested models) are used when model-averaging coefficients.  There's certainly plenty of literature where authors model-averaged with best subsets of models or only models containing the variable of interest, but some folks also consider that "cheating" in the context of a priori modeling.  Todd Arnold's 2010 JWM paper on model selection and assessment is a good supplement to Burnham and Anderson's book and existing knowledge on IT methods.

Matt

                                                                 
Matt Giovanni, Ph.D.
Research Wildlife Ecologist
402-617-3764
Research website:
http://sites.google.com/site/matthewgiovanni/

[Richard Chandler]

Hi Terri,

unmarked allows you to model-average *predictions* but not coefficients. See ?fitList for an example. Model-averaging regression coefficients is not implemented because there many issues with it. B&A discuss a few cases where it is not appropriate, and there are probably others. See sections 4.2.1 and 4.2.2 of B&A for more discussion. In spite of those issues, I will throw this out there because it might be useful in some cases.


modAvgCoefs <- function(ms) {

    theta.hat <- coef(ms)
    vars <- SE(ms)^2
    wts <- ms@Full$AICwt

    theta.hat[is.na(theta.hat)] <- vars[is.na(vars)] <- 0
    theta.hatbar <- colSums(theta.hat*wts)

    theta.hatbar.rk <- matrix(theta.hatbar, nrow(theta.hat), ncol(theta.hat),
        byrow=TRUE)
    theta.hatbar.se <- colSums(wts*sqrt(vars + (theta.hat-theta.hatbar.rk)^2))
    return(data.frame(Ests=theta.hatbar, SE=theta.hatbar.se))
    }
 

# This should return model-averaged coefficients and SEs
modAvgCoefs(ms.wp)


Let me know if works.

And following up on Matt's comment, Marc Mazerolle suggested to me that he was going to update his package to work with unmarked, but I don't know the status of that.


Richard