####DESCRIPTION #The code is composed by 3 steps (R-codes) #step 1: creating "rhobetamodel.bug" (put the file in the working directory) #step 2: (auxiliary) posterior probability "posteriorprob" function #step 3: running Bugs routine "rhohierarchicalmodel" # require(R2WinBUGS); # #step 1: building and saving the model--> in archive "rhobetamodel.bug" rhobetamodel<-function() { for(i in 1:I){ for(j in 1:J){ Y[i,j]~dnorm(a[i],tau.e); } a[i]~dnorm(mu,tau.a) } #priors(in exponential component we put parameter=1. Second parameter) #error precision tau.e~dgamma(alpha,1); #variance var.e<-pow(tau.e,-1); # #class precision tau.a~dgamma(beta,1); #variance var.a<-pow(tau.a,-1); # mu~dnorm(muo,tauo); #rho has the beta distribution with parameters: alpha and beta; #put alpha=first parameter of tau.e #put beta=first parameter of tau.a rho<-var.a/(var.a+var.e); } if (is.R()){ # for R ## some temporary filename: filename <- file.path(tempdir(), "rhobetamodel.bug") } else{ # for S-PLUS ## put the file in the working directory: filename <- "rhobetamodel.bug" } #} write.model(rhobetamodel, "rhobetamodel.bug") ## show the model... file.show("rhobetamodel.bug"); #end step 1 # #step 2: posterior probability: rho<=cust posteriorprob<-function(sample,cust){ outputp<-matrix(c(0),nrow=length(cust),ncol=2,dimnames=list(c(),c("cust","P(rho<=cust)"))); outputp[,1]<-cust; for(i in 1: length(cust)){ auxiliar<-matrix(c(0),nrow=length(sample),ncol=1); for(k in 1:length(sample)){ if(sample[k]<=cust[i]){auxiliar[k]<-1} } outputp[i,2]<-mean(auxiliar); } show(outputp) } #end step 2 # #### #step 3: to run Bugs rhohierarchicalmodel<-function(subs,cust,file,I,J,alpha,beta,muo,tauo) { if(missing(alpha)){alpha<-1}; if(missing(beta)){beta<-1}; # alpha=beta=1-->rho uniform(0,1) if(missing(muo)){muo<-0.0}; if(missing(tauo)){tauo<-1e-06};# muo=0, tauo=1e-06-->noninformative prior # #bugs parameters # J<-J;# replications (the same!) by class; I<-I; data1<-scan(file); x<-matrix(c(data1),nrow=J,ncol=I); Y<-t(x); #bugs function needs: I x J data matrix; if(subs=="null"){ I<-I # class (or sample) 1,..., I; Y<-Y; newsam<-c(1:I); # }else{newsam<-sample(c(1:I),subs); Y<-Y[newsam,]; I<-subs # class (or sample) 1,..., subs ; }; # # data2 <- list ("I", "J", "Y","alpha","beta","muo","tauo"); # #random initial values. Recom: Bayesian Data Analysis. Gelman et al (2004, 2 ed.) initsf<-function(){list(a=rnorm(I,0,100),tau.e=runif(1,0,100),tau.a=runif(1,0,100),mu=rnorm(1,0,100))}; inits<-list(inits1=initsf(),inits2=initsf(),inits3=initsf()); parameters <- c("a", "tau.e", "tau.a", "mu","var.e", "var.a", "rho"); rhobetamodel.sim <- bugs (data2, inits, parameters, "rhobetamodel.bug", n.chains=3, n.iter=500000,n.thin=50); windows(); plot(rhobetamodel.sim); ## # ex.cs1 <- expression(plain(pi(rho)), plain(pi(rho / data)))# ## windows(); Max1<-max(density(rhobetamodel.sim\$sims.list\$rho)\$y); if(alpha>1 & beta>1){modabeta<-(alpha-1)/(alpha+beta-2)}; if(alpha==1 & beta==1){modabeta<-0}; if(alpha<1 & beta>=1){modabeta<-0}; if(alpha>=1 & beta<1){modabeta<-1}; Max2<-dbeta(modabeta,alpha,beta); Max12<-max(Max1,Max2)+1; if(Max12=="Inf"){Max12<-Max1+1}; hist(rhobetamodel.sim\$sims.list\$rho,freq=FALSE,main=paste("Histogram (posterior sample) and densities [0,1]"),xlim=c(0,1),ylim=c(0,Max12),xlab=expression(rho),ylab=expression(density(rho))); lines(density(rhobetamodel.sim\$sims.list\$rho),col="red",lty=2); curve(dbeta(x,alpha,beta),0,1,add=TRUE,col="blue"); legend("topleft",ex.cs1,lty=1:2,col=c("blue","red")); print(rhobetamodel.sim); ## output<-list(data=Y,newsam=newsam,simulation=rhobetamodel.sim,postrhosample=rhobetamodel.sim\$sims.list\$rho, PProbrho=posteriorprob(rhobetamodel.sim\$sims.list\$rho,cust)); } #end step 3 # #usage: rhohierarchicalmodel(subs,cust,file,I,J,alpha,beta,muo,tauo) #arguments: #subs: #file: dataset (I x J data matrix) #I: samples #J: replications #alpha: parameter of beta distribution (rho) #beta: parameter of beta distribution (rho) #muo, tauo: location and precision parameters (a~dnorm(muo,tauo)) #### ####examples: #MMM<-rhohierarchicalmodel(subs=8,c(0.7,0.8,0.9),"datalimon",60,3,10,3,5,1); #MMM<-rhohierarchicalmodel(subs="null",c(0.7),"datalimon",60,3,2,3); #MMM<-rhohierarchicalmodel(subs=8,c(0.7,0.88),"datalimon",60,3); #MMM1<-rhohierarchicalmodel(subs=5,c(0.72,0.74,0.77,0.8,0.83,0.85,0.88,0.9),"datalimon",60,3); #MMM2<-rhohierarchicalmodel(subs=5,c(0.72,0.74,0.77,0.8,0.83,0.85,0.88,0.9),"datalimon",60,3,2,2); #MMM3<-rhohierarchicalmodel(subs=8,c(0.72,0.74,0.77,0.8,0.83,0.85,0.88,0.9),"datalimon",60,3,1,10); #MMM4<-rhohierarchicalmodel(subs=8,c(0.72,0.74,0.77,0.8,0.83,0.85,0.88,0.9),"datalimon",60,3,10,1); #MMM5<-rhohierarchicalmodel(subs=8,c(0.72,0.74,0.77,0.8,0.83,0.85,0.88,0.9),"datalimon",60,3,1,1,5,1); #MMM6<-rhohierarchicalmodel(subs=8,c(0.72,0.74,0.77,0.8,0.83,0.85,0.88,0.9),"datalimon",60,3,2,10); #MMM7<-rhohierarchicalmodel(subs=8,c(0.72,0.74,0.77,0.8,0.83,0.85,0.88,0.9),"datalimon",60,3,10,2); #rodando muestra paper #M10<-rhohierarchicalmodel(subs=10,c(0.72,0.74,0.77,0.8,0.83,0.85,0.88,0.9),"datalimon",60,3); # MM1022T07<-rhohierarchicalmodel(subs="null",c(0.05815,0.09535,0.25615,0.5086,0.7696,0.91135),"sample10",10,3,2,2); Read 30 items Inference for Bugs model at "rhobetamodel.bug", fit using WinBUGS, 3 chains, each with 5e+05 iterations (first 250000 discarded), n.thin = 50 n.sims = 15000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff a[1] 5.0 0.2 4.7 4.9 5.0 5.1 5.3 1 15000 a[2] 5.0 0.2 4.6 4.9 5.0 5.1 5.3 1 15000 a[3] 5.2 0.2 4.8 5.0 5.2 5.3 5.5 1 7500 a[4] 5.1 0.2 4.8 5.0 5.1 5.2 5.4 1 8700 a[5] 5.3 0.2 4.9 5.1 5.3 5.4 5.6 1 15000 a[6] 5.6 0.2 5.3 5.5 5.6 5.7 5.9 1 11000 a[7] 5.3 0.2 5.0 5.2 5.3 5.4 5.7 1 15000 a[8] 5.4 0.2 5.0 5.3 5.4 5.5 5.7 1 15000 a[9] 5.3 0.2 5.0 5.2 5.3 5.4 5.6 1 15000 a[10] 5.1 0.2 4.7 4.9 5.1 5.2 5.4 1 9600 tau.e 11.3 3.1 6.0 9.0 11.0 13.1 18.3 1 11000 tau.a 5.0 2.0 1.9 3.5 4.8 6.2 9.8 1 15000 mu 5.2 0.2 4.9 5.1 5.2 5.3 5.5 1 15000 var.e 0.1 0.0 0.1 0.1 0.1 0.1 0.2 1 11000 var.a 0.2 0.1 0.1 0.2 0.2 0.3 0.5 1 15000 rho 0.7 0.1 0.5 0.6 0.7 0.8 0.9 1 12000 deviance -4.8 8.6 -20.6 -10.9 -5.3 0.8 12.9 1 11000 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = Dbar-Dhat) pD = 10.0 and DIC = 5.2 DIC is an estimate of expected predictive error (lower deviance is better). cust P(rho<=cust) [1,] 0.05815 0.00000000 [2,] 0.09535 0.00000000 [3,] 0.25615 0.00000000 [4,] 0.50860 0.04713333 [5,] 0.76960 0.76606667 [6,] 0.91135 0.99733333 # MM10210T07<-rhohierarchicalmodel(subs="null",c(0.05815,0.09535,0.25615,0.5086,0.7696,0.91135),"sample10",10,3,2,10); Read 30 items Inference for Bugs model at "rhobetamodel.bug", fit using WinBUGS, 3 chains, each with 5e+05 iterations (first 250000 discarded), n.thin = 50 n.sims = 15000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff a[1] 5.0 0.2 4.7 4.9 5.0 5.1 5.3 1 15000 a[2] 5.0 0.2 4.7 4.9 5.0 5.1 5.3 1 15000 a[3] 5.2 0.2 4.9 5.1 5.2 5.3 5.5 1 15000 a[4] 5.1 0.2 4.8 5.0 5.1 5.2 5.4 1 9900 a[5] 5.2 0.2 4.9 5.1 5.2 5.3 5.5 1 14000 a[6] 5.5 0.2 5.2 5.4 5.6 5.7 5.9 1 15000 a[7] 5.3 0.2 5.0 5.2 5.3 5.4 5.6 1 15000 a[8] 5.4 0.2 5.0 5.3 5.4 5.5 5.7 1 7300 a[9] 5.3 0.2 5.0 5.2 5.3 5.4 5.6 1 5200 a[10] 5.1 0.2 4.8 5.0 5.1 5.2 5.4 1 15000 tau.e 11.5 3.2 6.1 9.2 11.2 13.4 18.6 1 12000 tau.a 11.9 3.2 6.4 9.5 11.5 13.9 19.0 1 15000 mu 5.2 0.1 5.0 5.1 5.2 5.3 5.4 1 15000 var.e 0.1 0.0 0.1 0.1 0.1 0.1 0.2 1 12000 var.a 0.1 0.0 0.1 0.1 0.1 0.1 0.2 1 15000 rho 0.5 0.1 0.3 0.4 0.5 0.6 0.7 1 15000 deviance -5.8 8.4 -21.2 -11.7 -6.1 -0.4 11.5 1 14000 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = Dbar-Dhat) pD = 8.7 and DIC = 2.9 DIC is an estimate of expected predictive error (lower deviance is better). cust P(rho<=cust) [1,] 0.05815 0.000000000 [2,] 0.09535 0.000000000 [3,] 0.25615 0.005533333 [4,] 0.50860 0.572266667 [5,] 0.76960 0.999266667 [6,] 0.91135 1.000000000 # MM10102T07<-rhohierarchicalmodel(subs="null",c(0.05815,0.09535,0.25615,0.5086,0.7696,0.91135),"sample10",10,3,10,2); Read 30 items Inference for Bugs model at "rhobetamodel.bug", fit using WinBUGS, 3 chains, each with 5e+05 iterations (first 250000 discarded), n.thin = 50 n.sims = 15000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff a[1] 5.0 0.1 4.7 4.9 5.0 5.1 5.2 1 15000 a[2] 5.0 0.1 4.7 4.9 5.0 5.0 5.2 1 15000 a[3] 5.2 0.1 4.9 5.1 5.2 5.2 5.4 1 15000 a[4] 5.1 0.1 4.8 5.0 5.1 5.2 5.3 1 15000 a[5] 5.3 0.1 5.0 5.2 5.3 5.3 5.5 1 15000 a[6] 5.6 0.1 5.4 5.5 5.6 5.7 5.9 1 15000 a[7] 5.3 0.1 5.1 5.2 5.3 5.4 5.6 1 15000 a[8] 5.4 0.1 5.1 5.3 5.4 5.5 5.6 1 13000 a[9] 5.3 0.1 5.0 5.2 5.3 5.4 5.6 1 15000 a[10] 5.0 0.1 4.8 5.0 5.0 5.1 5.3 1 15000 tau.e 18.4 4.1 11.3 15.6 18.1 21.0 27.3 1 15000 tau.a 5.1 2.1 2.0 3.6 4.8 6.3 9.9 1 7400 mu 5.2 0.2 4.9 5.1 5.2 5.3 5.5 1 15000 var.e 0.1 0.0 0.0 0.0 0.1 0.1 0.1 1 15000 var.a 0.2 0.1 0.1 0.2 0.2 0.3 0.5 1 7400 rho 0.8 0.1 0.6 0.7 0.8 0.8 0.9 1 7900 deviance -18.3 7.1 -31.2 -23.3 -18.6 -13.7 -3.3 1 15000 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = Dbar-Dhat) pD = 9.9 and DIC = -8.4 DIC is an estimate of expected predictive error (lower deviance is better). cust P(rho<=cust) [1,] 0.05815 0.000000000 [2,] 0.09535 0.000000000 [3,] 0.25615 0.000000000 [4,] 0.50860 0.001066667 [5,] 0.76960 0.399000000 [6,] 0.91135 0.977466667 >