[R] Discrete choice model maximum likelihood estimation

[Berend Hasselman]

See below.
On 14-05-2012, at 13:21, infinitehorizon wrote:

> Hello again,
> 
> I changed the name to tt. 
> and for a and tt actually I was getting them from data, I didn't put them
> here in the question. Now I restructured my code and below I copy the full
> code, I tried many things but still getting the same error, I don't
> understand where is the mistake.
> 
> I also defined one more variable to increase comprehension of the problem.
> Instead of data, I define three representative vectors in the beginning. If
> you run this code, you will see the error message.
> 
> # Variables, a: discrete choice variable-dependent, x and tt are independent
> variables, tt is binary
> a  <- c(1,1,2,3,2,3,1,2,2,3,1,1)
> x  <- c(23,26,12,27,10,30,13,20,23,44,17,15)
> tt <- c(1,0,0,1,1,1,0,0,1,0,1,1)
> 
> # First, just to see, the linear model
> 
> LM	<-lm(a ~ x+tt)
> coefLM	<- coefficients(LM)
> summary(LM)
> 
> # Probabilities for discrete choices for a=3, a=2 and a=1 respectively 
> P3 <- function(bx,b3,b,tt) { 
> P <- exp(bx*x+b3+b*(tt==1))/(1-exp(bx*x+b3+b*(tt==1))) 
> return(P) 
> } 
> P2 <- function(bx,b2,b,tt) { 
> P <- exp(bx*x+b2+b*(tt==1))/(1-exp(bx*x+b2+b*(tt==1))) 
> return(P) 
> } 
> P1 <- function(bx,b1,b,tt) { 
> P <- exp(bx*x+b1+b*(tt==1))/(1-exp(bx*x+b1+b*(tt==1))) 
> return(P) 
> }
> 
> # Likelihood functions for discrete choices for a=3, a=2 and a=1
> respectively
> 
> L3 <- function(bx,b1,b2,b3,b,tt) { 
> P11 <- P1(bx,b1,b,tt) 
> P22 <- P2(bx,b2,b,tt) 
> P33 <- P3(bx,b3,b,tt) 
> 
> L3l <- P11*P22*P33 
> return(L3l) 
> } 
> 
> L2 <- function(bx,b1,b2,b,tt) { 
> P11 <- P1(bx,b1,b,tt) 
> P22 <- P2(bx,b2,b,tt) 
> 
> L2l <- P11*P22 
> return(L2l) 
> } 
> 
> L1 <- function(bx,b1,b,tt) { 
> P11 <- P1(bx,b1,b,tt) 
> 
> L1l <- P11 
> return(L1l) 
> }
> 
> # Log-likelihood function 
> 
> llfn <- function(param) { 
> 
> bx <- param[1] 
> b1 <- param[2] 
> b2 <- param[3] 
> b3 <- param[4] 
> b <- param[5] 
> 
> lL1 <- log(L1(bx,b1,b2,b,tt)) 
> lL2 <- log(L2(bx,b1,b2,b3,b,tt)) 
> lL3 <- log(L3(bx,b1,b2,b3,b,tt)) 
> 
> llfn <- (a==1)*lL1+(a==2)*lL2+(a==3)*lL3 
> } 
> start.par <- c(0,0,0,0,0) 
> est <- optim(param=start.par,llfn,x=x, a=a, tt=tt, method =
> c("CG"),control=list(trace=2,maxit=2000), hessian=TRUE)
> 

Due to urgent matters, I can only answer briefly.

optim doesn't have an argument param. It does have an argument par so you should have written par=start.par.

If you do that you will get other error messages. You are calling L1, L2 and L3 with too many arguments.

And then you will find that the return value of llfn is a vector and not a scalar.

Finally why are passing x=x in the  optim call to llfn? It is not used anywhere.

You need to rethink your approach. And certainly read the help for optim.

Berend