[R-sig-ME] glmmTMB: how to calculate posterior prob. of structural zero?
Ben Bolker
bbolker @ending from gm@il@com
Wed Oct 3 21:01:09 CEST 2018
If you're fitting a zero-inflated Poisson model, it would be
especially easy because the zero probability is simply exp(-lambda), so
I believe the "posterior"(ish) probability that the zero is due to the
structural rather than conditional part of the model would be
P(zero|structural)/P(zero) =
P(zero|structural)/(P(zero_structural)+P(zero_conditional))
or
prob_struc <- function(model) {
strucprob <- predict(model, type="zprob")
condprob <- exp(-predict(model, type="conditional"))
return(strucprob/(strucprob+condprob))
}
For a negative binomial ("nbinom2", or quadratic parameterization) the
zero probability is (k/(k+mu))^k (I think ... e.g.
dnbinom(0,mu=0.5,size=0.5) is sqrt(0.5)). sigma(model) returns the
overdispersion parameter k, so a function as above but with
k <- sigma(model)
condprob <- (k/(k+predict(model,type="conditional")))^k
inserted should do the trick. (Please test these yourself and make
sure they are sensible before proceeding!)
On 2018-10-03 02:18 PM, Aaron Mackey wrote:
> I'm happily using glmmTMB to fit zero-inflated count models with my data,
> but I'd like to also know which zeroes in my data are more likely (or not)
> to be structural vs. expected from the conditional distribution. I know how
> to use Bayes formula to calculate the posterior, and predict(zinb,
> type="zprob") gives me the prior probabilites for each data point being
> structural or not (respecting the zero inflation part of the model), and
> the likelihoods for the structural components are 1 (if the data point is a
> zero) or 0 (if the data point is not a zero) -- but is there a way to
> extract the likelihood for each zero data point with respect to the
> conditional part of the model?
>
> thanks,
> -Aaron
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list