Dear All,
how I can translate this winBug model for studying the treatment effect to rstanarm interface, using rstanarm package taking into account all the parameters defined below:
and thanks in advance:
> mydata = data.frame(treat = c( rep(0, 120), rep(1, 150), rep(0, 40), rep(1, 30)),
+ death = c( rep(0, 120), rep(0, 150), rep(1, 40), rep(1, 30) ) )
> # nc-yc nt-yt yc # yt
>
> str(mydata)
'data.frame': 340 obs. of 2 variables:
$ treat: num 0 0 0 0 0 0 0 0 0 0 ...
$ death: num 0 0 0 0 0 0 0 0 0 0 ...
>
>
> winBugModel <- function(){
+ yt ~ dbin(p1, nt)
+ yc ~ dbin(pc, nc)
+ pc ~ dunif(0,1)
+ eta <- log(pc)
+ log(p1) <- eta + coef
+ coef ~ dnorm(meanPrior, sigmaPrior)
+ RR <- exp(coef)
+ }
>
> initial <- function(){
+ list(pc = runif(1,0,1),
+ coef = runif(1, -2,0))
+ }
>
>
>
> dataTrail <- list(yt = 30,
+ nt = 180,
+ yc = 40,
+ nc = 160,
+ meanPrior = 0,
+ sigmaPrior = 1/200
+ )
>
> output <- jags(data = dataTrail,
+ inits = initial,
+ parameters.to.save = c("RR"),
+ model.file = winBugModel,
+ n.chains = 4,
+ n.iter = 20000,
+ n.thin = 2)
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 2
Unobserved stochastic nodes: 2
Total graph size: 14
Initializing model
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|**************************************************| 100%
>
>
> OUTPUTS <- rbind(output$BUGSoutput$sims.array[,1,],
+ output$BUGSoutput$sims.array[,2,],
+ output$BUGSoutput$sims.array[,3,],
+ output$BUGSoutput$sims.array[,4,])
>
> RR = mean(OUTPUTS[,1])
>
> c(RR = mean(OUTPUTS[,1]), quantile(OUTPUTS[,1], 0.025), quantile(OUTPUTS[,1],0.975))
RR 2.5% 97.5%
0.6669388 0.4223434 0.9914351
>
> 100*mean( (1-OUTPUTS[,1])>0 )
[1] 97.675
>
> 100*mean( (0.5-OUTPUTS[,1])>0 )
[1] 11.215