Make lsmeans work when called inside a function with a nb.glm object

Question

I need to make a function (fit.function) that calls lsmeans with different formulas and data based on a negative binomial model fit from MASS (nb.glm).

I get the following error when I try to call lsmeans inside the function:

Error in terms.formula(formula, data = data) : 
'data' argument is of the wrong type
Error in ref.grid(object = list(coefficients = c(1.69377906086784, 
2.30790181649084,  : 
Perhaps a 'data' or 'params' argument is needed

It seems like the error has something to do with the environment of the ref.grid function.

Could anyone help me to fix the error? Any idea for a workaround?

My code is as follows:

library(lsmeans)
library(MASS)
df1 <-data.frame(y=rnbinom(100,size=0.75,mu =5 ), x="A")
df2 <-data.frame(y=rnbinom(100,size=0.75,mu =50 ), x="B")
df3 <-data.frame(y=rnbinom(100,size=0.75,mu =500 ), x="C")

df <- rbind(df1,df2,df3)

nb.fit<-function(formula,data){
  glm.nb(formula,data=data)
}

fit.function <- function(formula, data){
  lsmeans(glm.nb(formula, data = data), "x", adjust = "tuckey")
}

# lsmeans are calcultated when both lsmeans and glm.nb are explicitly called
main.fit <- lsmeans(glm.nb(y ~ x,data=df), "x", adjust = "tuckey")
main.fit

CLD <- cld(main.fit, type= "response")
plot(CLD)

# no problem wrapping glm.nb into nb.fit
class(glm.nb(y ~ x,df))
nb.model <-nb.fit(y ~ x,df)
class(nb.model)

# The Error appears once I wrap lsmeans into fit.function

func.fit <- fit.function(y ~ x,df)