R help - Nov 2012 - Invalid 'times' argument three-category ordered probit with maximum likelihood

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=",")
attach(econrate)

#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.

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