# B2prop.r
# written by JuG
# April 19 2021
#' Bayesian comparison of two proportions
#' @author JuG
#' @description Cod to compare two proportions using bayesian model with beta priors
#' @param var1 variable 1
#' @param var2 = variable 2
#' @param n.iter number of iteration
#' @param n.burnin number of iteration used in the update of the model
#' @param n.thin n thin
#' @param prior.beta beta priors (vector of 4 parameters, two first for var1)
#' @details
#' @examples
#' dtf<-data.frame(grp=c(rep("A",15),rep("B", 15)),
#' res=c(rbinom(n = 15, size = 1, prob = .5),
#' rbinom(n = 15, size = 1, prob = .49)) )
#'
#' res <- B2prop(var1 = dtf$res[dtf$grp=="A"], var2 = dtf$res[dtf$grp=="B"])
#' res
#' @return list with Stats (descriptive stats) and MCMC (inference)
#' @export
B2prop<-function(var1,var2,n.iter=11000,n.burnin=1000,n.thin=1,prior.beta=c(1,1,1,1)){
library(rjags)
#--------------------------------------------------------------------------------------------------
mod1_string = " model {
for (i in 1:N1) { X[i] ~ dbin(mu.x,1) }
for (j in 1:N2) { Y[j] ~ dbin(mu.y,1) }
mu.x~dbeta(a,b)
mu.y~dbeta(c,d)
di<-mu.x-mu.y
p<-step(mu.x-mu.y)
or<-mu.x*(1-mu.y)/(mu.y*(1-mu.x))
}
"
donnees<-list( N1 =length(var1),
N2 =length(var2),
X =as.numeric(var1),
Y =as.numeric(var2),
a =prior.beta[1],
b =prior.beta[2],
c =prior.beta[3],
d =prior.beta[4] )
parametre <-c("mu.x","mu.y","di","p","or")
inits1 = function() {
inits = list("mu.x"= .5,
"mu.y" = .5)
}
mod1 = jags.model(textConnection(mod1_string), data=donnees, inits = inits1, n.chains=3, quiet = T)
update(mod1, n.burnin) # burn-in
mod1_sim = coda.samples(model=mod1,
variable.names=parametre,
n.iter=n.iter, thin = n.thin)
mod1_csim = do.call(rbind, mod1_sim)
mod1_csim2 <- mod1_csim
nom1 <-deparse(substitute(var1)) ### Autres STATS a faire apparaitre ?galement
nom2 <-deparse(substitute(var2))
indexDiff <-grep ("di",colnames(mod1_csim2))
indexMux <-grep ("mu.x",colnames(mod1_csim2))
indexMuy <-grep ("mu.y",colnames(mod1_csim2))
indexProb <-grep ("p",colnames(mod1_csim2))
indexOR <-grep ("or",colnames(mod1_csim2))
RES<-mod1_csim2[,c(indexMux,indexMuy,indexDiff,indexProb,indexOR)]
moyennes <-round(apply(RES,2,mean,na.rm=T) ,3)
variances <-round(apply(RES,2,var,na.rm=T),3)
sd <-round(apply(RES,2,sd,na.rm=T),2)
q0.025 <-round(apply(RES,2,quantile,na.rm=T, prob=.025),2)
q0.250 <-round(apply(RES,2,quantile,na.rm=T, prob=.25),2)
q0.500 <-round(apply(RES,2,quantile,na.rm=T, prob=.5),2)
q0.750 <-round(apply(RES,2,quantile,na.rm=T, prob=.75),2)
q0.975 <-round(apply(RES,2,quantile,na.rm=T, prob=.975),2)
RES <-cbind(moyennes,variances,sd,q0.025,q0.250,q0.500,q0.750,q0.975)
rownames(RES) <-c( paste("Prop",nom1),paste("Prop",nom2),"difference",paste("P(",nom1,">",nom2,")")[1],"Rapport_cotes")
diff.signif<-c("",""," ***","","")
if(RES[3,4]<0 & RES[3,8]>0){diff.signif<-c("","","","","")}
RES<-noquote(cbind(RES,diff.signif))
colnames(RES) <-c("Moy","Var","Sd","2.5%","25%","50%","75%","97.5%","")
RES<-RES[,-2]
var11<-var1
var22<-var2
var1<-var1[!is.na(var1)]
var2<-var2[!is.na(var2)]
succes1 <-sum(var1,na.rm=TRUE)
echec1 <-length(var1)-succes1
succes2 <-sum(var2,na.rm=TRUE)
echec2 <-length(var2)-succes2
p.succes1 <-round(succes1/length(var1),2)
p.succes2 <-round(succes2/length(var2),2)
p.echec1 <-round(echec1/length(var1),2)
p.echec2 <-round(echec2/length(var2),2)
NA1 <-round(sum(is.na(var11)),2)
NA2 <-round(sum(is.na(var22)),2)
pNA1 <-round(NA1/length(var11),2)*100
pNA2 <-round(NA2/length(var22),2)*100
N1 <-length(var11)
N2 <-length(var22)
RES1<-noquote(cbind(c(N1,N2),c(NA1,NA2),c(pNA1,pNA2),c(succes1,succes2),c(echec1,echec2),c(p.succes1,p.succes2),c(p.echec1,p.echec2)))
colnames(RES1) <-c("N","NA","%NA","#1","#0","%1","%0")
rownames(RES1) <-c(nom1,nom2)
RESULTAT <-list(Stats=RES1,MCMC=RES)
return(RESULTAT)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.