best subset selection on random effects model

I don't know of any package that will do it (or if violating the
marginality principle by having non-nested models even makes sense)
but you could always build your own search through all possible
models. Problem is for a large number of fixed effects this can take a
good bit of time...
Here is a piece of code I had from a few of years back. There may be
more elegant solutions that will reduce run time such as using
update.lme. Basically my code finds all possible combinations of model
terms with ?combn and than fits them.

require(nlme)
fm <- lme(distance ~age*Sex, random = ~ 1|Subject,data=Orthodont,method='ML')
Terms <- terms(fm)
todrop <-  1:length(attr(Terms,'term.labels'))
subs <- unlist(sapply(todrop,function(p)
combn(todrop,p,simplify=F)),recursive =F)
fm.subList <- lapply(subs[-length(subs)],function(s,...){
  newf<- formula(drop.terms(terms(fm),s,keep.response = TRUE))
  update(fm,newf)
})
names(fm.subList) <- sapply(fm.subList, function(x) paste('fm',attr(
terms(x),'term.labels'),sep='.'))
sort(sapply(fm.subList,BIC))               # 1st is best based on BIC

Cheers