Credit of this answer goes to @Joshua who gave an awesome answer when I posted this question to the R and Statistics community on Google+. I am simply pasting his answer below.
For running regression (without latent variable modeling), please read my notes typed after the quoted text.
Handling missing data with Maximum Likelihood on all available data
(so-called FIML) is a very useful technique. However, there are a
number of complications that make it challenging to implement in a
general way. Consider a simple linear regression model, predicting
some continuous outcome from say age, sex, and occupation type. In
OLS, you do not worry about the distribution of age, sex, and
occupation, only the outcome. Typically for categorical predictors,
they are dummy coded (0/1). To use ML, distributional assumptions are
required for all variables with missingness. By far the easiest
approach is multivariate normal (MVN). This is what for example Mplus
will do by default if you do not go out for your way to declare the
type of variable (e.g., categorical). In the simple example I gave,
you would probably want to assume, normal for age, Bernoulli for sex,
and multinomal for job type. The latter is tricky because what you
actually have are several binary variables, but you do not want to
treat them as Bernoulli. This means you do not want to work with the
dummy coded variables, you need to work with the actual categorical
variable so the ML estimators can properly use a multinomial, but this
in turn means that the dummy coding process needs to be built into the
model, not the data. Again complicating life. Further, the joint
distribution of continuous and categorical variables is nontrivial to
compute (when I run into problems like this in Mplus, it pretty
quickly starts to break down and struggle). Finally, you really
ideally specify the missing data mechanism. In SEM style, FIML, all
variables are essentially conditioned on all others, but this is not
necessarily correct. For example, perhaps age is missing as a
function not of gender and occupation type, but their interaction.
The interaction may not be important for the focal outcome, but if it
is important for missingness on age, then it must also be in the
model, not necessarily the substantive model of interest but the
missing data model.
lavaan will use ML for MVN, but presently I believe the categorical
data options are limited (again coming from the SEM field, this is
standard). Multiple imputation seems less elegant at first because it
makes explicit many hidden assumptions behind FIML (like
distributional assumptions for every variable and the predictive model
assumed for missingness on every variable). However, it gives you a
lot of control and explicitly thinking about the distribution of each
variable, and the optimal missing data mechanism for each is valuable.
I am becoming more and more convinced that Bayesian models are the way
to handle missing data. The reason is that they are very flexible at
including distributions for each variable, allowing many different
types of distributions, and can easily incorporate the variability
introduced by missing data on predictors, into the overall model
estimates (which is the trick with multiple imputation where you then
have to somehow combine results). Of course, these methods are not
the easiest and can take a lot of training and time to use.
So that doesn't really answer your question, but explains a bit of why
completely general frameworks for dealing with missingness are tricky.
In my semutils package for the covariance matrices, I use lavaan
underneath to use ML. I do that because I assume for a variance
covariance matrix that you are using continuous variables anyway so
that I assume my users are already assuming MVN for their data.
This means that if all variables with missingness are continuous, lavaan, a structural equation modelling (SEM) package is a nice one to use for FIML in R.
Now going back to my initial question. My intention was to have a magic fix for missingness when running linear regression. All my variables with missing were nice and continuous. So I proceeded to run my analyses in two styles:
- The usual way with multiple imputation
- In SEM style with lavaan using FIML.
I was missing a lot of things by doing regression in SEM style. Both styles gave similar coefficients and R squares, but in SEM style I didn't get the significance testing of the regression (the typical F values with df), instead I got fit indices that were not helpful as I had used up all my degrees of freedom. Also when one model had a larger R2 than another, I couldn’t find a way to compare whether the difference was significant. Additionally, doing regression the usual way gives access to a bunch of testing for regression assumptions that are invaluable. For a more detailed answer on this issue see my other question that was nicely answered by @StasK.
So the conclusion seems to be that lavaan is a decent package for FIML in R, yet the use of FIML depends on statistical assumptions and the type of analysis one is conducting. As far as regression (without latent variable modeling) goes, keeping it out of SEM programs and using multiple imputation is probably a wise move.