I'm running the following code:
Data_Frame <- structure(list(Predictor_Variable = c(20.5273283761926, 11.4039835403673, 5.31423215288669, 12.5192134582903, 15.8487716107629, 7.16369490255602, 3.91227747895755, 8.90797802712768, 14.1797202872112, 4.82832778943703, 14.3634092528373, 18.3354590029921, 1.56595673179254, 15.27424925589, 2.01471757609397, 17.8340036130976, 24.4356162613258, 24.7056286665611, 8.86376699781977, 3.95776731311344, 1.70528281596489, 24.1929906886071, 7.84694050671533, 2.11047385819256, 24.8220994370058, 14.8339046107139, 19.9926545028575, 19.0412578289397, 4.78973534191027, 18.7796145037282, 19.3868586677127, 14.6739382355008, 20.0653701729607, 19.6233123773709, 7.19542927108705, 16.7649424169213, 9.41006114007905, 19.2996966594364, 21.0871973482426, 13.0715743871406, 17.3938042367809, 18.3959823509213, 9.12836091592908, 2.27788038901053, 10.3573848318774, 4.74680115585215, 9.6620995842386, 16.9264123018365, 11.6498990275431, 20.4607627529185, 24.9407043796964, 21.8109163164627, 11.807474057423, 0.530748104210943, 9.57337225554511, 16.5205421217252, 2.07126556197181, 14.2894889519084, 5.47687077778392, 17.44882510975, 14.0162174822763, 11.7680913361255, 1.04056366835721, 1.95715780719183, 21.3338757341262, 19.7388443630189, 18.9714497013483, 19.5131896412931, 14.1728786868043, 10.4211826983374, 0.625250476878136, 17.7159008861054, 17.5371950492263, 3.91660461900756, 20.3449436288793, 1.5704334480688, 13.397057261318, 20.7767494546715, 7.27919469354674, 23.8504881446715, 19.7164187382441, 5.95191353349946, 20.7321806927212, 10.5640658119228, 16.616114001954, 15.9614237083588, 3.23146634036675, 16.6699483001139, 14.8473161039874, 18.7451402307488, 3.54583212174475, 5.97735904157162, 8.23308551334776, 3.3043225936126, 9.93713270290755, 5.86233736248687, 2.55345033365302, 11.350765300449, 7.88409785600379, 16.5159953408875), Response_Variable = c(1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L)), class = "data.frame", row.names = c(NA, -100L))
sink("Generalised Logistic Function Model.txt")
cat("model {
# Priors
Intercept ~ dnorm(0, 0.001)
Slope ~ dnorm(0, 0.001)
Exponent ~ dlnorm(0, 1)
Sigma ~ dlnorm(0, 1)
Tau <- (1 / (Sigma ^ 2))
# Likelihood and Model Fit
for (i in 1:Number_of_Observations) {
Response[i] ~ dnorm(Mean[i], Tau)
Mean[i] <- ((1 + exp(-(Intercept + (Slope * Predictor[i])))) ^ (-Exponent))
Actual_Squared_Residual[i] <- ((Response[i] - Mean[i]) ^ 2)
Simulated_Response[i] ~ dbern(Mean[i])
Simulated_Squared_Residual[i] <- ((Simulated_Response[i] - Mean[i]) ^ 2)
}
Bayesian_p_Value <- step((sum((Simulated_Squared_Residual[]) ^ 2)) / (sum((Actual_Squared_Residual[]) ^ 2)) - 1)
}", fill = T)
sink()
Parameters <- c("Intercept", "Slope", "Exponent", "Bayesian_p_Value")
Data <- list(Predictor = Data_Frame$Predictor_Variable, Response = Data_Frame$Response_Variable, Number_of_Observations = 100)
Initial_Values <- function () {
list()
}
(Output_1 <- R2jags::jags(Data, Initial_Values, Parameters, "Generalised Logistic Function Model.txt", n.chains = 3, n.thin = 1, n.iter = 1000, n.burnin = 100, working.directory = getwd()))
When I run this Bayesian model, I get the following output:
Inference for Bugs model at "Generalised Logistic Function Model.txt", fit using jags,
3 chains, each with 1000 iterations (first 100 discarded)
n.sims = 2700 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
Bayesian_p_Value 0.436 0.496 0.000 0.000 0.000 1.000 1.000 1.003 870
Exponent 1.122 1.002 0.149 0.431 0.811 1.412 3.900 1.189 16
Intercept -3.550 3.022 -11.786 -4.783 -2.785 -1.452 0.148 1.299 14
Slope 0.259 0.124 0.117 0.173 0.221 0.300 0.597 1.197 18
deviance 112.166 2.789 108.793 110.218 111.477 113.362 119.313 1.008 270
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = var(deviance)/2)
pD = 3.9 and DIC = 116.0
DIC is an estimate of expected predictive error (lower deviance is better).
Then when I try using the R2jags::autojags() function, I get the following output:
(Output_2 <- R2jags::autojags(Generalised_Logistic_Function_Model_Output, Rhat = 1.05))
Inference for Bugs model at "Generalised Logistic Function Model.txt", fit using jags,
3 chains, each with 1000 iterations (first 0 discarded)
n.sims = 3000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
Bayesian_p_Value 0.377 0.485 0.000 0.000 0.000 1.000 1.000 1.020 110
Exponent 0.889 0.831 0.037 0.339 0.670 1.186 3.195 1.757 6
Intercept -8.189 12.251 -47.264 -5.649 -3.372 -1.833 -0.135 1.831 6
Slope 0.473 0.560 0.126 0.191 0.241 0.350 2.130 1.871 6
deviance 112.520 2.919 108.919 110.427 111.930 113.896 119.283 1.103 28
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = var(deviance)/2)
pD = 4.0 and DIC = 116.5
DIC is an estimate of expected predictive error (lower deviance is better).
I was under the impression that this function would continue to run until it converges (until Rhat reached 1.1). This isn't the case - the second model actually has higher Rhat values. I've tried playing around with the Rhat argument in the autojags() function but that hasn't helped. I can see that the second model didn't discard any of the initial iterations, but since it started with a model that was already decent (the model I got using the jags() function), it shouldn't matter. What is going on?
This is what the code for autojags() looks like. We'll go through to see how it works:
n.burnin = object$n.iter
n.thin.auto <- max(1, floor((n.iter - n.burnin)/1000))
n.thin <- ifelse(n.thin > n.thin.auto, n.thin, n.thin.auto)
if (any(!class(object) %in% c("rjags", "rjags.parallel")))
stop("model must be a rjags object")
object <- update(object, n.iter = n.iter, n.thin = n.thin,
refresh = refresh, progress.bar = progress.bar, ...)
check <- any(object$BUGSoutput$summary[, "Rhat"] > Rhat)
if (check) {
n.update (which defaults to 2) and while at least one of the R-hat values is bigger than the one indicated in the function call, do what is in the braces. count <- 1
while (check & (count < n.update)) {
object <- update(object, n.iter = n.iter, n.thin = n.thin,
refresh = refresh, progress.bar = progress.bar,
...)
count counter and check convergence again. count <- count + 1
check <- any(object$BUGSoutput$summary[, "Rhat"] >
Rhat)
}
}
return(object)
In your call to autojags(), you're using the default n.iter (1000) and default n.update (2). This means that the function will stop at a maximum of 3000 iterations (1000 for the first run through and 1000 additional for each of the 2 updates) regardless of convergence. Bayesian models are only guaranteed to converge as the number of draws from the posterior approaches infinity. The stopping rules are there so your computer won't try to run for infinite time. It's also worth noting that the checking is not done iteration-by-iteration, but n.iter-by-n.iter chunks.