Comparing a mixed model (subject as random effect) to a simple linear model (subject as a fixed effect) I am finishing up some analysis on a large set of data. I would like to take the linear model used in the first part of the work and re-fit it using an linear mixed model (LME). The LME would be very similar with the exception that one of the variables used in the model would be used as a random effect. This data comes from many observations (>1000) in a small group of subjects (~10) and I know that modeling the effect of subject is better done as a random effect (this is a variable that I want to shift). The R code would look like:
my_modelB <- lm(formula = A ~ B + C + D)    
lme_model <- lme(fixed=A ~ B + C, random=~1|D, data=my_data, method='REML')

Everything runs fine and the results are vastly similar. It would be nice if I could use something like RLRsim or an AIC/BIC to compare these two models and decide which is the most appropriate. My colleagues don't want to report the LME because there isn't an easily accessible way of choosing which is "better", even though I think the LME is the more appropriate model. Any suggestions?
 A: This is to add to @ocram's answer because it is too long to post as a comment. I would treat A ~ B + C as your null model so you can assess the statistical significance of a D-level random intercept in a nested model setup. As ocram pointed out, regularity conditions are violated when $H_0: \sigma^2 = 0$, and the likelihood ratio test statistic (LRT) will not necessarily be asymptotically distributed $\chi^2$. The solution was I taught was to bootstrap the LRT (whose bootstrap distribution will likely not be $\chi^2$) parametrically and compute a bootstrap p-value like this:
library(lme4)
my_modelB <- lm(formula = A ~ B + C)
lme_model <- lmer(y ~ B + C + (1|D), data=my_data, REML=F)
lrt.observed <- as.numeric(2*(logLik(lme_model) - logLik(my_modelB)))
nsim <- 999
lrt.sim <- numeric(nsim)
for (i in 1:nsim) {
    y <- unlist(simulate(mymodlB))
    nullmod <- lm(y ~ B + C)
    altmod <- lmer(y ~ B + C + (1|D), data=my_data, REML=F)
    lrt.sim[i] <- as.numeric(2*(logLik(altmod) - logLik(nullmod)))
}
mean(lrt.sim > lrt.observed) #pvalue

The proportion of bootstrapped LRTs more extreme that the observed LRT is the p-value. 
A: I am not totally sure to figure out what model is fitted when you use the lme function. (I guess the random effect is supposed to follow a normal distribution with zero mean?). However, the linear model is a special case of the mixed model when the variance of the random effect is zero. Although some technical difficulties exist (because $0$ is in the boundary of the parameter space for the variance) it should be possible to test $H_0:variance = 0$ vs $H_1: variance > 0$... 
EDIT
In order to avoid confusion: The test mentioned above is sometimes used to decide whether or not the random effect is significant... but not to decide whether or not it should be transformed into a fixed effect.
