Maximum Likelihood in R for a log function

I am having some issues on using the mle function in R. The model I have is, log(Y)~log(K)+log(L), and when I input this model into R using I keep on getting error message about missing the function minuslog1. How do I resolve this issue using the model I have listed above?

Below is the code and a small set of data.

Thank you.

> require(stats4)
Loading required package: stats4
> prod.mle<-mle(log(Y)~log(K)+log(L),)  # log version
Error in minuslogl() : could not find function "minuslogl"
In addition: Warning messages:
1: In formals(fun) : argument is not a function
2: In formals(fun) : argument is not a function

DF <- structure(list(Y = c(26971.71, 330252.5, 127345.3, 3626843, 37192.73
), K = c(32.46371, 28.42238, 5.199048, 327.807, 16.01538), L = c(3013256.014, 
135261574.9, 39168414.92, 1118363069, 9621912.503)), 
class = "data.frame", row.names = c(NA, -5L))

Solution

stats4::mle does not have a formula interface. You have to add a function minuslogl which is to be used to "to calculate negative log-likelihood". See ?mle for examples, in particular there is an example that starts "## Linear regression using MLE" which provides a way to write a linear regression, which I am assuming you want.

Therefore the log-likelihood can be written (with the error log-transformed to keep it positive) as:

nll <- function(b0, b1, logsd) {    
  mu <- cbind(1, log(DF$K)) %*% c(b0, b1) ;   
  -sum(dnorm(log(DF$Y), mu, exp(logsd), log=TRUE)) 
  }

and estimated with

stats4::mle(minuslog=nll, start=c(0,0,1))

The bblme package offers a formula notation which you may find easier to use.

library(bbmle)
mod2 <- mle2(Y ~ dnorm(mean=X %*% c(b0, b1), sd=exp(logsd)), 
           start=list(b0=0,b1=0,logsd=1), # use named list of parameters
           data=list(X=cbind(1, log(DF$K)), Y=log(DF$Y)))
summary(mod2)

It can also be used with a similar syntax to stats4::mle but allows you to pass vectors of parameters and data arguments, which can make the code a bit cleaner.

nll2 <- function(par) {
   mu <- X %*% par[1:2] ;
  -sum(dnorm(Y, mu, exp(par[3]), log=TRUE))
}

# set the parameter names & set `vecpar` to TRUE
parnames(nll2) <- c("b0", "b1", "logsd")
mod3 <- mle2(nll2,
         start=list(b0=0,b1=0,logsd=1),
         data=list(X=cbind(1, log(DF$K)), Y=log(DF$Y)), vecpar=TRUE)
summary(mod3)