[R-sig-ME] Teaching Mixed Effects

[Raubertas, Richard]

The accumulation of quoted posts in this thread is quite long, so I 
have trimmed out all but some paragraphs from Doug Bates that prompt my 
comments.

> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org 
> [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf 
> Of Douglas Bates
> Sent: Thursday, January 22, 2009 5:40 PM
> To: Bolker,Benjamin Michael
> Cc: R Models Mixed
> Subject: Re: [R-sig-ME] Teaching Mixed Effects
> 
   [ ... ]
> 
> If we want to perform a hypothesis test related to a fixed-effects
> term in a mixed model (and, for the moment, I will not go into the
> question of whether statistical inferences should always be phrased as
> the result of hypothesis tests) I would claim we should start at the
> beginning, which is considering two models for the data at hand, one
> model being a special case of the other.  We need to decide how we
> measure the quality of the fit of the general model relative to the
> more specific model and how we measure the additional cost of the
> general model.  Then we need to formulate a test statistic.  If we are
> incredibly lucky, the null distribution of this test statistic will be
> well-defined (that is, it will not depend on the values of other,
> unknown parameters) and we can evaluate probabilities associated with
> it.  That does happen in the case of the Gaussian linear model.  I
> personally don't think it will be possible to possible to provide a
> general approach that isolates the effect of a fixed-effect term in a
> linear mixed model using a statistic that does not depend on the
> values of other parameters.  I would be delighted if someone can do it
> but I think there is too much that goes right in the case of the
> Gaussian linear model to expect that the same incredible
> simplifications will apply to other models.
> 
> I don't feel that holy grail of inference in mixed effects models
> should be a magical formula for degrees of freedom to be associated
> with some ratio that looks like a t or an F statistic (despite the
> strongly held beliefs of those in the First Church of the
> Kenward-Roger Approximation). Certainly there has been a lot of
> statistical research related to approximating a difficult distribution
> by a more common distribution but I view this approach as belonging to
> an earlier era.  It is the approach embodied in software like SAS
> whose purpose often seems to be to evaluate difficult formulas and
> provide reams of output including every number that could possibly be
> of interest.  I think we should use the power of current and future
> computers to interact with and explore data and models for the data.
> MCMC is one way to do this.  In nonlinear regression Don and I
> advocated profiling the sum of squares function with respect to the
> values of individual parameters as another way of assessing the actual
> behavior of the model versus trying to formulate an approximation.
> I'm sure that creative people will come up with many other ways to use
> the power of computers to this end.  The point is to explore the
> actual behavior of the model/data combination, not to do just one fit,
> calculate a bunch of summary statistics, apply approximate
> distributions to get p-values and go home.

But a goal of this exploration of the model/data is still to come up 
with an approximation to a sampling distribution (of a test statistic 
or parameter estimate), right?  Ultimately, once all the exploration is 
done, the scientific researcher still wants to be able to tell his/her 
colleagues "Based on these data my conclusion about the effect of 
treatment X is ______ and my confidence in this conclusion is _____."  
That is, the researcher wants to make *inferences* from the data.

> 
> If we want to generalize methods of inference we should consider the
> whole chain of reasoning that leads us to the result rather than
> concentrating only on the last step, which is "now that I have the
> value of a statistic how do I convert it to a p-value?" or, even more
> specifically, "I have calculated something that looks like a t-ratio
> so I am going to assume that its distribution is indeed a
> t-distribution which leaves me with only one question and that is on
> how many degrees of freedom".
> 
> I appreciate that this is inconvenient to those applying such models
> to their data.  Philosophical discussions of the fundamentals of
> statistical inference are all well and good but when the referees on
> your paper say you have to provide a p-value for a particular term or
> tests, it is a practical matter, not an academic, theoretical debate.
> Those with sufficient energy and skill plus a stellar reputation as a
> researcher may be able to convince editors that p-values are not the
> "be all and end all" of data analysis - Reinhold Kleigl has been able
> to do this in some cases - but that is definitely not the path of
> least resistance.  The sad reality is that p-values have taken on a
> role as the coin of the realm in science that is far beyond what any
> statistician would have imagined.  (Apparently the default
> "significance level" of 5%, which is considered in some disciplines to
> be carved in stone, resulted from a casual comment by Fisher to the
> effect that he might regard an outcome that would be expected to occur
> less than, say, 1 time in 20 as "significant".)

This paragraph brings to mind some comments about p-values and
hypothesis 
tests that seem popular on this list.  Among users of lme4 a theme seems
to 
be that "this ignorant editor/referee is insisting that I demonstrate
that 
my discovery of the effect of treatment X is not a false positive.  How
can 
I get around this?"

Every field of science needs to protect its literature from being
overwhelmed 
by false "discoveries".  Since all the professional and economic rewards
go to 
those who make discoveries, there is enormous incentive to claim that
one 
has found an effect or association, and so there needs to be a reality
check.  
The conventional way of judging these claims is via p-values and
hypothesis 
tests.  Granting all of their faults and limitations, what do people
propose 
as a better way?  

Some will probably say we should focus on (interval) estimation of
parameters 
rather than testing.  But this doesn't solve the problem at hand.
Remember 
how this discussion started:  the difficulty of calculating reliable
p-values 
for fixed effects.  Because of the duality between hypothesis tests and
confidence 
intervals, if you can't get reliable p-values, then essentially by
definition 
you can't get reliable confidence intervals either.  So the issue of
testing-
versus-estimation is a red herring with respect to the deeper problem of

inference for fixed effects in mixed models.

Others may try to dodge the problem by claiming to be Bayesians,
calculating 
credible intervals from posterior distributions.  But this won't hold
water 
if they are using lme4 and mcmcsamp.  No honest Bayesian can use an 
"uninformative" prior (if such a thing even exists), or at least not
more 
than once in any area of research:  after analysis of one data set, the 
posterior from that analysis should inform the prior for the next, but
lme4 
has its priors hard-coded.  I think the real rationale for mcmcsamp is
the 
hope that it will produce results with good frequentist properties.  I
am 
not aware that this has been demonstrated for mixed models in the
peer-reviewed 
statistics literature.

> 
> It is unhelpful of me not to provide p-values in the lmer summaries
> but I develop the software out of interest in doing it as well as I
> possibly can and not because someone assigns me a task to compute
> something.  I really don't know of a good way of assigning p-values to
> t-ratios or F-ratios so I don't.  I still report the ratio of the
> estimate divided by it standard error, and even call it a t-ratio,
> because I think it is informative.
> 

Doug, can you elaborate on that last clause?  In what way is the
(absolute) 
ratio |T| informative that the monotonic transformation 2*(1 -
pnorm(abs(T)))
is not?  In other words, if a p-value (based in this case on a standard 
normal) is not reliable for inference, what inferential value does T
have?
Less formally, if T = 2, for example, what exactly do you conclude about

the parameter, and what is your confidence in that conclusion?
 
   [ ... ]
> 
> I feel that the likelihood ratio is a perfectly reasonable way of
> comparing two model fits where one is a special case of the other.  In
> fact, if the models have been fit by maximum likelihood, the
> likelihood ratio would, I think, be the first candidate for a test
> statistic.  The problem with likelihood ratio tests is not the
> likelihood ratio, per se -- it is converting the likelihood ratio to a
> p-value.  You need to be able to evaluate the distribution of the
> likelihood ratio under the null hypothesis.  The chi-square
> approximation to the distribution is exactly that - an approximation -
> and its validity depends on not testing at the boundary and on having
> a large sample, in some sense of the sample size.  If I were really
> interested in evaluating a p-value for the likelihood ratio I would
> probably try a parametric bootstrap to get a reference distribution.
> 

Yes, but the problem is that (as you noted earlier in your post) the
null 
distribution depends on the values of the nuisance parameters, except in

very special cases.  A simple parametric bootstrap conditions on the 
estimated values of those parameters, as if they were known.
Nevertheless, 
a function to do this might be a useful addition to lme4.

Rich Raubertas
Merck & Co.
Notice:  This e-mail message, together with any attachme...{{dropped:12}}