[R] optim() for ordered logit model with parallel regression assumption

[R. Michael Weylandt]

On Wed, Aug 1, 2012 at 4:34 PM, Xu Jun <junxu.r at gmail.com> wrote:
> Thanks Michael. Now I switched my approach after doing some google.
> Following are my new codes:
>
> ###################################################
>  library(foreign)
>  readin <- read.dta("ordfile.dta", convert.factors=FALSE)
>  myvars <- c("depvar", "x1", "x2", "x3")
>  mydta <- readin[myvars]
>  # remove all missings
> mydta <- na.omit(mydta)
>
> ologit.lf <- function(theta, y, X) {
>   X <- as.matrix(X)
>   y <- as.matrix(y)
>   n <- nrow(X)
>   k <- ncol(X)
>   b <- theta[1:k]
>   tau1 <- as.numeric(theta [k+1])
>   tau2 <- as.numeric(theta [k+2])
>   tau3 <- as.numeric(theta [k+3])
>
>   p1 <- (1/(1+exp( - tau1 + X %*% b)))
>   p2 <- (1/(1+exp( - tau2 + X %*% b))) - (1/(1+exp( - tau1 + X %*% b)))
>   p3 <- (1/(1+exp( - tau3 + X %*% b))) - (1/(1+exp( - tau2 + X %*% b)))
>   p4 <- 1 - (1/(1+exp( - tau3 + X %*% b)))
>
>   # create indicator variables
>   i  <- rep(1, nrow(X))
>   i1 <- as.numeric(  i == y)
>   i2 <- as.numeric(2*i == y)
>   i3 <- as.numeric(3*i == y)
>   i4 <- as.numeric(4*i == y)
>
>   llik <- sum(i1*log(p1) + i2*log(p2) + i3*log(p3) + i4*log(p4))
>   return(-llik)
> }
>
> y <- (mydta$depvar)
> X <- cbind(mydta$x1, mydta$x2, mydta$x3)
>
> if (FALSE) { # providing start values in optim
>   xint <-  cbind(rep(1,nrow(X)),X)
>   bols <- (solve(t(xint) %*% xint)) %*% (t(xint) %*% y)
>   startval <- rbind(as.matrix(bols[2:nrow(bols)]),0,0,0)
> }
>
> runopt <- optim(rep(0, ncol(X)+3), ologit.lf, method="BFGS",
> hessian=T, y=y, X=X)
> #########################################################################
>
> But then I got the following error message;
> Error in optim(rep(0, ncol(X) + 3), ologit.lf, method = "BFGS", hessian = T,  :
>   initial value in 'vmmin' is not finite
>
> I know there is something wrong with the way that I set up my initial
> values, but just couldn't figure out how. Any help would be greatly
> appreciated!
>

Without having data I can't reproduce [1] but I think your problem is
likely that p2 or p3 == 0 so the log is non-finite.

Also:  I'd imagine this has already been implemented in some package,
but I'm not familiar enough with the domain to give you any pointers.

Best,
Michael

[1]  http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example