rmleoptim

Problem specifying arguments when using "optim" for MLE

I am trying to implement Maximum likelihood estimation for the following 2-parameter Beta-Poisson model

enter image description here

Working through other solutions in StackOverflow I came to the conclusion - possibly incorrectly - that the "optim" function might be the best solution. Accordingly, I have made an attempt at coding this up:

Log_Lik_BP <- function(alpha_beta, D, Y, N){
  alpha <- alpha_beta[1]
  beta <- alpha_beta[2]
  sum(Y * log(1 - (1 + (D/beta))^(-alpha))) - alpha*sum(N-Y)*log(1+(D/beta))
}

optim(par = c(0,0), fn = Log_Lik_BP, D = df$D, N = df$N, Y = df$Y)

But I haven't been very successful and I am receiving the following error message, with "length 22" referring to the 22 data points in the dataframe

Error in optim(par = c(0, 0), fn = Log_Lik_BP, D = df$D, N = df$N, Y = df$Y) : 
  objective function in optim evaluates to length 22 not 1

The dataset I am applying this MLE estimation is as follows:

D <- c(4.2,8,10.6,11.9,12.2,12.6,13,13.2,28.3,28.7,68.2,69,71.1,83,91.7,95.5,106.5,136.5,171.2,186.9,309.3,1557.1)
N <- c(1,1,1,1,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
Y <- c(0,0,1,1,2,1,0,1,0,1,1,0,1,1,1,1,1,1,1,1,1,1)
df <- data.frame(D, N, Y)

Thanks in advance for advice and corrections to my code

Solution

  • There are a bunch of issues here, I'm going to walk through the debugging process.

    When troubleshooting optimization problems, it always helps to start by evaluating the function at the starting parameters (e.g. Log_Lik_BP(c(0,0), D, N, Y)) to remove one layer of complexity. Furthermore, it helps to either (1) apply debug() or debugonce() to function to step through it, or (2) run the code by bits and pieces outside of the function scope to see what's going on.

    1. Run Log_Lik_BP(c(0,0), D, N, Y)). Result: a vector of 22 NaN values. The first thing we need to do is see why it's length 22 rather than 1 as it should be (we'll come back to the NaN problem later).

    2. Set alpha <- beta <- 0 and evaluate the individual components sum(Y * log(1 - (1 + (D/beta))^(-alpha))) and alpha*sum(N-Y)*log(1+(D/beta)). Result: the first component is length-1, as it should be, the second component is length 22. Aha! There's incorrect grouping in the formula, we need the parentheses to include the the log(1+(D/beta)) term as well, i.e. alpha*sum((N-Y)*log(1+(D/beta))).

    3. Fix Log_Lik_BP accordingly. Result: If we evaluate the function at c(0,0), we again get NaN; if we try optim(), we get "function cannot be evaluated at initial parameters".

    4. Looking at the formula, c(0,0) seems like a questionable starting point, since we need to divide by beta. What happens if we start at c(1,1) instead? Result: optim() runs without an error, but we get funny-looking answers - the default Nelder-Mead algorithm runs out of iterations, and it looks like beta is converging to zero.

    5. A slight digression, but I tried setting using an alternate method that allows for setting bounds (method = "L-BFGS-B") and setting lower bounds on alpha and beta (lower = c(0.001, 0.001)) (setting the lower bounds to exactly 0 runs into trouble when beta=0, stops with "L-BFGS-B needs finite values of 'fn'"). Result: this reaches a solution fairly quickly, but the result still seems wonky (alpha is huge and beta hits the boundary). (This step may be unnecessary, see next step for the real solution to the remaining problems.)

    6. Realize that optim() minimizes the objective function by default, and you have defined a log-likelihood function that should be maximized. Solution: set control = list(fnscale = -1) instead to maximize. Results: reasonable-looking parameter estimates and negative log-likelihood. You should still check that the results make sense in context!