How to use dround

[Jonas Lindeløv]

Although the dround "distribution" seems conceptually simple, I get an initialization error message that makes me wonder if I'm using it right. I'll illustrate with a minimal example: given an observed integer x we want to infer the real value truescore which is know to be in the interval 0-100.

    truescore ~ dunif(0, 100)
    x ~ dround(truescore, 0)

This works for x=50 (or more generally for x in the middle of dunif) but not for any other values (e.g. 51, 2, 65 etc.). It says Error in node x, Observed node inconsistent with unobserved parents at initialization. It seems strange that initialization values is necessary in such a simple case. Am I doing it wrong?

Furthermore, for x=50, the posterior on "truescore" behaves as if the observation was rounded to nearest integer for digits parameter below 2. (E.g. for 1, 0, -1, -100 etc. the posterior is x +/- 0.5). I would expect the posterior to broaden as the measurement resolution decreases and for digits=1 to give an estimate of +/- 0.05 rather than +/- 0.5.

Best,
Jonas

[Gwern Branwen]

I think I might have run into the same or maybe a similar issue here, using JAGS 4.2.0.

For my case, I have a normally distributed data variable 1-5, mean 3, SD 0.5 or so, but it is rounded to 1/2/3/4/5. I thought it would be nice to include the rounding and saw that the JAGS manual says:

Rounded data: dround
The dround distribution represents rounded data. It has two parameters: t the original continuous variable and d, the number of decimal places to which the measurements are rounded. Thus if t = 1.2345 and d = 2 then the rounded value is 1.23. Note that d can be negative: if d = −2 then the data are rounded to the nearest 100.

There aren't any examples in the manual or online anywhere that I could find, but my assumption is that I should be able to just do this normally as dround(x[i], 0) etc. But this doesn't seem to work for me, even when an almost identical use of dinterval() does:

data_rounded <- round(rnorm(100, mean=3, sd=0.5)); data_rounded
#  [1] 2 3 3 3 3 3 3 3 4 4 4 3 2 3 4 3 2 3 3 4 3 3 2 3 3 2 3 3 2 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 2 2 3 3 3 2 3 3 3 3 2 2 4 3 3 3 3 3 3 4 4 3 2 2 3 3 3 4 2 2 3 3 3 3 3 3 3 3 2 3 3 3 3 2
# [93] 3 4 4 3 2 3 2 3

## Model treating it as rounded:
library(runjags)
model1 <- "model {
    for (i in 1:n) {
           y[i] ~ dround(true[i], 1)
           true[i] ~ dnorm(mu, tau)
       }

    mu ~ dnorm(3, pow(1, -2))
    tau ~ dgamma(0.001, 0.001)
    }"
data <- list(y=data_rounded, n=length(data_rounded))
params <- c("mu", "tau")
j1 <- run.jags(model=model1, monitor=params, data=data); j1
# Compiling rjags model...
# Error: The following error occured when compiling and adapting the model using rjags:
#  Error in rjags::jags.model(model, data = dataenv, n.chains = length(runjags.object$end.state),  : 
#   Error in node y[1]
# Node inconsistent with parents
# 
# 
# 
# It may help to use failed.jags(c('model','data','inits')) to see model/data/inits syntax with line numbers
# In addition: Warning messages:
# 1: No initial value blocks found and n.chains not specified: 2 chains were used 
# 2: No initial values were provided - JAGS will use the same initial values for all chains 
failed.jags(c('model','data','inits'))
# 
# JAGS model syntax:
# 
# 1  |  model{
# 2  |  for (i in 1:n) {
# 3  |  y[i] ~ dround(true[i], 1)
# 4  |  true[i] ~ dnorm(mu, tau)
# 5  |  }
# 6  |  mu ~ dnorm(3, pow(1, -2))
# 7  |  tau ~ dgamma(0.001, 0.001)
# 8  |  }
# 
# 
# JAGS data:
# 
# 1  |  "y" <- c(2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 3, 2, 3, 4, 3, 2, 3, 3, 4, 3, 3, 2, 3, 3, 2, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 3, 3, 3, 2, 3, 3, 3, 3, 2, 2, 4, 3, 3, 3, 3, 3, 3, 4, 4, 3, 2, 2, 3, 3, 3, 4, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 2, 3, 4, 4, 3, 2, 3, 2, 3)
# 2  |  "n" <- 100
# 3  |  
# 
# 
# JAGS chains initial values / end states:
# 
# Chain 1:
# 
# 
# 
# Chain 2:

## Model using the dinterval to express discretizing:
breaks <- 1:5
model2 <- "model {
    for (i in 1:n) {
           y[i] ~ dinterval(true[i], breaks)
           true[i] ~ dnorm(mu, tau)
       }

    mu ~ dnorm(3, pow(1, -2))
    tau ~ dgamma(0.001, 0.001)
    }"
data <- list(breaks=breaks, y=data_rounded, n=length(data_rounded))
j2 <- run.jags(model=model2, monitor=params, data=data); j2
# Compiling rjags model...
# Calling the simulation using the rjags method...
# Note: the model did not require adaptation
# Burning in the model for 4000 iterations...
#   |**************************************************| 100%
# Running the model for 10000 iterations...
#   |**************************************************| 100%
# Simulation complete
# Calculating summary statistics...
# Calculating the Gelman-Rubin statistic for 2 variables....
# Finished running the simulation
# Warning messages:
# 1: No initial value blocks found and n.chains not specified: 2 chains were used 
# 2: No initial values were provided - JAGS will use the same initial values for all chains 
# 
# JAGS model summary statistics from 20000 samples (chains = 2; adapt+burnin = 5000):
#                                                                                            
#     Lower95 Median Upper95   Mean       SD   Mode      MCerr MC%ofSD SSeff     AC.10   psrf
# mu   3.2999 3.4059  3.5188 3.4063 0.055658 3.4078 0.00050657     0.9 12072 0.0064681 1.0001
# tau  2.9369 4.3424  5.9953 4.3968  0.79776 4.2138  0.0085467     1.1  8713 0.0046073      1
# 
# Total time taken: 1.2 seconds

[Martyn Plummer]

Observable functions like dround impose strong a posteriori constraints on the parameters. JAGS requires you to start the chain in a state where those constraints are already satisfied. One way to do this is to set true to be equal to the observed value of y at initialization,

inits <- list(true=y)
j2 <- run.jags(model=model2, monitor=params, data=data, inits=inits)

A more refined approach is to set true to be equal to y plus a small random perturbation eps so that true + eps still rounds to to y to the desired accuracy.

inits <- list(true=y+runif(length(y),-0.49,0.49))

Another remark is that your data are rounded to 0 digits, so in your model you should have

y[i] ~ dround(true[i], 0)

[Gwern Branwen]

For reference, here's a simple dround example on my system:

set.seed(2016-11-24)

n <- 10000
Y <- rnorm(n)
X <- round(Y+rnorm(n), digits=1)

model <- 'model {
 for (i in 1:n) {
     roundX[i] ~ dround(X[i], 1)
     mean[i] <- mu + beta*roundX[i]
     Y[i] ~ dnorm(mean[i], tau)
     }
 sd   ~ dgamma(0.01, 0.01)
 tau <- 1/sqrt(sd)

 mu ~ dnorm(0, 100)
 beta ~ dnorm(0, 100)
 }'

library(runjags)
model <- run.jags(model, data = list(n=n, X=X, Y=Y), monitor=c("beta")); model
# ...JAGS model summary statistics from 20000 samples (chains = 2; adapt+burnin = 5000):
#                                                                                                     
#      Lower95  Median Upper95    Mean        SD    Mode       MCerr MC%ofSD SSeff       AC.10    psrf
# beta 0.49143 0.50113  0.5112 0.50113 0.0050271 0.50126 0.000035547     0.7 20000 -0.00032104 0.99997
summary(lm(Y ~ X))
# ...Residuals:
#         Min          1Q      Median          3Q         Max 
# -3.01522501 -0.48632717  0.00091433  0.48105198  2.95046643 
# 
# Coefficients:
#                Estimate  Std. Error  t value Pr(>|t|)
# (Intercept) 0.002548507 0.007142824  0.35679  0.72125
# X           0.502410476 0.005069677 99.10108  < 2e-16
# 
# Residual standard error: 0.7141207 on 9998 degrees of freedom
# Multiple R-squared:  0.4955352,   Adjusted R-squared:  0.4954848 
# F-statistic: 9821.025 on 1 and 9998 DF,  p-value: < 2.2204e-16

[William Beatty]

I think I have run into a similar issue with ~dround. In this example, I have binomially distributed data y[i]. We have an estimate for binomial distribution parameter p = 0.8, and we consider this as data in the model. We would like to estimate n[i] , and we assume that n[i] is distributed with a Pareto distribution with shape parameter alpha = alpha to be estimated and scale parameter c = 1 as data. The Pareto distribution has support [c, ∞). Thus, this seems like an appropriate place to use ~dround. I placed constraints on the initial values for n[i] and tmp[i]. However, I still get the error “Failed check for discrete-valued parameters in distribution dbin” with or without providing initial values as shown in the script below. An example in R is below:``

library(EnvStats)
library(runjags)

## Simulate data ###
nsamps <- 1000
n <- round(rpareto(nsamps, 1, 1.25))  # Simulate the true group size with scale = 1 and alpha = 1.25
p <- 0.8                           # Success probability 
y <- rbinom(nsamps, n, p)         # Simulate observed data y

# p provided as data 
# y[i] represents observed number of animals for group i.
# n[i] represents the true number of animals for group i.
# We want to estimate alpha from the Pareto distribution
cat( "model {
     alpha ~ dgamma(0.01, 0.01) T(1.00001, ) # Mean for Pareto distribution   for alpha < 1 is inf so use trunc
     for (i in 1:nsamps) {
        tmp[i] ~ dpar(alpha, 1)
        n[i] ~ dround(tmp[i], 0)
        y[i] ~ dbin(p, n[i])
     }
     }", fill = TRUE, file = "ParetoSim.jags")
jags.data <- list(y = y, nsamps = nsamps, p = p)
parameters.to.save <- c("alpha")
# Set intial values to meet constraints of dround
inits <- list(n = y,
              tmp = y)
# Also tried this:
# inits <- list(n = round(y/0.8),
#               tmp = round(y/0.8))
model <- run.jags(model = "ParetoSim.jags", monitor = parameters.to.save,
                  data = jags.data, inits = inits, n.chains = 1)