Bioequivalence and Bioavailability Forum

❝ Well, that’s extreme. Did I really give such an example somewhere?

❝ Up to you. As long as you describe the procedure unambiguously in the SAP (which has to be accepted by the agency), fine. That’s my standard procedure for ages. No guideline is that specific.

❝ Are you talking to me?

❝ Tricky – IMHO, case by case (should be laid down in an SOP or the SAP, of course). IIRC, Health Canada had a rule that 1 (one!) concentration is sufficient for C_max and 3 (oh dear!) fo AUC. Gone with the wind. THX, HC.

❝ The EMA’s Method B evaluated by lmer() of package lmerTest. Kenward-Roger’s degrees of freedom ≥ Satterthwaite’s.

Users are often surprised and alarmed that the summary of a linear mixed model fit by lmer provides estimates of the fixed-effects para­meters, standard errors for these parameters and a t-ratio but no p-values. Similarly the output from anova applied to a single lmer model provides the sequential sums of squares for the terms in the fixed-effects specification and the corresponding numerator degrees of freedom but no denominator degrees of freedom and, again, no p-values.

Because they feel that the denominator degrees of freedom and the corresponding p-values can easily be calculated they conclude that failure to do this is a sign of inattention or, worse, incompetence on the part of the person who wrote lmer (i.e. me).

Perhaps I can try again to explain why I don’t quote p-values or, more to the point, why I do not take the “obviously correct” approach of attempting to reproduce the results provided by SAS. Let me just say that, although there are those who feel that the purpose of the R Project – indeed the purpose of any statistical computing whatsoever – is to reproduce the p-values provided by SAS, I am not a member of that group. If those people feel that I am a heretic for even suggesting that a p-value provided by SAS could be other than absolute truth and that I should be made to suffer a slow, painful death by being burned at the stake for my heresy, then I suppose that we will be able to look forward to an exciting finale to the conference dinner at UseR!2006 next month. (Well, I won’t be looking forward to such a finale but the rest of you can.)

As most of you know the t-statistic for a coefficient in the fixed-effects model matrix is the square root of an F statistic with 1 numerator degree of freedom so we can, without loss of generality, concentrate on the F statistics that were present in the anova output. Those who long ago took courses in “analysis of variance” or “experimental design” that concentrated on designs for agricultural experiments would have learned methods for estimating variance components based on observed and expected mean squares and methods of testing based on “error strata”. (If you weren’t forced to learn this, consider yourself lucky.) It is therefore natural to expect that the F statistics created from an lmer model (and also those created by SAS PROC MIXED) are based on error strata but that is not the case.

The parameter estimates calculated by lmer are the maximum likelihood or the REML (residual maximum likelihood) estimates and they are not based on observed and expected mean squares or on error strata. And that’s a good thing because lmer can handle unbalanced designs with multiple nested or fully crossed or partially crossed grouping factors for the random effects. This is important for analyzing data from large observational studies such as occur in psychometrics.

There are many aspects of the formulation of the model and the calculation of the parameter estimates that are very interesting to me and have occupied my attention for several years but let’s assume that the model has been specified, the data given and the parameter estimates obtained. How are the F statistics calculated? The sums of squares and degrees of freedom for the numerators are calculated as in a linear model. There is a slot in an lmer model that is similar to the “effects” component in a lm model and that, along with the “assign” attribute for the model matrix provides the numerator of the F ratio. The denominator is the penalized residual sum of squares divided by the REML degrees of freedom, which is n-p where n is the number of observations and p is the column rank of the model matrix for the fixed effects.

Now read that last sentence again and pay particular attention to the word “the” in the phrase “the penalized residual sum of squares”. All the F ratios use the same denominator. Let me repeat that – all the F ratios use the same denominator. This is why I have a problem with the assumption (sometimes stated as more that just an assumption – something on the order of “absolute truth” again) that the reference distribution for these F statistics should be an F distribution with a known numerator degrees of freedom but a variable denominator degrees of freedom and we can answer the question of how to calculate a p-value by coming up with a formula to assign different denominator degrees of freedom for each test. The denominator doesn’t change. Why should the degrees of freedom for the denominator change?

Most of the research on tests for the fixed-effects specification in a mixed model begin with the assumption that these statistics will have an F distribution with a known numerator degrees of freedom and the only purpose of the research is to decide how to obtain an approximate denominator degrees of freedom. I don’t agree.

There is one approach that I think may be fruitful and that I am currently pursuing. The penalized least squares formulation of a mixed-effects model shows that the residual sum of squares in actually a penalized residual sum of squares and there is a quantity that behaves like the degrees of freedom for a penalized least squares problem. I will insert code to calculate this and see how that behaves in simulation.

For the time being, I would recommend using a Markov Chain Monte Carlo sample (function mcmcsamp) to evaluate the properties of individual coefficients (use HPDinterval or just summary from the coda package). Evaluating entire terms is more difficult but you can always calculate the F ratio and put a lower bound on the denominator degrees of freedom.

If anyone wants to contribute code to calculate the “obviously correct” denominator degrees of freedom from SAS I will incorporate it. However, be warned that the penalized least squares approach and sparse matrix methods used in lmer will require considerable translation from the formulas which typically occur in papers on this topic. Generally those formulas involve the inverse of an n by n matrix where n is the number of observations and you really, really don’t want to try to do that when you have a couple of million observations.