0

I'm trying to perform 2-way bayesian ANOVA with jags, but there is an error I can not understand.

## data
set.seed(123)
n <- 30 
y <- log(rnorm(n, 3, 1))
x1 <- as.numeric(c(1, 2, 1, 2, 2, 1, 1, 1, 2, 1, 2, 1, 2, 2, 1, 2, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 2, 1, 1)) 
x2 <- as.numeric(c(2, 1, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 1, 2, 1, 2, 1, 1, 1, 1, 2, 1, 1, 2, 2, 1, 1, 1))
                                  
##
                                  
data <- list(n = length(y), 
             g = length(table(x1)),
             h = length(table(x2)),
             x1 = as.numeric(x1),
             x2 = as.numeric(x2),
             y = y)
                                  
## 
mean.x1 <- as.numeric(tapply(y,x1,mean))
mean.x2 <- as.numeric(tapply(y,x2,mean))
                                  
inits <- list(list( x1 = mean.x1, x2 = mean.x2, tau = c(2, 3)),
              list( x1 = mean.x1+10, x2 = mean.x2+10, tau = c(2, 3)*0.5),
              list( x1 = mean.x1-15, x2 = mean.x2-15, tau = c(2, 3)*2))
                                  
##
sink("model.txt")
cat('
model{
base ~ dnorm(0, 1.0E-6)
                                         
for(i in 1:g) {
x1[i] ~ dnorm(0.0, 1.0E-6)
}
                                         
for(i in 1:h) {
x2[i] ~ dnorm(0.0, 1.0E-6)
}
                                         
tau ~ dgamma(0.001, 0.001) 
sigma <- 1/sqrt(tau)
                                         
for (i in 1:n) {
mean[i] <- base + x1[g[i]] + x2[h[i]]
y[i] ~ dnorm(mean[i], tau)
}
}
')
sink()


##
library(jagsUI)
mcmc <- jags(data = data,
            inits = initis,
            parameters.to.save = c("x1","x2","sigma"), 
            model.file = "model.txt",
            n.chains = 3,n.adapt = 1000,
            n.iter = 12000,n.burnin = 2000,n.thin = 5)

Error in jags.model(file = model.file, data = data, inits = inits, n.chains = n.chains, : RUNTIME ERROR: Compilation error on line 17. Index out of range taking subset of g

Also, I would like to add an interaction term.

Any answer would be helpful. Thank you.

Lucas Rodrigues
  • 1
  • 2
  • Hi Lucas, welcome to Stack Overflow. In order to be able to help, we need a fully reproducible example of your code. As you will note in the error, it says the subset of `g` is out of range. In line 16, you subset `g[i]` with `i` in a loop between `1` and `n`. However `n` is not defined in the code you have provided. Could this be your problem? Otherwise, see [How to make a reproducible example](https://stackoverflow.com/questions/5963269/) for more info. – Ian Campbell Jul 12 '20 at 04:25
  • An ANOVA is just a linear model with a categorical predictor. I would just use `brms::brm()`. – Phil Jul 12 '20 at 04:57
  • Hi Ian thanks for the quick and cordial answer. I am dealing with bacterial community, based on metabarcoding. In that community there is approximately one thousand OTUs in 30 samples, which I could calculate Shannon's diversity index for each sample. I update the code, and I hope it can be useful now. – Lucas Rodrigues Jul 12 '20 at 14:41

0 Answers0