An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20120918/80fb85b8/attachment.pl>
calculating number of parameters for AICc
3 messages · Jude Phillips, Paul Johnson, Ben Bolker
Might be of interest... Heredity (2012) 109, 235-245; doi:10.1038/hdy.2012.35; published online 18 July 2012 Bayesian adaptive Markov chain Monte Carlo estimation of genetic parameters Mathew et al. "a new fast adaptive Markov chain Monte Carlo (MCMC) sampling algorithm for the estimation of genetic parameters in the linear mixed model with several random effects." Cheers, Paul The University of Glasgow, charity number SC004401
Jude Phillips <birdlists at ...> writes:
Hi Listers, I am wondering about how 'K' (number of parameters in AICc equation) is calculated for mixed effects models. I am using aictab in the AICcmodavg library, to produce an AICc table for a set of models. In the table, K seems to equal number of fixed effects + number of random effects + 2. Is that the correct calculation? In Vaida and Blanchard (2005), they suggest that for the 'population' focus (ie when you are interested in the fixed effects, not the random effects, which is the case for my study), K is the number of fixed parameters counting mean parameters and variance components. I'm not sure how you calculate the number of variance components, so I'm not sure if this gives the same value for K as is supplied in aictab? So I'm wondering if the K given in aictab is really correct for a 'population focus', and if not, how should K be calculated?
Well, you can just look in the AICc.mer() function:
K <- attr(logLik(mod), "df")
this leads you (obscurely, I admit) to getMethods("logLik",sig="mer"):
attr(val, "df") <- dims[["p"]] + dims[["np"]] +
as.logical(dims[["useSc"]])
To get past here, you would have to dig quite a bit deeper, but
basically the answer is that dims[["p"]] is the same as
length(fixef(model)) -- i.e., the number of fixed-effect coefficients --
and dims[["np"]] is the same as length(getME(model,"theta")) -- the
number of variance parameters. For each random term in the model with
q components, it has q*(q+1)/2 parameters -- for example, a term of
the form (slope|group) has 3 parameters (intercept variance, slope
variance, correlation between intercept and slope). The last term
says whether the model uses a scale parameter or not (yes for
linear mixed models, no for typical GLMMs like binomial or Poisson).
Your statement of "number of fixed effects + number of random effects + 2"
doesn't seem correct, but perhaps if you gave an example ...
Whether to add nuisance parameters or not, such as the residual
variance parameter that is estimated based on the residual variance,
is as far as I know an open question. In the classic AIC context
it doesn't matter as long as one is consistent. In the AICc context,
I don't think anyone really knows the answer ... adding +1 for
the residual variance parameter (as lme4 does) would make the
model selection process slightly more conservative.