Five Category Ordered Probit Log Likelihood Function in R Help : Survey Data Research

#1
I am currently working with some survey data that has a 5 category scale. I have used a ML method before for a 4 category scale and I am trying to adapt the code for the new survey in R.

The original code that worked was:


#Load Data
mathcomp <- read.csv("RRawComposite11.csv", header=TRUE, sep=",")
attach(mathcomp)

#Ordered Probit Liklihood
llk.oprobit4 <- function(param, x, y) {
os <- rep(1, nrow(x))
x <- cbind(os, x)
b <- param[1:ncol(x)]
t2 <- param[(ncol(x)+1)]
t3 <- param[(ncol(x)+2)]
xb <- x%*%b
p1 <- log(pnorm(-xb))
if (t2<=0) p2 <- -(abs(t2)*10000)
else p2 <- log(pnorm(t2-xb)-pnorm(-xb))
if (t3<=t2) p3 <- -((t2-t3)*10000)
else p3 <- log(pnorm(t3-xb)-pnorm(t2-xb))
p4 <- log(1-pnorm(t3-xb))
-sum(cbind(y==1,y==2,y==3,y==4) * cbind(p1,p2,p3,p4))
}

#Define Data
y <- MORE
x <- cbind(GENDER, FREQCAL, FREQPC, INT, EXT, POSS, USE)
model <- (MORE ~ GENDER + FREQCAL + FREQPC + INT + EXT + POSS + USE)

#Use optim directly
ls.result <- lm(y~x)
stval <- c(ls.result$coefficients,1,2)
oprobit.result <- optim(stval, llk.oprobit4, method="BFGS", x=x, y=y, hessian=T)
pe <- oprobit.result$par
vc <- solve(oprobit.result$hessian)
se <- sqrt(diag(vc))
ll <- -oprobit.result$value

The modified code I am trying to use is:


llk.oprobit5 <- function(param, x, y) {
os <- rep(1, nrow(x))
x <- cbind(os, x)
b <- param[1:ncol(x)]
t2 <- param[(ncol(x)+1)]
t3 <- param[(ncol(x)+2)]
t4 <- param[(ncol(x)+3)]
xb <- x%*%b
p1 <- log(pnorm(-xb))
if (t2<=0) p2 <- -(abs(t2)*10000)
else p2 <- log(pnorm(t2-xb)-pnorm(-xb))
if (t3<=t2) p3 <- -((t2-t3)*10000)
else p3 <- log(pnorm(t3-xb)-pnorm(t2-xb))
if (t4<=t3) p4 <- -((t3-t4)*10000)
else p4 <- log(pnorm(t4-xb)-pnorm(t3-xb))
p5 <- log(1-pnorm(t4-xb))
-sum(cbind(y==1,y==2,y==3,y==4,y==5) * cbind(p1,p2,p3,p4,p5))
}

y <- q9
x <- cbind(Level,PIAP,Postgraduate,Sex,Age,Region)
model <- (q9 ~ Level + PIAP + Postgraduate + Sex + Age + Region)

#Use optim directly
ls.result <- lm(y~x)
stval <- c(ls.result$coefficients,1,2)
oprobit.result <- optim(stval, llk.oprobit5, method="BFGS", x=x, y=y, hessian=T)
pe <- oprobit.result$par
vc <- solve(oprobit.result$hessian)
se <- sqrt(diag(vc))
ll <- -oprobit.result$value

I am getting an error message saying:

Error in if (t4 <= t3) p4 <- -((t3 - t4) * 10000) else p4 <- log(pnorm(t4 - :
missing value where TRUE/FALSE needed

Does anyone know how to fix the likelihood function to make this work? Any help would be appreciated. Thanks.