[R] Invalid 'times' argument three-category ordered probit with maximum likelihood

[David Winsemius]

On Nov 12, 2012, at 2:16 PM, jackle1o4 wrote:

> Hello,
> First time poster here so let me know if you need any more information. I am
> trying to run an ordered probit with maximum likelihood model in R with a
> very simple model (model <- econ3 ~ partyid). Everything looks ok until i
> try to run the optim() command and that's when I get " Error in rep(1,
> nrow(x)) : invalid 'times' argument". I had to adapt the code from a 4
> category likelihood and I have a suspicion that the problem is in there. The
> data set consists of two variables with 558 observations. Econ3 is a 1-3
> rating and partyid is a range from -3(Strong Democrat) to 3 (Strong
> Republican). Here is the code I used:
> 
> setwd("C:/Users/Terry/Desktop/Terry/School/Fall 2012/ML/HW")
> 
> #Load Libraries
> library(MASS)
> library(tile)
> library(simcf)
> 
> #Load Data
> econrate <- read.csv("hw4econ3.csv", header=TRUE, sep=",")

It's somewhat unusual to see a posting with an effort at posting code go uncommented for 18 hours, but in your case I suspect it is because people got to this point an were noticing the sprong resemblance to a homework assignment, which is discouraged on Rhelp. Speaking of a different sort of likelihood ... Help with assignments are not completely refused, but the chances of getting a reply would increase if you had indicated what academic institution you were studying at and what their policies are for soliciting help with academic assignments. Further increases in the likelihood of a response would occur if the h4econ3.csv file were made available at an URL or attached as a text file or included in the body of the message using the dput function.


> attach(econrate)

Generally a bad idea to attach() data objects, and that is especially so when you have not described the structure of the dataframe.
> 
> #Ordered Probit Liklihood
> llk.oprobit3 <- function(param, x, y) {
> os <- rep(1, nrow(x))
>  x <- cbind(os, x)  
>  b <- param[1:ncol(x)]
>  t2 <- param[(ncol(x)+1)]
> xb <- x%*%b
>  p1 <- log(pnorm(-xb))
>  if (t2<=0)  p2 <- -(abs(t2)*10000)
> else p2 <- log(pnorm(t2-xb)-pnorm(-xb))  
>  p3 <- log(1-pnorm(t2-xb))
> -sum(cbind(y==1,y==2,y==3) * cbind(p1,p2,p3))
> }
> 
> #Define Data
> y <- econ3
> x <- partyid
> model <- (econ3 ~ partyid)
> 
> #Use optim directly
> ls.result <- lm(y~x)
> stval <- c(ls.result$coefficients, 1)
> oprobit.result <- optim(stval, llk.oprobit3, method="BFGS", x=x, y=y,
> hessian=T) ###Here is where it all breaks down
> pe <- oprobit.result$par
> vc <- solve(oprobit.result$hessian)
> se <- sqrt(diag(vc)) 
> ll <- -oprobit.result$value  
> 
> Any help would be greatly appreciated. Thanks.
> 
> 
> 
> --
> View this message in context: http://r.789695.n4.nabble.com/Invalid-times-argument-three-category-ordered-probit-with-maximum-likelihood-tp4649341.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> ______________________________________________
> R-help at 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.

David Winsemius, MD
Alameda, CA, USA