0

I'm trying to run this script in R2jags following the instructions provided in "Lahoz-Monfort JJ, Guillera-Arroita G, Tingley R (2015) Statistical approaches to account for false positive errors in environmental DNA samples. Molecular Ecology Resources, 16, 673–685." and it seems that it worked ok, but I can't figure out the command to see the results... could anyone please help?

cat("model {
    # Priors
    psi ~ dunif(0,1)
    p11 ~ dunif(0,1)
    p10 ~ dunif(0,p10_max)
    
    # Likelihood 
    for (i in 1:S){
    z[i] ~ dbern(psi)
    p[i] <- z[i]*p11 + (1-z[i])*p10
    for (j in 1:K){
    Y[i,j] ~ dbern(p[i])
    }
    }
    } ",fill=TRUE)
sink()


Bayesian <- function(psi,p11,p10,S,K,nsims=100,doprint=TRUE,p10_max=0.05,
                                     ni=100000,nt=2,nc=1,nb=50000,myparallel=TRUE) {   
  psihat<-p11hat<-p10hat<-rep(nsims)
  modelSummaries<-list()
  for(ii in 1:nsims){
    if (doprint) cat("\r", ii, "of", nsims,"   ")    
    hh<-genSimData(psi,r11=0,p11,p10,S,K1=0,K2=K)
    # fit the model    
    jags.inits <-function()(list(psi=runif(1,0.05,0.95),p11=runif(1,p10_max,1),p10=runif(1,0,p10_max)))
    jags.data  <-list(Y=hh,S=S,K=K,p10_max=p10_max)
    jags.params<-c("psi","p11","p10")
    Thoropa_model<-jags(data = jags.data, inits = jags.inits, parameters.to.save= jags.params, 
                model.file= "Thoropa.txt", n.thin= nt, n.chains= nc, 
                n.iter= ni, n.burnin = nb, parallel=myparallel)  #, working.directory= getwd()
    # extract results (medians of the marginal posteriors)
    psihat[ii] <- model$summary["psi","50%"]
    p11hat[ii] <- model$summary["p11","50%"]
    p10hat[ii] <- model$summary["p10","50%"]    
    modelSummaries[[ii]]<-model$summary
  }
  if (doprint){
    printsummres(psihat,thename="estimated psi")
    printsummres(p11hat,thename="estimated p11")
    printsummres(p10hat,thename="estimated p10")
  }
  return(list(psihat=psihat,p11hat=p11hat,p10hat=p10hat,modelSummaries=modelSummaries))
  }

The file "Thoropa.txt" is a presence/absence matrix as follows:

PCR1    PCR2    PCR3    PCR4    PCR5    PCR6    PCR7    PCR8    PCR9    PCR10   PCR11   PCR12
1   0   0   0   0   0   0   0   0   0   0   0
0   0   0   1   0   0   0   1   0   0   0   0
0   0   0   0   0   0   1   0   0   0   0   0
1   1   1   1   1   1   1   1   1   1   1   1
0   0   0   0   0   0   0   0   0   0   0   0
1   0   1   0   1   1   1   1   1   1   1   1
0   0   1   0   0   0   1   0   0   0   0   0
1   0   1   0   0   0   0   0   1   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0
1   0   0   0   1   1   0   0   0   1   0   0
1   1   0   1   0   1   0   1   0   0   1   0
1   1   0   0   0   0   0   1   0   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0
1   1   1   1   0   1   1   1   1   0   1   1
1   1   1   1   1   1   1   1   1   1   1   1
0   0   1   0   1   0   0   0   0   0   1   1
1   1   1   1   1   1   1   1   1   1   1   1

Following the comment from Limey (thank you!) I changed the script to:


sink("Thoropa2.txt")
cat("model {
    # Priors
    psi ~ dunif(0,1)
    p11 ~ dunif(0,1)
    p10 ~ dunif(0,p10_max)
    
    # Likelihood 
    for (i in 1:S){
    z[i] ~ dbern(psi)
    p[i] <- z[i]*p11 + (1-z[i])*p10
    for (j in 1:K){
    Y[i,j] ~ dbern(p[i])
    }
    }
    } ",fill=TRUE)
sink()

y=Thoropa# the detection/non detection table
S=nrow(y)
K=ncol(y)

psi ~ dunif(0, 1)
p11 ~ dunif(0, 1)
p10 ~ dunif(0, p10_max)
p10_max=0.05

jags.data<-list(y=y, S=S, K=K, p10_max=p10_max)

jags.inits <-function()(list(psi=runif(0,1),p11=runif(0,1),p10=runif(0,p10_max)))

jags.params<-c("psi","p11","p10")

Thoropa_model<-jags.parallel(data = jags.data, inits = jags.inits, parameters.to.save= jags.params, model.file= "Thoropa2.txt", n.chains= 4, n.thin= 10, n.iter = 100000, n.burnin=50000, jags.seed = 333)

and the data file is as before. Now I am getting the error message:

"Error in checkForRemoteErrors(val) : 4 nodes produced errors; first error: Indexing outside the bounds"

Could anyone help identify the error in my script? I'm no expert and I'm learning by myself, so sorry if it is a stupid question... (maybe there is something wrong with the format of the data...?)

Thank you all!

Paula Cabral Eterovick
  • 23
  • 5
  • So are you calling the `Bayesian()` function somewhere? Or how exactly are you running this? – MrFlick Oct 02 '21 at 23:23
  • No, I just used "Bayesian" as a name, there is nothing else in the script other than what I pasted in the message... I used the function jags() for the model – Paula Cabral Eterovick Oct 02 '21 at 23:31
  • 1. How do you expect the `Bayesian` function to execute if you don't call it? 2. `jags` expects your model file to be passed in its `model.file` argument and the observed data in its `data` parameter. You have them the other way round. I'm surprised you don't have any errors reported. – Limey Oct 03 '21 at 07:26
  • @Limey I believe I fixed it as you said (yes, it was also surprising that I got no error messages at all... I also identified other mistakes). I did the following script but now I am getting the message "4 nodes produced errors; first error: indexing outside the bounds" - would you be so nice to suggest how I can fix this? I'm no expert and I'm learning by myself, so sorry if it is a stupid question... (I will edit the question to add the new script and data file - maybe there is something wrong with the format of the data...?) – Paula Cabral Eterovick Oct 05 '21 at 23:17

1 Answers1

0

Your model is not working because some syntax errors in your R script. Note that the R syntax is different of the jags syntax, even if you are running the jags inside de R. These are the errors:

  1. The symbol "~" is not used for sampling in R. Delete the lines:

    psi ~ dunif(0, 1)
p11 ~ dunif(0, 1)
p10 ~ dunif(0, p10_max)

  2. The Y variable in the jags model is capitalized, so you must correct the syntax in jags.data.

    jags.data<-list(Y=y, S=S, K=K, p10_max=p10_max)

  3. In jags.inits, a) the body of the function must be inside of curly braces, and not parentheses, and b) the function runif takes 3 arguments: n (number of values you want to sample), min and max. The correct syntax is the following:

     jags.inits <-function(){list(psi=runif(1,0,1),p11=runif(1,0,1),p10=runif(1,0,p10_max))}

Fixing those errors, your model should run without errors. After run the model, you can extract the median of the parameters "psi" using one of these two options:

   Thoropa_model$BUGSoutput$median$psi

   Thoropa_model$BUGSoutput$summary["psi","50%"]
Renato Fernandes
  • 107
  • 5