On Wed, Aug 1, 2012 at 4:34 PM, Xu Jun <junx...@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 ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.