Dmulti distribution producing posterior samples <0

[EN Pete]

Hi I have a hierarchical model that assumes a multinomial likelihood function and hierarchical model multinomial probabilities. I assumed the dmulti() likelihood in nimble would ensure the the multinomial probability posterior samples could not go below 0, but Im getting negative posterior samples. How is that possible? Im not able to share the data due to data sharing resitrctions. 

##################################################################
# Misclassification Data Model
##################################################################
for(j in 1:N_j){
 x_multi_j[j,1:5] ~ dmulti(rho_stb[gets_j[j], gett_j[j], 1:5], tot_j[j])
}


for(s in 1:Nstates){
for(t in 1:Nyears){

rho_stb[s,t,tpindex] <- sens_st[s,t] * (truePM_st[s,t] - omega_st[s,t] * (1-rho_env_st[s,t] ))
rho_stb[s,t,fnindex] <-  (1-sens_st[s,t]) * (truePM_st[s,t] - omega_st[s,t] * (1-rho_env_st[s,t] ))
rho_stb[s,t,upindex] <- truePM_st[s,t] - (rho_stb[s,t,tpindex] + rho_stb[s,t,fnindex])
rho_stb[s,t,fpindex] <-  (1-truePM_st[s,t]) - (1- omega_st[s,t])*(1- rho_env_st[s,t])
rho_stb[s,t,unindex] <- 1- (rho_stb[s,t,tpindex] + rho_stb[s,t,fpindex] + rho_stb[s,t,fnindex] + rho_stb[s,t,upindex] )

# truePM_st[s,t] ~ dunif(omega_st[s,t] * (1- rho_env_st[s,t] ), rho_env_st[s,t]  + (omega_st[s,t] * (1-rho_env_st[s,t] )))

truePM_st[s,t] ~ dunif(0,1)

rho_env_st[s,t] ~ dunif(0, 1)
# omega_st[s,t] ~ dunif(0, 1)

# sens_st[s,t] ~ dunif(0.2,1) phi(eta_st[s,t])
 sens_st[s,t] <-  1/(1 + exp(-eta_stk[s,t,1]))
 omega_st[s,t]<-  1/(1 + exp(-eta_stk[s,t,2]))


 eta_stk[s,t,1] <- mean_eta_sk[s,1] + phi_s[s] + inprod(Dcomb[t,1:(Nyears-1)], delta_smk[s,1:(Nyears-1),1])
 eta_stk[s,t,2] <- mean_eta_sk[s,2] + phi_s[s] + inprod(Dcomb[t,1:(Nyears-1)], delta_smk[s,1:(Nyears-1),2])
}

#Deltas are deviations away from mean
for(tindex in  1:(Nyears-1)){
delta_smk[s,tindex,1] ~ dnorm(0, sd = sigma_delta_se)
delta_smk[s,tindex,2] ~ dnorm(0, sd = sigma_delta_om)
}


mean_eta_sk[s,1] ~ dnorm(eta_us_k[1], sd = sigma_se)
mean_eta_sk[s,2] ~ dnorm(eta_us_k[2], sd = sigma_om)

}

# ICAR
phi_s[1:Nstates] ~ dcar_normal(adj[1:L], weights[1:L], num[1:Nstates], tau = 1, zero_mean = 1) # its scaled so tau = 1 zero_mean = 1 written before.


eta_us_k[1] <- log(sens_us/(1 - sens_us))
eta_us_k[2] <- log(omega_us/(1 - omega_us))
sens_us ~ dunif(0,1)
omega_us ~ dunif(0,1)
sigma_delta_se ~ T(dnorm(0,1), 0, )
sigma_delta_om ~ T(dnorm(0,1), 0, )
sigma_se ~ T(dnorm(0,1), 0, )
sigma_om ~ T(dnorm(0,1), 0, )
rho_se ~ dunif(-1,1)
rho_om ~ dunif(-1,1)
})

[Daniel Turek]

Yes, that should not be happening.  We will look into this ASAP.

I understand you cannot share the data, and that's no problem.  But, can you please indicate to which *which nodes* in your model are data?  Is it some subset of the x_multi_j nodes?  Or otherwise?

Thank you,

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/b6cddc53-006b-4868-a922-e722e0942564n%40googlegroups.com.

[Perry de Valpine]

Sorry for the delay approving these messages. To fight spam recently we had to turn on the moderation feature for new list members. We try to be diligent but sometimes miss something and a few days go by.

On Mon, Mar 25, 2024 at 4:30 PM EN Pete <epeters...@gmail.com> wrote:

Gosh Im so sorry for emailing you so many times Daniel!

But I do think I figured out what may be happening. In the model code, when I let the truePM_st parameter be a draw from a U(0,1) then the multinomial probabilities go negative. If I put a bound on the truePM_st parameter shown below then they stay positive. 

# truePM_st[s,t] ~ dunif(omega_st[s,t] * (1- rho_env_st[s,t] ), rho_env_st[s,t]  + (omega_st[s,t] * (1-rho_env_st[s,t] )))

Again this doesnt make exact sense to me, because in theory they should never be positive and so an unknown parameter (truePM_st) should be automatically bounded by the multinomial probabilities, I would assume. But i think this is where the issue is occuring, in the allowed bounds for truePM_st. 

Promise this is the last annoying email. 

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/5cf8e5a9-4b21-4079-bd5e-ec0602283f97n%40googlegroups.com.

[Daniel Turek]

Thanks for the additional information.

I admit upfront I have not completely processed and understood your model.  That said, following are some statements, which, with varying degrees of confidence, I believe are true.  If any of these are false, please let me know.

- You are providing observed numbers for "x_multi_j", and these numbers are being passed into the model as "data".

- None of the numbers being provided for "x_multi_j" are NA.

- Therefore, no samplers are being assigned to any elements of "x_multi_j".  That is, "x_multi_j" does not appear anywhere in the Samplers section of the printed MCMC configuration object.

- You are providing values for "tot_j" through either the "constants" list, or the "data" list, when constructing the model.

- When using the U(0, 1) prior distribution for truePM_st[s,t], some values in the interval (0,1) can produce negative values of rho_stb[...]

- When using the more constrained prior distribution:

       dunif(omega_st[s,t] * (1- rho_env_st[s,t] ), rho_env_st[s,t]  + (omega_st[s,t] * (1-rho_env_st[s,t] )))

for truePM_st[s,t], the values of rho_stb[...] cannot go negative.

- In the uncertainty plots you sent, the plots shown for quantities such as rhoFN, rhoTP, rhoUN, etc, represent values from the model of rho_stb[s,t,index], for different values of "index".

- Your original concern has been resolved and addressed.

Thanks!

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CABeTmqihoD%2BmpdF4wODf92cfd5X7saYDC5VkPoHc6Z1Cyb5VYQ%40mail.gmail.com.

[EN Pete]

Hi Daniel 

Here are some clarifications requested. See my responses below

THanks so much for any help you can give.

On Monday, March 25, 2024 at 10:04:50 PM UTC-4 Daniel Turek wrote:

Thanks for the additional information.

I admit upfront I have not completely processed and understood your model.  That said, following are some statements, which, with varying degrees of confidence, I believe are true.  If any of these are false, please let me know.

- You are providing observed numbers for "x_multi_j", and these numbers are being passed into the model as "data".

## Correct. these are multinomial counts

 

- None of the numbers being provided for "x_multi_j" are NA.

##Correct none are NA

 

- Therefore, no samplers are being assigned to any elements of "x_multi_j".  That is, "x_multi_j" does not appear anywhere in the Samplers section of the printed MCMC configuration object.

## Correct not anywhere in the MCMC object

 

- You are providing values for "tot_j" through either the "constants" list, or the "data" list, when constructing the model.

## Correct the sum of x_multi_j is equal to the total sum in the multinomial tot_j 

- When using the U(0, 1) prior distribution for truePM_st[s,t], some values in the interval (0,1) can produce negative values of rho_stb[...]

## Yes Im getting negative rho values, but if I have parametrized correctly this should not be the case since the rho's represent the multinomial probabilities and therefore should also be between 0 and 1.

 

- When using the more constrained prior distribution:

       dunif(omega_st[s,t] * (1- rho_env_st[s,t] ), rho_env_st[s,t]  + (omega_st[s,t] * (1-rho_env_st[s,t] )))

for truePM_st[s,t], the values of rho_stb[...] cannot go negative.

## For most cases. I still found a few samples that are negative. 

 

- In the uncertainty plots you sent, the plots shown for quantities such as rhoFN, rhoTP, rhoUN, etc, represent values from the model of rho_stb[s,t,index], for different values of "index".

##Correct 

- Your original concern has been resolved and addressed.

## Not yet because as states above, Im still getting a few cases of negative values for the rho's. 

## Any help from you is so much appreciated since I have to get this resolved.

[Daniel Turek]

Thanks for the response, and helping answer my questions.

I'd be keen to try and investigate this, but I'd truly need an example demonstrating the problem, in order to investigate deeper.  If you're able to send me an example (even privately, off list), then I'll look into it.

Thanks again,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/30ab9ddd-d6f0-4d4f-a881-ce911fd6255dn%40googlegroups.com.