[R-sig-ME] [R] understanding I() in lmer formula

[Emmanuel Curis]

Replies in the text. Please read about maximum likelihood methods, the
objective function is the likelihood of the data, and write it down
will give you what you want. In practice, it might be slightly
different if you used restricted maximum likelihood (REML) instead of
maximum likelihood (ML), but that does not change the interpretation
of the different models and outputs, just the objective function.

On Thu, Jun 15, 2017 at 02:49:50PM +0000, Don Cohen wrote:
« 
«           (    V(A)     cov( A, B ) )   ( f   g )
«   Sigma = (                         ) = (       )
«           ( cov( A, B )     V(B)    )   ( g   h )
« 
« 
«   so the aim of lme4 (or any other software) is to find the << best >>
«   values for all of these five parameters.
« 
« So when I run lmer I should be able to recover these 5 values, right?
« How do I do that?

All of them are in summary( lmer( formula, data = ... ) )

For µA and µB estimtions, you can have them with fixef()
For f, g and h, you can have them with VarCorr()

« When you say it finds the best values, that means to me that the
« objective function (whatever we're trying to minimize) includes 
« something that depends on those parameters.
« What is that function? 

The opposite of the log likelihood, -log L, of the data, more
precisely assuming independant data points

L = \prod_{i=1}^n \prod_{,j=1}^{n_i} F'(Yij=y_ij|Ai=ai,Bi=bi) G'(Ai=ai, Bi=Bi)

where F' is the density of a Gaussian of mean ai + bi * days_{i,j} and
variance sigma (the variance of epsilon, the residual) and G' the
joint density of a binormal Gaussian vector of expectation (µA, µB)
and of covariance matrix Sigma = ( (f, g ), ( g, h ) ) and you
minimize -log L. i the subject index, j the measure index within the subject.

« I guess you have to pick some particular model in order to show its
« objective function, and I suggest the simplest possible model, in
« this case something like
«  reaction ~ 1+days + (1+days|subject)
« and similarly the || version, which I gather differs only in that
« some g terms are left out.
« 
«   This is done by trying to
«   find the values that allow to say << I obtained my data because they
«   were the more likely to occur >>. This leads to a complex function that
«   often reduce to least-squares, but not always, and in Gaussian
«   mixed-effects models are not linear least-squares because of the f, g
«   and h parameters.
« 
«   Traditionnally, you fix uA = uB = 0 because should they have other
«   values, they could not be distinguished from the fixed part of the
«   model (but you could also say that there is no distinct fixed part and
«   that you try to find their values, it's the same).
« 
«   When you use (1|Day), you allow lmer to fit a model with f, g and h
«   variable.
« 
«   When you use (1||Day), you constrain lmer to fit a model with g = 0
«   and only f and h can be fitted.
« 
« All of this makes sense to me, but in order to really understand what's
« going on I want to know the objective function.

Just write it down using the hints above... Or use the formulas in
Douglas Bates' book (or other books). They _are_ the objective
function. By e-mail, that would be almost impossible to write them in
a readable way.

«   Note that with lmer, f, h (and g if allowed) are not obtained by first
«   computing slope and intercept for all subjects, then doing usual
«   descriptive statistics ...
« 
« I understand, but however the software actually works, I should be able
« to see that the objective function is minimized by simply computing it
« on the output and then also on various perturbations of the output.

Once you wrote it, you can use dnorm to compute it, and use the logLik
function to check the results.

« That would tell me WHAT the software is doing.  I could worry later about
« HOW it does that.

It maximises the likelihood (L) — or, equivalently, minimises -log L —
as a function of (µA, µB, f, g, h).

«   Last note: do not confuse the correlation between A and B, the random
«   effects in the population, given by g, and the correlation between the
«   estimators of the (mean) slope [uA] and the estimator of the (mean)
«   intercept [uB], M_A and M_B, which may be what you had in mind when
«   saying that A and B values were correlated << before >> (it exists in
«   usual linear regression). They correspond to different ideas:
« 
«   g means that there is in the population a kind of constraint
«    between A and B ;
« 
«   the correlation between the _estimators_ means that any error on
«    the estimation of the (mean) slope will lead to an error in the
«    estimation of the (mean) intercept and it is a property of the
«    method used to find them, not of the data/underlying world.
« 
« I had not even thought about the second of those.  
« But I think that is similar to one of the outputs of summary(lmer(...))
« where it says correlation of fixed effects.
« 
«   Hope this clarifies,
« 
« So far I don't think I've learned anything new, but I may be getting
« close to describing to you what I'm missing.

-- 
                                Emmanuel CURIS
                                emmanuel.curis at parisdescartes.fr

Page WWW: http://emmanuel.curis.online.fr/index.html