Well, lm() produces an OLS solution, which are also MLE solutions for the fixed
effects. I think this is an easy way, although maybe not the best.
BHHH is a numerical approximation that can be used when a closed form solution is not
available. It is less sophisticated than Newton-Raphson.
Is t
It is not obvious to me what parameters in what model you want to fit.
Function "optim" does very well with many different kinds of problems.
If you just want to estimate parameters of a probability distribution,
function "fitdistr" in library(MASS) will do that. A couple of days
ago, I needed
