Profile Confidence Intervals in R: mle2

Question

I am trying to use the mle2 command in the bbmle package. I look at p2 "Maximum Likelihood Assessment and Analysis with the bbmle " from Bolker. For some reason I can’t enter the correct starting values. Here's the reproducible code:

 l.lik.probit <-function(par, ivs, dv){ Y <- as.matrix(dv) X <- as.matrix(ivs) K <-ncol(X) b <- as.matrix(par[1:K]) phi <- pnorm(X %*% b) sum(Y * log(phi) + (1 - Y) * log(1 - phi)) } n=200 set.seed(1000) x1 <- rnorm(n) x2 <- rnorm(n) x3 <- rnorm(n) x4 <- rnorm(n) latentz<- 1 + 2.0 * x1 + 3.0 * x2 + 5.0 * x3 + 8.0 * x4 + rnorm(n,0,5) y <- latentz y[latentz < 1] <- 0 y[latentz >=1] <- 1 x <- cbind(1,x1,x2,x3,x4) values.start <-c(1,1,1,1,1) foo2<-mle2(l.lik.probit, start=list(dv=0,ivs=values.start),method="BFGS",optimizer="optim", data=list(Y=y,X=x)) 

And this is the error I get:

 Error in mle2(l.lik.probit, start = list(Y = 0, X = values.start), method = "BFGS", : some named arguments in 'start' are not arguments to the specified log-likelihood function 

Any idea why? Thank you for your help!

Answer 1

You missed a couple of things, but most importantly, by default, mle2 accepts a list of options; you can force it to take the parameter vector instead, but you will have to work a little harder.

I changed the code a bit in places. (I changed the log-likelihood function to a negative log-likelihood function, without which it would never work!)

 l.lik.probit <-function(par, ivs, dv){ K <- ncol(ivs) b <- as.matrix(par[1:K]) phi <- pnorm(ivs %*% b) -sum(dv * log(phi) + (1 - dv) * log(1 - phi)) } n <- 200 set.seed(1000) dat <- data.frame(x1=rnorm(n), x2=rnorm(n), x3=rnorm(n), x4=rnorm(n)) beta <- c(1,2,3,5,8) mm <- model.matrix(~x1+x2+x3+x4,data=dat) latentz<- rnorm(n,mean=mm%*%beta,sd=5) y <- latentz y[latentz < 1] <- 0 y[latentz >=1] <- 1 x <- mm values.start <- rep(1,5) 

Now we come. The main vecpar=TRUE is to specify vecpar=TRUE and use parnames so that mle2 knows the names of the elements in the parameter vector ...

 library("bbmle") names(values.start) <- parnames(l.lik.probit) <- paste0("b",0:4) m1 <- mle2(l.lik.probit, start=values.start, vecpar=TRUE, method="BFGS",optimizer="optim", data=list(dv=y,ivs=x)) 

As indicated above for this specific example, you just re-performed probit restriction (although I understand that now you want to expand this to somehow provide heterosedasticity ...)

 dat2 <- data.frame(dat,y) m2 <- glm(y~x1+x2+x3+x4,family=binomial(link="probit"), data=dat2) 

As a final note, I would say that you should check the arguments argument, which allows you to specify a sublinear model for any of the parameters and the formula interface:

 m3 <- mle2(y~dbinom(prob=pnorm(eta),size=1), parameters=list(eta~x1+x2+x3+x4), start=list(eta=0), data=dat2) 

PS confint(foo2) seems to work fine (setting CI profiles on demand) with this setting.

 ae <- function(x,y) all.equal(unname(coef(x)),unname(coef(y)),tol=5e-5) ae(m1,m2) && ae(m2,m3)