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