roptimizationregressionmle

Using `bbmle:mle2` with vector parameters (already works using `optim`)

I am having some trouble using the bbmle:mle2 function when trying to do regression. To illustrate my problem, I have come up with a toy example.

We define the minus log-likelihood for the Poisson distribution (or any custom distribution):

LL <- function(beta, z, x){
  -sum(stats::dpois(x, lambda = exp(z %*% beta), log = TRUE))
}

In the code above, beta is parameter vector I would like to estimate, z is a model/design matrix and x is my variable of interest.

I then generate some random data to work with:

set.seed(2)
age <- round(exp(rnorm(5000, mean = 2.37, sd = 0.78) - 1))
claim <- rpois(5000, lambda = 0.07

I can easily use optim for my regression. Here is the intercept only model:

z1 <- model.matrix(claim ~ 1)
optim(par = 0, fn = LL, z = z1, x = claim)

Here is the intercept + age model:

z2 <- model.matrix(claim ~ age)
optim(par = c(0, 0), fn = LL, z = z2, x = claim)

The way a large number of different models can be assessed is quite easy, one just has to specify the model matrix. How can this be made to work with the mle2 function from the bbmle package?

I can do it, if beta is one-dimensional:

mle2(minuslogl = function(beta){ LL(beta = beta, z = z1, x = claim) },
  start = list(beta = 0))

But if beta is a vector, then I run into problems:

mle2(
  minuslogl = function(beta){ LL(beta = beta, z = z2, x = claim) },
  start = list(beta = c(0, 0)),
  vecpar = T,
  parnames = colnames(z2)
  )

I cannot get the syntax right and I cannot find any examples in the documentation or vignettes to help me. The problem is surely that beta is a vector now. The documentation suggests using the vecpar = T argument is the way forward for "compatibility with optim". Any tips would be appreciated.

Also, is there a way to pass the z and x arguments in my log-likelihood function in a more elegant way to mle2 like I have done so in optim?

Solution

  • I think the main problem is that you need to provide start as an atomic vector (not a list).

     library(bbmle)
 LL2 <- function(beta) {
     LL(beta, z = z2, x = claim)
 }
 parnames(LL2) <- colnames(z2)
 mle2(
   minuslogl = LL2 ,
     start = setNames(c(0,0),colnames(z2)),
     vecpar = TRUE
  )

    It might help to know that you can implement something like Poisson regression much more easily in bbmle with the formula interface and the parameters argument:

    mle2(claim~dpois(exp(loglambda)),     ## use log link/exp inverse-link
     data=data.frame(claim,age),      ## need to specify as data frame
     parameters=list(loglambda~age),  ## linear model for loglambda
     start=list(loglambda=0))         ## start values for *intercept*