[R-sig-ME] Numerical methods used to compute correlation coefficients

[Douglas Bates]

On Wed, Apr 15, 2009 at 1:14 PM, H c <harlancampbell at gmail.com> wrote:

> Thank you.  This  is exactly what I wanted to know.  Perhaps one last
> question.  Is there any particular reason to prefer  this optimization
> routine over another, such as "optimize()"?

optimize() only handles scalar problems and you must define an
interval in which the optimum is to be found.

I have used other unconstrained optimizers such as nlm() and optim()
but had trouble with both of those.  The original version of the nlme
package for S and S-PLUS used the ms() optimizer in S which was based
on some Fortran code by David Gay.  Initially that code could not be
used in R because it was not released under an open source license.
Eventually it was released to the public domain and I wrote the nlminb
interface code for R, primarily so I could use it as the optimizer in
the nlme and lme4 packages.

Even this code by David Gay is far from ideal.  The development
version of the lme4 package factors out the evaluation of the deviance
function or the REML criterion from the optimization of these, so that
those who are interested can experiment with different optimizers and
see what works best in their problems.

> On Wed, Apr 15, 2009 at 12:59 PM, Douglas Bates <bates at stat.wisc.edu> wrote:
>>
>> On Tue, Apr 14, 2009 at 12:02 PM, H c <harlancampbell at gmail.com> wrote:
>> > Thank you for the quick response.
>> > when I refer to "correlation parameters", I mean the "generally small
>> > set of
>> > parameters \lambda" that parametrize the \Lambda_{i} Variance-Covariance
>> > matrix.
>> > For example, one has time series data such that every subject has been
>> > observed at 4 time points.  One wishes to model this using an AR(1)
>> > correlation structure within the mixed model.
>> > the AR(1) is parametrized by a fixed parameter, \phi :
>> > lme(y~X, random=~1|ID, method="ML", data=data,
>> > correlation=corAR1(0.5,form=~X,fixed=FALSE))
>>
>> > Since there is no closed form solution for the maximum-likelihood
>> > estimate
>> > of \phi.  what numerical methods are used to arrive at the given
>> > estimate?
>> > Hopefully this has clarified my question.
>>
>> In the sense that I know what you mean by the correlation parameters.
>> I'm not sure how to characterize the numerical methods other than to
>> say that the deviance (negative twice the log-likelihood) or the
>> corresponding version of the REML criterion is expressed with respect
>> to an unconstrained parameter and the resulting function optimized
>> using optimization software within R (nlminb).  You can check the
>> definition of the corAR family to determine exactly how the
>> parameterization is defined.
>>
>> > Thanks again,
>> > Harlan
>> >
>> > On Tue, Apr 14, 2009 at 12:37 PM, Douglas Bates <bates at stat.wisc.edu>
>> > wrote:
>> >>
>> >> On Tue, Apr 14, 2009 at 9:20 AM, H c <harlancampbell at gmail.com> wrote:
>> >> > I have have already posted the following question:  What numerical
>> >> > methods
>> >> > are used in nlme to estimate correlation parameters?
>> >> > I was referred to the Pinheiro and Bates book.  Unfortunately, on p.
>> >> > 202,
>> >> > section 5.1.1, under the title "Estimation and Computational
>> >> > Methods",
>> >> > no
>> >> > description on a numerical method is provided.  (When the data is
>> >> > transformed to work with the profiled likelihood(y->ystar), one needs
>> >> > the
>> >> > parameters that define Lambda.  How are these parameters estimated?)
>> >>
>> >> I'm not sure what you mean by "correlation parameters".  If you mean
>> >> the correlation parameters in the unconditional distribution of the
>> >> random effects then those are estimated by maximum likelihood (ML) or
>> >> residual maximum likelihood (REML).  The profiled deviance or the
>> >> profiled REML criterion is evaluated with respect to a transformed set
>> >> of parameters and this value is optimized.
>> >
>> >
>
>