Creating function arguments from a named list (with an appendix to stats4 :: mle)

Question

I have to start with what I'm trying to do: I want to use the mle function without rewriting my likelihood function every time I want to try a different model specification. Since mle expects a named list of initial values, you apparently cannot just write a log-likelihood function as a transfer of the parameter vector. A simple example:

Suppose I want to fit a linear regression model using maximum likelihood, and first I ignore one of my predictors:

n <- 100 df <- data.frame(x1 = runif(n), x2 = runif(n), y = runif(n)) Y <- df$y X <- model.matrix(lm(y ~ x1, data = df)) # define log-likelihood function ll <- function(beta0, beta1, sigma){ beta = matrix(NA, nrow=2, ncol=1) beta[,1] = c(beta0, beta1) -sum(log(dnorm(Y - X %*% beta, 0, sigma))) } library(stats4) mle(ll, start = list(beta0=.1, beta1=.2, sigma=1) 

Now, if I want to customize another model, say:

 m <- lm(y ~ x1 + x2, data = df) 

I cannot reuse the log-likelihood function - I had to rewrite it to have beta3 parameter. What I would like to do is something like:

 ll.flex <- function(theta){ # theta is a vector that I can use directly ... } 

if I could somehow adjust the initial argument in mle to account for my now logarithmic likelihood function with a vector input or banning what has a function that builds a logarithmic likelihood function at runtime, say by building a named argument list, and then using it to define a function, for example, something like this:

 X <- model.matrix(lm(y ~ x1 + x2, data = df)) arguments <- rep(NA, dim(X)[2]) names(arguments) <- colnames(X) ll.magic <- function(bring.this.to.life.as.function.arguments(arguments)){...} 

Update:

I ended up writing a helper function that can add an arbitrary number of named arguments x1, x2, x3 ... to the passed function f.

 add.arguments <- function(f,n){ # adds n arguments to a function f; returns that new function t = paste("arg <- alist(", paste(sapply(1:n, function(i) paste("x",i, "=",sep="")), collapse=","), ")", sep="") formals(f) <- eval(parse(text=t)) f } 

This is ugly, but he did his job, allowing me to reuse the log likelihood function on the fly.

Answer 1

You can use the mle2 function from the mle2 package, which allows you to pass vectors as parameters. Here is a sample code.

 # REDEFINE LOG LIKELIHOOD ll2 = function(params){ beta = matrix(NA, nrow = length(params) - 1, ncol = 1) beta[,1] = params[-length(params)] sigma = params[[length(params)]] minusll = -sum(log(dnorm(Y - X %*% beta, 0, sigma))) return(minusll) } # REGRESS Y ON X1 X <- model.matrix(lm(y ~ x1, data = df)) mle2(ll2, start = c(beta0 = 0.1, beta1 = 0.2, sigma = 1), vecpar = TRUE, parnames = c('beta0', 'beta1', 'sigma')) # REGRESS Y ON X1 + X2 X <- model.matrix(lm(y ~ x1 + x2, data = df)) mle2(ll2, start = c(beta0 = 0.1, beta1 = 0.2, beta2 = 0.1, sigma = 1), vecpar = TRUE, parnames = c('beta0', 'beta1', 'beta2', 'sigma')) 

It gives you

 Call: mle2(minuslogl = ll2, start = c(beta0 = 0.1, beta1 = 0.2, beta2 = 0.1, sigma = 1), vecpar = TRUE, parnames = c("beta0", "beta1", "beta2", "sigma")) Coefficients: beta0 beta1 beta2 sigma 0.5526946 -0.2374106 0.1277266 0.2861055 

Answer 2

It may be easier to use optim directly; which is mle anyway.

 ll2 <- function(par, X, Y){ beta <- matrix(c(par[-1]), ncol=1) -sum(log(dnorm(Y - X %*% beta, 0, par[1]))) } getp <- function(X, sigma=1, beta=0.1) { p <- c(sigma, rep(beta, ncol(X))) names(p) <- c("sigma", paste("beta", 0:(ncol(X)-1), sep="")) p } set.seed(5) n <- 100 df <- data.frame(x1 = runif(n), x2 = runif(n), y = runif(n)) Y <- df$y X1 <- model.matrix(y ~ x1, data = df) X2 <- model.matrix(y ~ x1 + x2, data = df) optim(getp(X1), ll2, X=X1, Y=Y)$par optim(getp(X2), ll2, X=X2, Y=Y)$par 

With exit

 > optim(getp(X1), ll2, X=X1, Y=Y)$par sigma beta0 beta1 0.30506139 0.47607747 -0.04478441 > optim(getp(X2), ll2, X=X2, Y=Y)$par sigma beta0 beta1 beta2 0.30114079 0.39452726 -0.06418481 0.17950760 

Answer 3

This may not be what you are looking for, but I would do it like this:

 mle2(y ~ dnorm(mu, sigma),parameters=list(mu~x1 + x2), data = df, start = list(mu = 1,sigma = 1)) mle2(y ~ dnorm(mu,sigma), parameters = list(mu ~ x1), data = df, start = list(mu=1,sigma=1)) 

You may be able to adapt this formulation for polynomials, although dmultinom may not work - you may need to write Dmultinom() , which takes a matrix from multidimensional samples and returns probability (log).

Answer 4

The R code provided by Ramnath can also be applied to the optimization function, since it also accepts vectors.