Hi Doug & Ted, The multivariate Aitken accelerator suggested by Ted is numerically ill-conditioned. I have written a globally-convergent, general-purpose EM accelerator that works well. It is quite simple to implement for any EM-type algorithm (e.g. ECM, ECME which are all monotone in likelihood). My paper on that should be coming out soon in Scandinavian J of Stats. I would be interested in helping with its implementation for EM acceleration in large data sets with non-nested random effects.
Best, Ravi. ---------------------------------------------------------------------------- ------- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html ---------------------------------------------------------------------------- -------- -----Original Message----- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Ted Harding Sent: Monday, February 11, 2008 11:19 AM To: r-help@r-project.org Subject: Re: [R] [OT] good reference for mixed models and EM algorithm On 11-Feb-08 15:07:37, Douglas Bates wrote: > [...] > Except that Doug Bates doesn't use the EM algorithm for fitting mixed > models any more. The lme4 package previously had an option for > starting with EM (actually ECME, which is a variant of EM) iterations > but I have since removed it. For large data sets and especially for > models with non-nested random effects, the EM iterations just slowed > things down relative to direct optimisation of the log-likelihood. > [...] The "raw" EM Algorithm can be slow. I have had good success using Aitken Acceleration for it. The basic principle is that, once an interative algorithm gets to a stage where (approximately) [A] (X[n+1] - X) = k*(X[n] -X) where X[n] is the result at the n-th iteration, -1 < k < 1, and X is the limit, then you can use recent results to predict the limit. Taking the above equation literally, along with its analogue for the next step: [B] (X[n+2] - X) = k*(X[n+1] -X) from which k = (X[n+2] - X[[n+1])/(X[n+1] - X[n]) and then [C] X = (X[n+1] - X[n])/(1 - k). If X is multidimensional (say dimension = p), then k is a pxp matrix, and you want all its eigenvalues to be less than 1 in modulus. Then you use the matrix analogues of the above equations, based it on (p+1) successive iterations X[n], X[n+1], ... , X[n+p+1]), i.e. on the p-vector c(X[n+1]-X[n], X[n+1]-X[n+1], ... , X[n+p+1]-X[n+p]) I have had good experience with this too! The best method of proceeding is: Stage 1: Monitor the sequence {X[n]} until it seems that equation [A] is beginning to be approximately true; Stage 2: Apply equations [A], [B], [C] to estimate X. Stage 3: Starting at this X, run a few more iterations so that you get a better (later) estimate of k, and then apply [A], [B], [C] aqain to re-estimate X. Repeat stage 3 until happy (or bored). The EM Algorithm, in most cases, falls into the class of procedures to which Aitken Acceleration is applicable. Best wishes to all, Ted. -------------------------------------------------------------------- E-Mail: (Ted Harding) <[EMAIL PROTECTED]> Fax-to-email: +44 (0)870 094 0861 Date: 11-Feb-08 Time: 16:18:45 ------------------------------ XFMail ------------------------------ ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
