getWAICdetails() producing empty vectors for fields "..._partialMC" and "..._elements"

[Chris Carleton]

Hi All,

I'm trying to get pWAIC_elements and lppd_elements from my mcmc object and/or mcmc output. Unfortunately, the relevant fields in the return from $getWAICdetails() are just empty vectors, "numeric(0)". I've tried with different models and I've tried removing and resinstalling from Cran. I'm using Ubuntu, R version 4.4.1 and Nimble 1.2.1. I tried the code in the help doc (help(waic)), with the same result ultimately but there are other bugs in the example that I didn't have time to troubleshoot. So, below I have a simplified reproducible example along with relevant output in comment lines (also attached). Any ideas would be appreciated.

Sincerely,

Chris

EXAMPLE:

library(nimble)

library(coda)

code <- nimbleCode({

for(n in 1:N){

y[n] ~ dnorm(mean = mu, sd = s)

}

mu ~ dnorm(mean = 0, sd = 1)

s ~ dexp(rate = 1)

})

mu <- 10

s <- 2

N <- 100

y <- rnorm(N, mean = mu, sd = s)

Rmodel <- nimbleModel(code,

constants = list(N = N),

data = list(y = y),

inits = list(mu = 0, s = 1))

conf <- configureMCMC(Rmodel,

enableWAIC = TRUE)

Rmcmc <- buildMCMC(conf)

Cmodel <- compileNimble(Rmodel)

Cmcmc <- compileNimble(Rmcmc,

project = Rmodel)

output <- runMCMC(Cmcmc,

niter = 1000,

WAIC = TRUE,

samplesAsCodaMCMC = TRUE)

output$WAIC

# nimbleList object of type waicNimbleList

# Field "WAIC":

# [1] 430.6459

# Field "lppd":

# [1] -209.3143

# Field "pWAIC":

# [1] 6.0086

Cmcmc$getWAIC()

# [Warning] There are 4 individual pWAIC values that are greater than 0.4.

# This may indicate that the WAIC estimate is unstable (Vehtari et al., 2017),

# at least in cases without grouping of data nodes or multivariate data nodes.

# nimbleList object of type waicNimbleList

# Field "WAIC":

# [1] 430.6459

# Field "lppd":

# [1] -209.3143

# Field "pWAIC":

# [1] 6.0086

Cmcmc$getWAICdetails()

# nimbleList object of type waicDetailsNimbleList

# Field "marginal":

# [1] FALSE

# Field "niterMarginal":

# [1] 0

# Field "thin":

# [1] FALSE

# Field "online":

# [1] TRUE

# Field "nburnin_extra":

# [1] 0

# Field "WAIC_partialMC":

# numeric(0)

# Field "lppd_partialMC":

# numeric(0)

# Field "pWAIC_partialMC":

# numeric(0)

# Field "niterMarginal_partialMC":

# numeric(0)

# Field "WAIC_elements":

# numeric(0)

# Field "lppd_elements":

# numeric(0)

# Field "pWAIC_elements":

# numeric(0)

summary(output$samples)

[Daniel Turek]

Chris, I took a look, and I think everything is functioning correctly.  The first several fields of the WAIC details are in fact scalars (marginal, niterMarginal, thin, online, nburnin_extra), which are indicating details of the WAIC calculation.  For example, marginal = FALSE indicates this is the "conditional" WAIC calculation (not the "marginal" calculation); see help(WAIC) for more information on that.

The next several elements (WAIC_partialMC, lppd_partialMC, pWAIC_partialMC, niterMarginal_partialMC) are quantities specific only for the "marginal" WAIC calculation.  The fact that those are empty vectors (numeric(0)) is consistent with this calculation *not* being the "marginal" WAIC, thus these fields are empty.

The final 3 fields (WAIC_elements, lppd_elements, pWAIC_elements) are indeed relevant to the calculations for online, conditional WAIC, as you're doing.  However, these elements are suppressed by default, since depending on the model, they can be somewhat large.  To provide those final elements, use:

Cmcmc$getWAICdetails(returnElements = TRUE)

Using that, those final fields will be populated with vectors of numbers.  And again, more information about the numbers you'll see is available in help(WAIC), or in the references therein.

I hope this helps.  Keep in touch with other questions.

Daniel

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/a17027a1-9246-4303-afa8-572ea17a33c9n%40googlegroups.com.

[Chris Carleton]

Thanks Daniel that worked. How can I relate the pWAIC values to my data nodes? I want to know which nodes (observations) are triggering the warning about some pWAIC values being > 0.4 in order to help me evaluate my models.

Best,

Chris

[Daniel Turek]

Chris, unless you specify the "dataGroups" element of the (named list) "controlWAIC" argument to configureMCMC (or to buildMCMC), then the data node groups are simply each of the individual data nodes.  So, in the example you gave, there are 100 data nodes (the 100 y[i] nodes), which themselves are 100 "groups" each containing a single data node, and you can find the element(s) of the pWAIC_elements vector (of 100 values) which exceed 0.4.  Those values of pWAIC_elements which exceed 0.4 are triggering the warning, and are in direct correspondence with the 100 data nodes (the y[i]'s).

Let me know if this doesn't make sense.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/1f327655-8f3a-44e1-8846-ff09cfe6a969n%40googlegroups.com.

[Chris Carleton]

Hi Daniel,

That's clear with the simple model, and thanks for the confirmation. I should have provided more context. I'm running another pair of models (see below) and want to compare these. I'm getting warnings about pWAIC > 0.4 for about 5% of the dataset in each case (data is the same for both models, naturally). So, following Vethari et al., I'm going to look into pointwise lpd comparisons and I wanted to be able to identify the data points that have extreme values (both lpd and pWAIC). I've not configured any dataGroups, so I was expecting to see length(lppd_elements) == length(pWAIC) == N where N is the number of data points in my model. But, when I extract the vectors and check the lengths, I'm getting more elements in each vector than N, despite having not passed anything for the dataGroups argument. It happens to be the case that my models are hierarchical (and a bit complex), so I thought maybe this grouping was being done explicitly by Nimble, but length(..._elements) != K, where K is the number of data groups in my models. So, I'm not sure what's going on exactly or how to determine which elements of lppd_elements or pWAIC_elements correpsond to my data nodes (assuming there are extra elements having to do with lifted nodes or other parameters or something). From the docs (help(waic) for instance) I assumed that length(..._elements) would be either N or K.

Here're the two models (hopefully it's obvious to you what's going on):

# power law

scalingCode <- nimbleCode({

# monument scaling params

intercept0 ~ dnorm(mean = 0, sd = 5)

sd0 ~ dexp(rate = 0.5)

scaling0 ~ dnorm(mean = 0, sd = 5)

sd1 ~ dexp(rate = 0.5)

for(k in 1:K){

intercept_raw[k] ~ dnorm(0, 1)

intercept[k] <- intercept0 + intercept_raw[k] * sd0

scaling_raw[k] ~ dnorm(0, 1)

scaling[k] <- scaling0 + scaling_raw[k] * sd1

}

# prior for negbinom size parameter (dispersion)

size ~ dexp(rate = 0.5)

# population--area linking params

b0 ~ dnorm(mean = 0, sd = 10)

b1 ~ dnorm(mean = 0, sd = 10)

sigma_pop ~ dexp(rate = 1)

# main model

for(n in 1:N){

pop_mu[n] <- b0 + b1 * x[n] # pop_mu, x are on log scale

pop[n] ~ dnorm(mean = pop_mu[n], sd = sigma_pop) # pop is on log scale

log_mu[n] <- intercept[v[n]] + scaling[v[n]] * pop[n] # log_mu, is on

# log scale

mu[n] <- exp(log_mu[n]) # mu is now on original, non-log scale

p[n] <- size / (size + mu[n])

y[n] ~ dnegbin(prob = p[n], size = size) # original scale for count data,

# but power-law model for mean

y_hat[n] <- (size * (1 - p[n])) / p[n]

}

})

# linear-log

scalingCode_linear_log <- nimbleCode({

# monument scaling params

intercept0 ~ dnorm(mean = 0, sd = 5)

sd0 ~ dexp(rate = 0.5)

scaling0 ~ dnorm(mean = 0, sd = 5)

sd1 ~ dexp(rate = 0.5)

for(k in 1:K){

intercept_raw[k] ~ dnorm(0, 1)

intercept[k] <- intercept0 + intercept_raw[k] * sd0

scaling_raw[k] ~ dnorm(0, 1)

scaling[k] <- scaling0 + scaling_raw[k] * sd1

}

# prior for negbinom size parameter (dispersion)

size ~ dexp(rate = 0.5)

# population--area linking params

b0 ~ dnorm(mean = 0, sd = 10)

b1 ~ dnorm(mean = 0, sd = 10)

sigma_pop ~ dexp(rate = 1)

# main model

for(n in 1:N){

pop_mu[n] <- b0 + b1 * x[n] # pop_mu, x are on log scale

pop[n] ~ dnorm(mean = pop_mu[n], sd = sigma_pop) # pop is on log scale

mu[n] <- intercept[v[n]] + scaling[v[n]] * pop[n] # mu is log scale

p[n] <- size / (size + mu[n])

y[n] ~ dnegbin(prob = p[n], size = size) # original scale for count data,

# but linear-log model for mean

# (b\c mu is still on log-scale)

y_hat[n] <- (size * (1 - p[n])) / p[n]

}

})

[Chris Carleton]

Sorry for the double post here.

I think I figured out where the extra elements are coming from. I have N=885, but 39 of these don't have measurements for 'x[n]'. So, they are imputed on the basis of the model for that variable, 'pop[n] ~...". When I use length(waic_details$lppd_elements) I get 924, which is 885+39. So, I guess these 39 nodes have additional entries in the elements list because of the structure of my model. If I'm correct, I have new questions:

1. How is this affecting pWAIC and WAIC, if at all?

2. Should I be grouping data somehow for the WAIC calculations to account for this?

3. How can I identify the corresponding elements in lppd_elements and pWAIC_elements? Presumably, the 39 'extra' entries are just at the top or bottom of the list, so without having to comb through the source related to the waic calculations in Nimble, maybe you know off the top of your head how these elements are ordered with respect to the model's nodes list?

Any information or ideas would be greatly appreciated!

Best,

Chris

[Chris Carleton]

Sorry for the double post. I just figured out (I think) where the extra elements are coming from. In my analysis, I have N = 885 observations, but only 39 of these have measurements for the variable 'x' in the model, the rest are NA and imputed using the model for pop[n]. When I extract the lppd_elements vector, I get length(lppd_elements) = 924. 924 - 885 = 39. So, I think the extra 39 elements must correspond to these 39 complete observations. The questions are now, then:

1. How is this doubling-up of those 39 data nodes affecting pWAIC/WAIC, if at all?

2. Is there something I should be doing with the configuration (maybe grouping somehow?) to account for this? And,

3. How can I identify the corresponding elements of lppd_elements for both the N data nodes and this subset of 39 observations?

Any ideas here would be greatly appreciated!

Best,

Chris

On Friday, August 23, 2024 at 11:38:37 AM UTC+2 Chris Carleton wrote:

[Chris Paciorek]

hi Chris,

885 of the elements should come from the 885 `y[n]` observations. There shouldn't be any additional entries because there aren't any more data nodes in the model. 

I am wondering if the issue relates to what you provide in `data` vs. `inits` vs `constants`. 

Can you provide the full code for how you create the model, including data/inits/constants values (these don't have to be your real data values if you don't want, but the NA/non-NA patterns need to be what you are actually using)? Then I should be able to reproduce what is happening.

As a side note in terms of the WAIC issue, but likely important in terms of your modeling, if `x[n]` doesn't have a prior, then nimble will not impute the `x[n]` values that are given as NA as part of the MCMC. The values will just stay at whatever you provide for them via `inits`.

-chris

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/7b86f0d9-d9b9-44ea-85fb-6f9ceaf72cdcn%40googlegroups.com.