rjagsr2jagsjags.parallel

R2JAGS output being placed in model text file instead of R object


I have been learning JAGS through Marc Kery's "Introduction to WinBUGS for Ecologists". When adapting the WinBUGS code from here for the first ch 5 example, I have been getting the odd result that my output is saved in my model.txt file instead of my designated R object. Why is this happening and how can I fix it? Here is my progress so far.

setWD("Insert info here")
library(R2jags)
library(mcmcplots)

#### Ch 5: Running WinBUGS from R via R2WinBUGS ####
### 5.2. Data generation
# Generate two samples of body mass measurements of male peregrines
y10 <- rnorm(n = 10, mean = 600, sd = 30) # Sample of 10 birds
y1000 <- rnorm(n = 1000, mean = 600, sd = 30) # Sample of 1000 birds

# Plot data
xlim = c(450, 750)
par(mfrow = c(2,1))
hist(y10, col = 'grey ', xlim = xlim, main = 'Body mass (g) of 10 male peregrines')
hist(y1000, col = 'grey', xlim = xlim, main = ' Body mass (g) of 1000 male peregrines')

### 5.3. Analysis using R
summary(lm(y1000 ~ 1))

### 5.4. Analysis using WinBUGS
sink("model.txt")
cat("
model {

# Priors
 population.mean ~ dunif(0,5000)        # Normal parameterized by precision
 precision <- 1/population.variance # Precision = 1/variance
 population.variance <- population.sd * population.sd
 population.sd ~ dunif(0,100)

# Likelihood
 for(i in 1:nobs){
    mass[i] ~ dnorm(population.mean, precision)
 }
}
",fill=TRUE)
sink()

# Package all the stuff to be handed over to WinBUGS
# Bundle data
jags.data <- list(mass = y1000, nobs = length(y1000))

# Function to generate starting values
inits <- function(){
  list(population.mean = rnorm(1,600),
       population.sd = runif(1, 1, 30))
}

# Parameters to be monitored (= to estimate)
params <- c("population.mean", "population.sd", "population.variance")

# Running model
out <- jags.parallel(data = jags.data,
                 inits = inits,
                 parameters.to.save = params,
                 model.file="model.txt",
                 n.thin = 1,
                 n.chains = 3,
                 n.burnin = 1,
                 n.iter = 1000)

out

Solution

  • I don't see the behaviour that you describe. What probably happened is that you opened a connection to model.txt but never closed it, because an error occurred. Then you opened it again and closed it properly, but the first version was still open and was now the default for output.

    The code you're using is dangerous because of this sort of thing. A much safer way to create model.txt is to use writeLines(), rather than using sink("model.txt") and cat( ... ). For example,

    model <- "model {
    
    # Priors
     population.mean ~ dunif(0,5000)        # Normal parameterized by precision
     precision <- 1/population.variance # Precision = 1/variance
     population.variance <- population.sd * population.sd
     population.sd ~ dunif(0,100)
    
    # Likelihood
     for(i in 1:nobs){
        mass[i] ~ dnorm(population.mean, precision)
     }
    }
    "
    writeLines(model, "model.txt")