Polyagamma negbin

[paul-vdb]

This is in partly response to: #1536 to solve the strict error message and not the way the user defined the target prior distribution. It also does not include multinomial yet as we discussed that might make more sense as a different sampler.

Some of these changes are discussed here: https://github.com/perrydv/nimbleCoreTeam/issues/568 in the context of reading about the multinomial sampler.

1. To solve the basic problem from the google groups. I've made sure that inflationNodes is initialized even when check = FALSE.
2. Updated so that we can now have a negative binomial.
3. Logic shift in calculating size and observations in code.

The shift in code for actually running occurs in the two functions setDataParams and updateDataParams. The theory for the negative binomial changes what the input data is. For this reason, I've set up some logic to set the data and size the first time, and then to update on each iteration just the stochastic values. Nothing else changes after updating these functions that replaces setSizeParam function. See some of the basic theory here: https://github.com/perrydv/nimbleCoreTeam/issues/568#issuecomment-2744576713

I've added conjCheckAll with default FALSE and worked through the logic of checking the sampler. When set to true, it loops through each observation and does a conjugacy check (when check = TRUE). In the case when check = FALSE, it still does check if size is stochastic for all observations.

Tests have been updated to check the negative binomial and a mixed example of dnegbin and dbinom.

[paul-vdb]

@paciorek I was just trying to use the polyagamma sampler with coefficients that switch on and off, not quite a RJMCMC case.

logit(p[i]) <- beta[1] + beta[2]*age[i]*z + beta[5]*len[i]*(1-z) + beta[3]*(len[i] - mu[age[i], sex[i]]) + beta[4]*(sex[i]-1)
where z ~ dbern(0.5)

This is not allowed but I think is a use case we need to account for. Is there a linkage with RJMCMC? Flagging here as you review some of the logic. Anything that is a constant in the linear relationship is valid, we just need to identify it for when we sample the polyagamma random variable.

[paul-vdb]

@paciorek I found another error in polyagamma.
If a designMatrix is passed by the user, it is overwritten via: defined below the check for control$designMatrix
initializeX <- TRUE

Additionally, I've been trying to use this for an actual project and noticed that if the design matrix is not normalized, which results in the design matrix not building correctly.
Then set beta = c(0,1,0,0,0)
if X[i,2] <- 41
then p[i] <- expit(41) leads to logit(model$getParam(yNodes[i], "prob")) = Inf
Given the current setting of initializeX <- TRUE I am unable to get around this by passing my own design matrix, except by normalizing all predictors.

I've updated the easy fix and maybe just mention in documentation that if a design matrix isn't provided then the user should use standardized predictors?

[paciorek]

@paul-vdb I think this might be stale in light of PR #1569.

But this has negative binomial stuff that we will want in the future, right? (But my thought is that we might not be including that for 1.4.0.)

So should we close this without merging or leave it open for later?

[paul-vdb]

@paciorek Yeah I think this is stale and we might want to hold off here. There is more testing to do and it doesn't seem like anyone is chomping on the bit for this. Given I've written code that allows Poisson too let's not rush this one. Will close.