Solved – Bayesian parameter estimation with negative binomial data. No model fit achieved!

I edited this question to a (hopefully) less confusing structure. The full R code is still available at the end of the post.

My problem is a bad fit of the parameter n to my data. I will try to shortly explain the steps I have already taken on solving the problem.

My question is, do you suggest that this function I am using is unsuitable for my data, or do I miss anything at the implementation of the model?

The rjags model I am working on is this: model {  for (i in 1:N) {  y[i] ~ dnegbin (p[i], r) p[i] <- r / (r + mu[i])  mu[i] <- b * n * (1 - exp(-x[i] / n))   # biological meaningful function    }  r ~ dunif(0,500)  logit(b) <- logit.b logit.b ~ dnorm(0.5, 0.00001)  log(n) <- log.n     # prior for n, parameterized as a normal distribution log.n ~ dnorm(0, 0.00001)       # #n ~ dpois(30000)         # alternative prior with poisson distribution  } 

I run the sampler with 3 chains, 10000 adapt steps, 50000 burn-in steps and 300000 iterations.

I am under the impression, that I am already using the reparameterization by parameters m and r, as J. Kruschke suggest in his blog post from 7th of April 2012. J. Kruschke blog

I also tried to get a better prior on n:

· I narrowed the mean and variance parameters of the dnorm distribution on the prior for n.

· I used a truncated normal distribution

· I tried a poisson distribution

Every time I narrowed the parameter space below a critical point, the posterior distribution of n was highly dependent on the prior. Every time I used an uninformed, broad prior like
log.n ~ dnorm(0,0.00001), the posterior distribution becomes bimodal and contains huge values.

The function I am using has a biological meaning. The b parameter is a proportion of suitable habitats and the n parameter the total number of habitats. It is obvious that the parameter n should not expand to incredibly high numbers.

This is the old, complete text, with some more details and the data

I got 36 observations of count data. The data is negativ binomial distributed.
I am using rjags so I used the JAGS code.
The negbin parameter p is derived from the mean of the model (mu) I want to fit to the data and r, which has a uniform prior.
For the model (mu) the parameters b and n use transformed values (logit for b, its a rate) and n (log) to match the normal prior distributions I choose.

The model I choose has a certain meaning in my research context, so I would avoid changing the model, if possible.

My problem is a bad fit of the parameter n!
As you can see in the last plot, the estimated model line does not fit the data at all. I am not shure if I am doing it right by looking at the posterior distributions and HPD for each parameter on it's own. I am under the impression that there has to be a "best fit" joint probability for both parameter together?

library(rjags) library(coda) library(dplyr) library(ggplot2) library(lattice) library(faraway) setwd("R:/desktop") y<-c(129,1,441,649,121,1,337,5,489,209,1,181,237,273,181,5,1,209,9,221,273,9,77,325,473,269,1,177,233,21,209,93,1,357,89,1) x <- c(310,10,1210,2410,610,10,1210,10,2410,310,10,610,1210,2410,610,10,10,310,10,610,2410,10,310,1210,2410,1210,10,310,610,10,1210,610,10,2410,310,10) #--------------------------------------------------------------- modelString =" model { for (i in 1:N) {  y[i] ~ dnegbin (p[i], r)  # Likelihood function for data(y), 0 < p < 1,  r natural number without zero, density: ((x+r-1)/x)*(p^r)*(1-p)^x  p[i] <- r / (r + mu[i])  # parametarisation for p[i] (probability of single-success in negbin) mu[i] <- b * n * (1 - exp(-x[i] / n))  # model function, mu (mean) parameter of p[i]  }  r ~ dunif(0, 500)  # uninformed prior for parameter r (number of success until abort) of the Likelihood function ( a < b; lower = a, upper = b), density: 1/(b-a) logit(b) <- logit.b logit.b ~ dnorm(0, 0.00001)  # prior distribution for parameter b, dnorm(mu,tau) tau is the precission = (1/(sd)²) log(n) <- log.n log.n ~ dnorm(0, 0.00001)  # prior distribution for n  } " # close quote for model string  # write modelString to a text file writeLines( modelString, con = "model1.txt")  #--------------------------------------------------------------  # Set MCMC chain parameters  dataList <- list( y = y , x = x ,  N = length(y)) initsList <- function() list(logit.b = rnorm(1), log.n = rnorm(1), r =runif(1))  # initial values for each parameter and chain. Can be omitted  nChains <- 3  # Number of chains to run the MCMC  adaptSteps <- 10000  # make dependet on number of iterations!!!!!!  burnInSteps <- 50000  # independet of the adaptive phase, but not saved!!!!!!  numSavedSteps <- 300000  # number of saved steps to the coda.samples MCMC output  thinSteps <- 1  # thinning interval for monitor  nIter <- ceiling( (numSavedSteps * thinSteps) / nChains)  # number of iterations to monitor  #--------------------------------------------------------------  # Run the chains  parameters = c("logit.b", "log.n", "r")  # parameter to be monitored  # using the jags.model function to run the chains  jagsModel = jags.model("model1.txt", data = dataList, inits = initsList, n.chains = nChains, n.adapt = adaptSteps)  # Burn-in: (this is a number of not-saved itereations, independet of the adaptive phase)  cat("Burning in the MCMC chain")  # text output update(jagsModel, n.iter=burnInSteps)  # n.iter is set above  # Now the saved MCMC chain is running  cat("Sampling final MCMC chain")  # text output codaSamples = coda.samples(jagsModel, variable.names = parameters, n.iter = nIter)  #-------------------------------------------------------------- # #Backtransformation of the mcmc.list object  codaResults <- codaSamples for(i in 1:3){ codaResults[[i]][,1] <-    exp(codaResults[[i]][,1]) codaResults[[i]][,2] <- ilogit(codaResults[[i]][,2]) dimnames(codaResults[[i]])[[2]][1] <- "n" dimnames(codaResults[[i]])[[2]][2] <- "b" } coda_sum <- summary(codaResults)  #-------------------------------------------------------------  #Plot results  gelman.plot(codaResults)  densplot(codaResults)   data_plot <- data.frame(y,x) xx <- 1:2600 b1 = coda_sum$statistics[2,1]       #mean for parameter b n1 = coda_sum$statistics[1,1]       #mean for parameter n curveP <- data.frame(xx = xx, y1 = b1 * n1 * (1 - exp(-xx / n1)))  plot_P <- ggplot(data_plot, aes(x = x)) + ggtitle(" curve fit to count data") +              geom_point(aes(y = y), colour = "black", size = 1.5, shape = 21, fill ="black") + geom_line(data = curveP, aes(y = y1, x = xx), lwd= 0.7, lty= 1) + labs(x = "density", y= "numbers") + xlim(0,2600) + ylim(0,1300) +            #(0,1150) oder (0,200) theme_bw() + theme(axis.text=element_text(size=8),     plot.title = element_text(color="black", size=10, face="plain"),     axis.title=element_text(size=10,face="plain"),     panel.grid.major = element_blank(),     panel.grid.minor = element_blank(),     panel.border = element_blank(),     panel.background = element_blank())+  theme(axis.line.x = element_line(color="black", size = .5),     axis.line.y = element_line(color="black", size = .5)) plot_P 

PS: I have in total 24 of those kind of data sets. This one is not the only one with such a confusing outcome.

# Original script from http://stats.stackexchange.com/questions/257742/bayesian-parameter-estimation-with-negative-binomial-data-no-model-fit-achieved # Extensively modified by John K. Kruschke, 10-March-2017.  graphics.off() # This closes all of R's graphics windows. rm(list=ls())  # Careful! This clears all of R's memory!  library(rjags)  y <- c(129,1,441,649,121,1,337,5,489,209,1,181,237,273,181,5,1,209,9,221,273,9,77,325,473,269,1,177,233,21,209,93,1,357,89,1)  x <- c(310,10,1210,2410,610,10,1210,10,2410,310,10,610,1210,2410,610,10,10,310,10,610,2410,10,310,1210,2410,1210,10,310,610,10,1210,610,10,2410,310,10)  #----------------------------------------------------------------------------- # Find Poisson MLE fit: fitTrend = function( paramVec , y=y , x=x ) {   b = paramVec[1]   n = paramVec[2]   mu = b * n * ( 1 - exp( -x / n ) )     neglogpy = -sum( dpois( y , mu , log=TRUE ) )   return( neglogpy ) } b = 1.0 ; n = 200 # initial values for search optimInfo = optim( c(b,n) , fitTrend , x=x , y=y ) bOpt = optimInfo$par[1] nOpt = optimInfo$par[2] # Plot data with best fitting trend: plot( y~x ) xcomb = seq( min(x) , max(x) , length=501 ) lines( xcomb , bOpt * nOpt * ( 1 - exp( -xcomb / nOpt ) )  , col="skyblue" ) title( main=bquote(list(b==.(bOpt),n==.(nOpt))) ) # Implication: Make sure that prior on b,n is appropriate for best fitting # values. My understanding is that both b and n must be positive, so I'll use # gamma priors with modes near the best fit values, and wide standard # deviations. # Function for finding dgamma parameters from DBDA2E-utilities.R, see # https://sites.google.com/site/doingbayesiandataanalysis/software-installation gammaShRaFromModeSD = function( mode , sd ) {   if ( mode <=0 ) stop("mode must be > 0")   if ( sd <=0 ) stop("sd must be > 0")   rate = ( mode + sqrt( mode^2 + 4 * sd^2 ) ) / ( 2 * sd^2 )   shape = 1 + mode * rate   return( list( shape=shape , rate=rate ) ) } nGammaPrior = gammaShRaFromModeSD( mode=nOpt , sd= nOpt/2 ) nSh = nGammaPrior$shape nRa = nGammaPrior$rate bGammaPrior = gammaShRaFromModeSD( mode=bOpt , sd= bOpt/2 ) bSh = bGammaPrior$shape bRa = bGammaPrior$rate #--------------------------------------------------------------- modelString ="   model {   for (i in 1:N) {     # Likelihood      y[i] ~ dnegbin( p[i] , r )     # re-parametarization for p[i]: See     # http://doingbayesiandataanalysis.blogspot.com/2012/04/negative-binomial-reparameterization.html     p[i] <- r / ( r + mu[i] )       # Trend function:     mu[i] <- b * n * ( 1 - exp( -x[i] / n ) )     }   # Prior:   r ~ dgamma( 0.1 , 0.1 )   b ~ dgamma( bSh , bRa )   n ~ dgamma( nSh , nRa ) } " # close quote for model string # write modelString to a text file writeLines( modelString, con = "model1.txt") #-------------------------------------------------------------- dataList <- list(   y = y ,   x = x ,    N = length(y),   bSh = bSh ,   bRa = bRa ,   nSh = nSh ,   nRa = nRa ) #-------------------------------------------------------------- nChains <- 3  # Number of chains to run the MCMC adaptSteps <- 1000  # make dependet on number of iterations!!!!!! burnInSteps <- 1000  # independet of the adaptive phase, but not saved!!!!!! numSavedSteps <- 10000  # number of saved steps to the coda.samples MCMC output thinSteps <- 1  # thinning interval for monitor nIter <- ceiling( (numSavedSteps * thinSteps) / nChains)  # number of iterations to monitor parameters = c("b", "n", "r")  # parameter to be monitored jagsModel = jags.model( "model1.txt", data=dataList,                          n.chains=nChains, n.adapt=adaptSteps ) cat("Burning in the MCMC chain")  # text output update(jagsModel, n.iter=burnInSteps)  # n.iter is set above cat("Sampling final MCMC chain")  # text output codaSamples = coda.samples(jagsModel, variable.names = parameters, n.iter = nIter) mcmcMat = as.matrix( codaSamples )  # IMPORTANT: Inspect chain diagnostics here... # **  #----------------------- # Plot fit over data: plot( x , y ) xcomb = seq( min(x) , max(x) , length=501 ) for ( stepIdx in floor( seq(1,nrow(mcmcMat),length=20) ) ) {   lines( xcomb ,           ( mcmcMat[stepIdx,"b"] * mcmcMat[stepIdx,"n"]           * ( 1 - exp( -xcomb / mcmcMat[stepIdx,"n"] ) ) ) ,          col="skyblue" ) } #-------------------------------------------------------------- 

Similar Posts:

• Solved – Estimation of the shape parameter of gamma distribution in Winbugs
• Solved – Combining posterior distributions
• Solved – Combining posterior distributions
• Solved – Using R2WinBUGS, how to extract information from each chain
• Solved – Using R2WinBUGS, how to extract information from each chain