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]]
______________________________________________
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.
