code problem with the optim() function

Please forgive me if someone has seen this duplicated email, because I
sent
the email to the wrong address just now.

Hi, all,

I'm estimating an inter-rater coefficient, Aickin (1990)'s alpha. It's
been
a long time for me to figure out how to make it work, but it's still of
no
avail. The problem seems to be my code with the optim() function. Can
anyone help me take a look at the code?

Thanks in advance.
Best,
Charles

x3<-c(rep(1:2,10),rep(1,5),rep(2,5))
x4<-c(rep(1:2,10),rep(2,5),rep(1,5))
ratings<-as.data.frame(cbind(x3,x4))
ratings <- as.matrix(na.omit(ratings))
ns <- nrow(ratings)
nr <- ncol(ratings)
r1 <- ratings[, 1]
r2 <- ratings[, 2]

if (!is.factor(r1))
r1 <- factor(r1)
if (!is.factor(r2))
r2 <- factor(r2)
ifelse (length(levels(r1)) >= length(levels(r2)),lev <- c(levels(r1),
levels(r2)), lev <- c(levels(r2), levels(r1)))
lev <- lev[!duplicated(lev)]

r1 <- factor(ratings[, 1], levels = lev)
r2 <- factor(ratings[, 2], levels = lev)
ttab <- table(r1, r2)
total<-margin.table(ttab)
pa<-(margin.table(ttab,1)/total)
pb<-(margin.table(ttab,2)/total)

ncate<-length(levels(as.factor(as.matrix(ratings))))
pp<-prop.table(ttab)
sumagr<-sum(diag(ttab))/ns
alpha<-(sumagr-sum(pa*pb))/(1-sum(pa*pb))


aickin <- function(theta,ratings) {
alphanew<-theta[1]
pah<-theta[2:3]
pbh<-theta[4:5]
ratings <- as.matrix(na.omit(ratings))
ns <- nrow(ratings)
nr <- ncol(ratings)
r1 <- ratings[, 1]
r2 <- ratings[, 2]
if (!is.factor(r1)) r1 <- factor(r1)
if (!is.factor(r2)) r2 <- factor(r2)
ifelse (length(levels(r1)) >= length(levels(r2)),lev <- c(levels(r1),
levels(r2)), lev <- c(levels(r2), levels(r1)))

lev <- lev[!duplicated(lev)]
r1 <- factor(ratings[, 1], levels = lev)
r2 <- factor(ratings[, 2], levels = lev)
ttab <- table(r1, r2)
total<-margin.table(ttab)
#if(total==0) return(0);
pa<-(margin.table(ttab,1)/total)
pb<-(margin.table(ttab,2)/total)

#pa[1]*pb[1]+pa[2]*pb[2]+pa[3]*pb[3]+pa[4]*pb[4]

pp<-prop.table(ttab)
sumagr<-sum(diag(ttab))/ns
alpha<-(sumagr-sum(pa*pb))/(1-sum(pa*pb))
for (i in 1:nrow(ttab)){
for (j in 1:ncol(ttab)){
d<-array(NA,dim=c(nrow(pp),ncol(pp)))
d[i,j]<-ifelse(i==j,1,0)
po<-array(NA,dim=c(nrow(pp),ncol(pp)))
 po[i,j]<-d[i,j]*pp[i,j]
for (i in 1:nrow(ttab)){
for (j in 1:ncol(ttab)){
d[i,j]<-ifelse(i==j,1,0)
po[i,j]<-d[i,j]*pp[i,j]
}
}
     }
}
      pseudo<-1
     total<-total+1
     ttabnew<-ttab+pseudo/(ncol(ttab)*nrow(ttab))
     ttabnew<-prop.table(ttabnew)

 pah<-(margin.table(ttabnew,1))
pbh<-(margin.table(ttabnew,2))


s<- foreach (ii =1:4) %dopar% {
for (i in 1:ncol(ttabnew)){
for (j in 1:ncol(ttabnew)){
s<-rep(NA,2)
s[ii]<-d[i,j]*pah[i]*pbh[j]
}
ssum<-sum(as.vector(unlist(s)))
pbhh<-NA
pahh<-NA
pbhh[j]<- d[i,j]*pbh[j]
pahh[i]<-d[i,j]*pah[i]
sumpbh<-sum(pbhh)
     sumpah<-sum(pahh)
pah[i]<-pa[i]/(1-alphanew+(alphanew*sumpbh/ssum) )
pbh[j]<-pb[j]/(1-alphanew+(alphanew*sumpah/ssum) )
#  s<-sum(d[i,j]*pah[i]*pbh[j])
pp[i,j]<-(1-alphanew+(alphanew*d[i,j]/ssum ))*pah[i]*pbh[j]
logpp<-log(pp[i,j])
}
}
 return(-logpp)
}

optim(c(alpha,as.vector(pa),as.vector(pb)),aickin,ratings=ratings,method="BFGS")

[[alternative HTML version deleted]]