Error with custom density function definition for mle2 formula call

I want to define my own density function to be used in the formula call to mle2 from R's bbmle package. The parameters of the model are estimated but I cannot apply functions like residuals or predict on the returned mle2 object.

This is an example in which I define a function for a simple Poisson model.

library(bbmle)

set.seed(1)
hpoisson <- rpois(1000, 10)

myf <- function(x, lambda, log = FALSE) {
  pmf <- (lambda^x)*exp(-lambda)/factorial(x)
  if (log)
    log(pmf)
  else
    pmf
}

myfit <- mle2(hpoisson ~ myf(lambda), start = list(lambda=9), data=data.frame(hpoisson))
residuals(myfit)

In myfit, the lambda is estimated correctly but when I call residuals on myfit, I get an error which says:

Error in myf(9.77598906811668) : 
  argument "lambda" is missing, with no default

On the other hand, if I simply fit the model as follows using R's built-in dpois function, the residuals are computed:

myfit <- mle2(hpoisson ~ dpois(lambda), start = list(lambda=9), data=data.frame(hpoisson))
    residuals(myfit)

Could anyone please tell me what I am doing wrong in the function definition of myf?

Thanks

Solution

It's not very clearly explained in the documents, but there are a few prerequisites for using custom density functions:

the function's name must start with d, must have first argument x, and must have a named argument log. (The log argument has to do something sensible: in particular, mle2 will call the function with log=TRUE, and the function had better return the log-likelihood!) In general, although it's not required, it's more numerically sensible to compute the log-likelihood directly and then exponentiate if log=FALSE, rather than computing the likelihood and logging it if log=TRUE (there are cases, such as zero-inflated models, where this isn't really feasible). For example, compare my dmyf() definition with the myf() definition in the OP's code ...
in order to use additional methods such as predict you have to define an additional function whose name starts with s; it returns a list of moments, summary statistics, etc. for a specified parameter -- see example below, which is copied from bbmle::spois.

library("bbmle")
set.seed(1)
hpoisson <- rpois(1000, 10)

dmyf <- function(x, lambda, log = FALSE) {
    logpmf <- x*log(lambda)-lambda-lfactorial(x)
    if (log) return(logpmf)  else return(exp(logpmf))
}
smyf <- function(lambda) {
    list(title = "modified Poisson",
         lambda = lambda, mean = lambda,
         median = qpois(0.5, lambda),
         mode = NA, variance = lambda, sd = sqrt(lambda))
}
myfit <- mle2(hpoisson ~ dmyf(lambda),
              start = list(lambda=9), data=data.frame(hpoisson))
residuals(myfit)