I know this is a kinda old question, but I just wonder if this has a solution now.
I usually perform the mixed effects model with lme4 package with lmer function. However, I know this function does not allow me to include the negative variance components in the model.
I would really want to include the negative variances in model with R. Does anyone have any suggestions that what packages I would use? Or, does the new lme4 allow it?
New lme4 doesn't allow it, nor does nlme. It looks like ASReml might do it, if you set the IGU argument as described here -- but ASReml is commercial, so you'd have to buy a license.
The more statistically sound way to deal with this situation is typically to fit a compound symmetry variance structure, which will allow negative correlations within groups. You can do this via nlme, or in several somewhat more experimental frameworks (e.g. the "flexLambda" development branch of lme4).
More discussion from the mailing list here.