This is R-help mailing list, not the do-my-work-for-me mailing list. You need
to specify what you expect it to do and where it is going wrong, because wading
through a lot of code on a wild goose chase is a waste of anyone else's
time.
---------------------------------------------------------------------------
Jeff Newmiller The ..... ..... Go Live...
DCN:<jdnewmil at dcn.davis.ca.us> Basics: ##.#. ##.#. Live
Go...
Live: OO#.. Dead: OO#.. Playing
Research Engineer (Solar/Batteries O.O#. #.O#. with
/Software/Embedded Controllers) .OO#. .OO#. rocks...1k
---------------------------------------------------------------------------
Sent from my phone. Please excuse my brevity.
Charles <fffchao at hotmail.com> wrote:
>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 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.