The R package nlme will let you do a simple comparison with the anova command, but will report the wrong p-values.
The reason is that when testing for the inclusion of a random effect, you are testing if some variance is equal to zero, i.e $H_0:\sigma_A^2=0$ which is on the boundary of the parameter space (as variances must be larger than zero). Asymptotic theory for LR testing requires the hypothesized value to be on the inside of the parameter space.
The correct p-value is computed from a mixture of two $\chi^2$-distributions. When including only a random intercept the p-value is computed from a 50/50 mixture of the $\chi^2_0$ (all weight on zero) and $\chi^2_1$ distributions.
In your case, this code will do what you want:
library(nlme)
fit1 = gls(mpg ~ wt, data=mtcars,method="REML")
fit2 = lme(mpg ~ wt,random=~1|cyl,data=mtcars,method="REML" )
> anova(fit1,fit2)
Model df AIC BIC logLik Test L.Ratio p-value
fit1 1 3 164.8199 169.0235 -79.40996
fit2 2 4 160.7230 166.3278 -76.36149 1 vs 2 6.096934 0.0135
I use the gls command to estimate the normal linear model using REML, which makes the models comparable by LR.
Like I said though, the p-value is incorrect. A correct p-value is found by:
p_value = 0.5*(pchisq(6.096934,df=0,lower.tail=F)+pchisq(6.096934,df=1,lower.tail=F))
> p_value
[1] 0.006770832
The low p-value indicating that the inclusion is worthwhile, the variance in question is larger than zero.
For a general inclusion of mixed effects, the mixture of distributions will change, but I'm sure there are packages that will do this automatically for you.