ANOVA Practice Note the estimated effect per group and compare.

setwd("/Users/miche/Desktop")
library(rjags)
## Warning: package 'rjags' was built under R version 3.4.3
## Loading required package: coda
## Warning: package 'coda' was built under R version 3.4.3
## Linked to JAGS 4.3.0
## Loaded modules: basemod,bugs
model.1 <- 
  "model {
### data model
  for(i in 1:N){
   y[i] ~dnorm(mu[i], tau)
  mu[i] <- b.15 + b.20*lev.20[i] +b.25 *lev.25[i] + 
                           b.30*lev.30[i] +b.35 * lev.35[i] }

###   prior
         b.15 ~dnorm(0,0.0001); ## referent group     
         b.20 ~dnorm(0,0.0001);      
         b.25 ~dnorm(0,0.0001);       
         b.30 ~dnorm(0,0.0001);       
         b.35 ~dnorm(0,0.0001);     
         tau  ~dgamma(1,1)  
### difference in strength between level 3 (25%) and level 4 (30%)
   b.30.25 <- b.30-b.25

### estimated strength in groups (30%)
   strength[1] <- b.15  
   strength[2] <- strength[1]+b.20 
   strength[3] <- strength[1]+b.25 
   strength[4] <- strength[1]+b.30 
   strength[5] <- strength[1]+b.35 
    }
"

 data.1 <- list( N=25,
lev.20 =  c(0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
lev.25 =  c(0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0),
lev.30 =  c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0),
lev.35 =  c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1),
y = c(7,7,15,11,9,12,17,12,18,18, 14,18,18,18,19,19,
         19,25,22,19,23,7,10,11,15,11))

 model_mean <- jags.model(textConnection(model.1), data=data.1, n.adapt=1000)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 25
##    Unobserved stochastic nodes: 6
##    Total graph size: 154
## 
## Initializing model
 test_mean  <- coda.samples(model_mean, c('b.15','b.20','b.25','b.30','b.35', 'strength'), n.iter = 10000)
 summary(test_mean)
## 
## Iterations = 1:10000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##               Mean    SD Naive SE Time-series SE
## b.15         9.885 1.617  0.01617        0.04675
## b.20         5.517 2.325  0.02325        0.05398
## b.25         7.510 2.317  0.02317        0.05349
## b.30        10.902 2.324  0.02324        0.05448
## b.35         3.300 2.316  0.02316        0.05247
## strength[1]  9.885 1.617  0.01617        0.04675
## strength[2] 15.402 1.682  0.01682        0.01597
## strength[3] 17.394 1.677  0.01677        0.01677
## strength[4] 20.786 1.662  0.01662        0.01662
## strength[5] 13.184 1.677  0.01677        0.01677
## 
## 2. Quantiles for each variable:
## 
##                2.5%    25%    50%    75%  97.5%
## b.15         6.7910  8.787  9.866 10.950 13.130
## b.20         0.8499  3.993  5.540  7.051  9.974
## b.25         2.9062  5.982  7.529  9.044 11.977
## b.30         6.2739  9.361 10.914 12.455 15.425
## b.35        -1.4141  1.805  3.320  4.835  7.819
## strength[1]  6.7910  8.787  9.866 10.950 13.130
## strength[2] 12.0891 14.311 15.374 16.498 18.759
## strength[3] 14.0806 16.310 17.391 18.483 20.728
## strength[4] 17.5899 19.691 20.802 21.844 24.075
## strength[5]  9.8558 12.085 13.198 14.266 16.515
 plot(test_mean)

 autocorr.plot(test_mean)

 test_mean  <- coda.samples(model_mean, c('b.15','b.20','b.25','b.30','b.35', 'strength','b.30.25'), n.iter = 20000,thin=5)
 summary(test_mean)
## 
## Iterations = 10005:30000
## Thinning interval = 5 
## Number of chains = 1 
## Sample size per chain = 4000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##               Mean    SD Naive SE Time-series SE
## b.15         9.749 1.691  0.02674        0.03737
## b.20         5.652 2.367  0.03743        0.04887
## b.25         7.695 2.380  0.03763        0.04811
## b.30        11.006 2.351  0.03718        0.04782
## b.30.25      3.311 2.360  0.03732        0.03732
## b.35         3.452 2.372  0.03750        0.04773
## strength[1]  9.749 1.691  0.02674        0.03737
## strength[2] 15.402 1.650  0.02609        0.02609
## strength[3] 17.444 1.673  0.02645        0.02645
## strength[4] 20.755 1.671  0.02642        0.02642
## strength[5] 13.202 1.670  0.02641        0.02701
## 
## 2. Quantiles for each variable:
## 
##                2.5%    25%    50%    75%  97.5%
## b.15         6.3776  8.650  9.747 10.828 13.134
## b.20         0.9824  4.105  5.701  7.183 10.443
## b.25         2.9940  6.118  7.716  9.305 12.471
## b.30         6.3350  9.481 11.009 12.550 15.735
## b.30.25     -1.4043  1.762  3.337  4.850  7.935
## b.35        -1.0961  1.900  3.427  4.960  8.158
## strength[1]  6.3776  8.650  9.747 10.828 13.134
## strength[2] 12.1111 14.321 15.443 16.484 18.601
## strength[3] 14.2809 16.319 17.437 18.550 20.732
## strength[4] 17.4235 19.665 20.754 21.836 24.024
## strength[5]  9.9612 12.101 13.176 14.295 16.511
 plot(test_mean)

 autocorr.plot(test_mean)

 out <- as.matrix(test_mean)
 boxplot(out[,7:11])