I agree that Doug Bates has identified the proper course to follow for developing 
criteria to compare the performance of software  packages for parameter estimation 
in nonlinear mixed models. One should create a list of methods -- Laplace approximation, 
adaptive Gauss-Hermite, ...   where each method is well defined so that in principle 
at least any differences can be explained.

For example for the Laplace approximation the solution consists of the following "data":  

1.) The vectors x_hat and u_hat where u_hat is the value of the random effects 
    (the parameters we integrate over), x_hat is the maximizing value of the "other" 
    effects (I am not sure what the good generic term for these is as they can contain
    fixed effects and other parameters)
2.) L(x_hat,u_hat) i.e. the value of the log-likelihood at the maximum.
3.) L_uu(x_hat,u_hat), the Hessian wrt u.

A reliable algorithm must be able to calculate L_uu(x_hat,u_hat) with high 
accuracy especially when L_uu is almost singular. Such calculations must be done not 
only at x_hat but for other values of x since we must first find the value of x which 
maximizes the Laplace approximation using a nonlinear optimization method. 
It often happens that during his optimization, values of x will be chosen for 
which the calculation of the Hessian is numerically difficult. The entire 
optimization scheme can "hang" at this point if the most accurate methods 
are not used to calculate quantities of interest.

My point is that is it not sufficient for a package to return the correct values for 
say the Laplace approximation at the maximum they must also be able to find this 
maximum starting from different initial values of the parameters. So the set of 
starting values must be added to the specification of the problem.

Now a word from our sponsor

We have thought a lot about the problem of comparing different estimation software
for nonlinear mixed models because we produce and sell such software.
I believe that our software is the most numerically stable and computationally
efficient software for nonlinear mixed models. It has been difficult to
arrange for comparisons. Of course one can do the comparison oneself,
but there are two difficulties. How to convince people that the results are 
not "cooked" in some way.  In addition if one is not a recognized expert in the
other software it will always be claimed that the model simply was not
formulated optimally for that package -- and maybe it wasn't.

As an example of the performance of our software and the kind of difficulties
one encounters, we have compared our software with SAS NLMIXED
for a negative binomial mixed model. The results are for  Gauss-Hermite
integration. ithis was a good example for us because the
SAS NLMIXED results are published. However the authors of the paper
would proably not qualify as SAS NLMIXED expert as that in addition 
we obtained SAS code from Dale Mclerran. 

The SAS and AD Model Builder RE (the name of our software)
results were essentially identical. SAS shows some mild sensitivity to
the starting values -AD Model Builder does not. The SAS
runs take approximately 10-14 minutes. The AD Model Builder runs take
20-30 seconds i.e. about 25 times faster. 
uthe difficultyin comparison I refereed to above is that
SAS auotmaticlly adjust the number of points in the integregation. 
In AD MOdel builder we prefer to kkep that fixed to that we can
monitor the effect on the solution as the number of points is increased.
So the number of points in the integration is an additional thing 
which should  be specified in the comparison of different
packages.

Everything needed to reproduce our results can be found at

  http:// xxxxxxxxx