[R-sig-ME] [ADMB Users] NLMM Model Selection
bbolker at gmail.com
Tue Feb 8 14:38:30 CET 2011
-----BEGIN PGP SIGNED MESSAGE-----
On 11-02-08 01:45 AM, Chris Gast wrote:
> Hello Dr. Bolker,
> I hope you don't mind a personal email--my followup question is more
> relevant to your TREE paper than to the ADMB discussion. I'm happy to
> post this elsewhere (probably r-sig-mixed) if you prefer.
> Within your TREE article, you discuss using LRTs to test whether
> random effect terms should be included in the model. In the
> supplementary material for you TREE article, you provide a thorough
> dataset analysis, using lmer() to fit models and QAIC for model
> selection for random and then fixed effects. Why the difference?
> Also, I cannot seem to find any documentation for how the penalty for
> AIC is computed from lmer() (a difficult issue that you discuss in Box
> 3 of the TREE article). As you seem to be active in such discussions,
> perhaps you can provide some insight?
[I'm going to go ahead and post this back to r-sig-mixed-models for
public enlightenment (or correction, or scorn, as the case may be).]
Since we had overdispersion in the data, we chose to do a
quasi-likelihood analysis (which I would replace by an
observational-level random effect, if I were doing it again). Thus,
likelihood-based analyses in the strict sense were not available to us.
In that case we had to choose between (1) 'quasi-likelihood' F tests
(what anova.glm() does for quasi-likelihood models) or (2) QAIC. We went
for #2 (don't remember why).
As for AIC, looking deep inside the code at where the degrees of
freedom are computed for the log-likelihood (which is in turn used by
AIC to compute differences in numbers of parameters):
attr(val, "df") <-
dims[["p"]] + dims[["np"]] + as.logical(dims[["useSc"]])
This is "# fixed effect parameters" plus "# random effects parameters"
independent variances and covariances + "scale parameter" (residual
variance term, if used).
This is the "marginal AIC", which is more often the one I want
(rather than the conditional AIC), although it doesn't account (as
Greven and Kneib point out) for the boundedness of the
You could also, probably, work this out experimentally by comparing
logLik(m) and AIC(m) for various models ...
> Chris Gast
> cmgast at gmail.com
> On Sat, Feb 5, 2011 at 1:55 PM, Ben Bolker <bbolker at gmail.com> wrote:
> On 11-02-05 01:02 PM, Chris Gast wrote:
>>>> This isn't precisely an ADMB topic, but it seems as though ADMB users
>>>> might be knowledgeable in this regard.
>>>> I've searched the archives and haven't found a lot of discussion
>>>> regarding model selection in nonlinear mixed models. For a given
>>>> dataset, I have a series of models which differ in combinations of
>>>> structure, number of effects considered random, and assumed distribution
>>>> of random effect components, and would like some (preferably
>>>> likelihood-based) method to rank them. Burnham and Anderson (Model
>>>> Selection and Multimodel Inference, 2002, page 310) describe a method
>>>> based on shrinkage estimators where the penalty term is computed
>>>> somewhere between 1 and the number of random components, but this
>>>> appears to require both a single random effect and a fit of the model
>>>> where each random component is considered a parameter; neither of these
>>>> is feasible with my models (or, I suspect, many others). I can't simply
>>>> use LRTs to decide between a mixed model and its fixed counterpart,
>>>> because the value of interest for the sigma parameter lies on the
>>>> boundary of its space, 0.
> Although you could, approximately, by doubling the p value (for a
> single random effect, the null distribution of the deviance is a 50/50
> mixture of chi^2 with df=0 and df=1; this is equivalent to halving the
> area in the tail of the distribution or equivalently doubling the p
> value. (See references in Bolker et al 2009 TREE article.)
>>>> I have found some instances where the problem is basically ignored
>>>> (Hall, D.B. and Clutter, M. 2004. Multivariate multilevel nonlinear
>>>> mixed effects models for timer yield predictions. Biometrics, 60:16-24).
>>>> To quote: "...the first-order approximate log likelihood is treated as
>>>> the true log likelihood, and standard errors for parameter estimates,
>>>> likelihood ratio tests for nested models, and model selection criteria
>>>> such as AIC and BIC are formed in the usual way. Although the formal
>>>> justification of this
approach to inference
>>>> is an open problem, it is commonly used in practice, and we adopt it for
>>>> our purposes in this article."
>>>> One simple method would be to choose the model that best reconstructs
>>>> the original data as measured by the chi-squared test statistic
>>>> sum((O-E)^2/E), but again, it would be nice to have something
>>>> likelihood-based such that the framework is a cohesive, and the
>>>> principle of parsimony is in effect.
>>>> One additional question: these models also may include covariates.
>>>> Holding all other model features of a mixed-model constant, LRTs should
>>>> be justified for model selection of covariates only, as they result from
>>>> a mathematical restriction of some beta=0, correct? I see plenty of
>>>> information about the LASSO for covariate selection in NLMMs, but
>>>> haven't yet found the time to learn this technique.
> A quite technical but useful recent paper is:
> Greven, Sonja, and Thomas Kneib. 2010. On the Behaviour of Marginal
> and Conditional
> Akaike Information Criteria in Linear Mixed
> Models. Biometrika 97, no. 4: 773-789.
> There is a fundamental distinction between the 'marginal AIC' (for
> population-level predictions, i.e. where you want to predict future
> values for a different set of random effects than those measured) and
> the 'conditional AIC' (for group-level predictions where you want to
> predict future values for the same random effects measured); see
> <http://glmm.wikidot.com/faq> (recently updated) for more information.
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/
-----END PGP SIGNATURE-----
More information about the R-sig-mixed-models