####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
>