setwd("D:/WZhao/WZhaoStJude/Paper Submission/Multi-reviewers/Code to publish") source("EM.R") source("likelihood.R") source("likelihood_o.R") ######### simulation study ############# ## define the true acumen credibility matrix p1=c(0.5, 0.5, 0.5, 0.5) p2=c(0.5, 0.5, 0.5, 0.5) p3=c(0.5, 0.5, 0.5, 0.5) p4=c(0.9, 0.9, 0.9, 0.9) p5=c(0.9, 0.9, 0.9, 0.9) p6=c(0.9, 0.9, 0.9, 0.9) pmatrix=cbind(p1, p2, p3, p4, p5, p6) m=6 ## number of reviewers K=4 ## number of categories no.sample=100 ## number of samples ratio=c(.7, .4, 0.2) ## 30% in cat 4, 30% in cat 3, 20% in cat 2, and 20% in cat 1. options(K=K) no.iter=1:1000 data.work=matrix(0, ncol=m, nrow=no.sample) no.iter=1000 for (jj in no.iter){ print(paste("simu", jj)) simu=jj ## simulate data for (i in 1:no.sample) { for (j in 1:m) { ctg=4; p=pmatrix[ctg,j]; p=c(rep((1-p)/(K-1), K-1),p) if(i<=no.sample*ratio[1]){ ctg=3; p=pmatrix[ctg,j]; p=c((1-p)/(K-1),(1-p)/(K-1),p,(1-p)/(K-1))} if(i<=no.sample*ratio[2]){ ctg=2; p=pmatrix[ctg,j]; p=c((1-p)/(K-1),p,(1-p)/(K-1), (1-p)/(K-1))} if(i<=no.sample*ratio[3]){ ctg=1; p=pmatrix[ctg,j]; p=c(p, rep((1-p)/(K-1), K-1))} data.work[i,j]=sample(1:K, 1, prob=p) } options(data.work=data.work) ## assign initial values ## all credibility is 0.5 ## true grade is uniform=1/K m=dim(data.work)[[2]] temp=matrix(0.1,ncol=K, nrow=K) diag(temp)=0.5 temp=temp[, -K] credmatrix=rep(as.vector(temp),m) options(credmatrix=credmatrix) pi1=matrix(1/K, ncol=K, nrow=dim(data.work)[[1]]) options(pi1=pi1) EM() credmatrix1=options()$credmatrix ## collect the credebility matrix for each simulation result=rbind(result, credmatrix1) }