stan-dev / projpred

Projection predictive variable selection
https://mc-stan.org/projpred/
Other
110 stars 26 forks source link

Guide on latent_ll_oscale for other families e.g. zero-inflated negative binomial, hurdle negbinom #494

Closed stanleyrazor closed 4 months ago

stanleyrazor commented 4 months ago

Greetings. Thanks for the excellent package

I have been working with arms + prospered, where my dependent variable is a proportion (with many zeros as well), so the options I have been playing with families such as: zero_inflated_beta, and at times I have also multiplied by 1000, to get prevalences (per 1000), and played around with hurdle/zero inflated negbinom and poisson. The best fits were obtained with hurdle negbinom models.

I have challenges running varsel for the hurdle and zero inflated negbinom models, they return NAs with the following warnings:

vs_neonate <- varsel(m2, nterms_max=35,
                     # cv_method = 'kfold', # loo ?
                     validate_search = F,
                     method = 'forward', # forward
                     refit_prj = F,
                     latent = T,
                     verbose = T,
                     resp = 'status1',
                     seed = 32)

summary(vs_neonate, resp_oscale = F) latent_ll_oscale returned only NAs, so ELPD, MLPD, and GMPD will also be NA as long as resp_oscale = TRUE.

(I have set validate_search, and refit_prj to false for speed)

I have read the documentation for latent projection, however, I am not in a position to know how to write a custom function latent_ll_oscale_nebin or latent_ppd_oscale_nebin

Is there any resources, or pointers you could help me with so I could solve this.

Thanks.

fweber144 commented 4 months ago

Hi Stanley,

does the latent vignette help (in particular the section with the negative binomial distribution example)?

stanleyrazor commented 4 months ago

Hi @fweber144 , Thanks for the quick response.

From the latent vignette, it is helpful, and here is the key part I need to configure for the zero inflated beta/negbinom and any other family:

refm_prec <- as.matrix(refm_fit_nebin)[, "reciprocal_dispersion", drop = FALSE]
latent_ll_oscale_nebin <- function(ilpreds, y_oscale,
                                   wobs = rep(1, length(y_oscale)), cl_ref,
                                   wdraws_ref = rep(1, length(cl_ref))) {
  y_oscale_mat <- matrix(y_oscale, nrow = nrow(ilpreds), ncol = ncol(ilpreds),
                         byrow = TRUE)
  wobs_mat <- matrix(wobs, nrow = nrow(ilpreds), ncol = ncol(ilpreds),
                     byrow = TRUE)
  refm_prec_agg <- cl_agg(refm_prec, cl = cl_ref, wdraws = wdraws_ref)
  ll_unw <- dnbinom(y_oscale_mat, size = refm_prec_agg, mu = ilpreds, log = TRUE)
  return(wobs_mat * ll_unw)
}
latent_ppd_oscale_nebin <- function(ilpreds_resamp, wobs, cl_ref,
                                    wdraws_ref = rep(1, length(cl_ref)),
                                    idxs_prjdraws) {
  refm_prec_agg <- cl_agg(refm_prec, cl = cl_ref, wdraws = wdraws_ref)
  refm_prec_agg_resamp <- refm_prec_agg[idxs_prjdraws, , drop = FALSE]
  ppd <- rnbinom(prod(dim(ilpreds_resamp)), size = refm_prec_agg_resamp,
                 mu = ilpreds_resamp)
  ppd <- matrix(ppd, nrow = nrow(ilpreds_resamp), ncol = ncol(ilpreds_resamp))
  return(ppd)
}
refm_nebin <- get_refmodel(refm_fit_nebin, latent = TRUE,
                           latent_ll_oscale = latent_ll_oscale_nebin,
                           latent_ppd_oscale = latent_ppd_oscale_nebin)

However, I'll have to admit I do not understand what the two functions are doing - as I'm a little inexperienced.

I think I'll need to know how to do a (rzinb) function for random number generating in zero inflated negbinom/beta and any other family I need in the future.

Is that all I need to change: the dnbinom and rnbinom part to a r-dist and d-dist (where dust is any future distribution I will need to use ?

Thanks

fweber144 commented 4 months ago

Is that all I need to change: the dnbinom and rnbinom part to a r-dist and d-dist (where dust is any future distribution I will need to use ?

Essentially yes, that's what you need to do (although, if you don't plan to use posterior_predict(), you could omit the ppd_oscale part). However, you need to handle any additional parameters appropriately. In the negative binomial latent-projection case, the precision parameter (reciprocal dispersion) is not projected onto any specific parameter of the (Gaussian) submodels. That's why the reference model's posterior draws for this precision parameter (refm_prec) are re-used (as they are) when calculating the response-scale log likelihood values (in the ll_oscale case) or drawing from the response-scale predictive distributions (in the ppd_oscale case) based on the latent predictors inferred from the projected posterior (and back-transformed to response scale, yielding ilpreds). In case you haven't read it yet, the documentation for extend_family() might also be helpful.

fweber144 commented 4 months ago

Closing now, as this issue doesn't seem to reflect a deficiency of projpred. Feel free to re-open if you think differently.

stanleyrazor commented 3 months ago

Thanks for the response.