- Operating System: Windows 10, Intel Xeon core
- brms Version: 2.10.0
Hello: this is my first time modeling with Stan and brms. This forum - the posts and answers within, especially @paul.buerkner - have been so helpful and are the reason I’ve made it this far. So thank you in advance!
I am attempting a zero-inflated negative binomial mixed-effect model, and I am also interested in predicting the shape (dispersion) parameter (so as to assess whether variance is different across years in my model).
I am predicting a count response: cars95, annual reproductive success
Predictors include:
Main effects:
age (discrete)
year (discrete)
mort (0/1)
bwev (0/1)
prev sire (0/1)
age x year interaction
Group-level effects:
age|id (where id is categorical for individual across years)
I am also interested in how the variance of my response (cars95) might vary across years, so thought to also model the dispersion parameter associated with a negative binomial distribution: shape ~ year
Thus, the model I would like to fit is:
model {
cars95 ~ age*year + mort + bwev + prevsire + (age|id), shape ~ year
However, I cannot set a prior on shape to constrain it as I can in a model where I do not predict shape:
model {
cars95 ~ age*year + mort + bwev + prevsire + (age|id)
When I try to constrain shape in the model where I predict shape, I am thrown the error:
Error: The following priors do not correspond to any model parameter:
shape ~ gamma(1,1)
When I set no priors and accept defaults, this model fails.
The priors on a well converged initial model where shape is not predicted include:
prior class coef group resp dpar nlpar bound
1 b
2 b age
3 b age:year
4 b bwev1
5 b mort1
6 b prevsire1
7 b year
8 student_t(3, -2, 10) Intercept
9 lkj_corr_cholesky(1) L
10 L id
11 student_t(3, 0, 10) sd
12 sd id
13 sd age id
14 sd Intercept id
15 gamma(1,1) shape
16 beta(1, 1) zi
The default priors on the non-converged model where shape is predicted include:
prior class coef group resp dpar nlpar bound
1 b
2 b age
3 b age:year
4 b bwev1
5 b mort1
6 b prevsire1
7 b year
8 student_t(3, -2, 10) Intercept
9 student_t(3, 0, 10) Intercept shape
10 lkj_corr_cholesky(1) L
11 L id
12 student_t(3, 0, 10) sd
13 student_t(3, 0, 10) sd shape
14 sd id
15 sd age id
16 sd Intercept id
17 sd year shape
18 sd Intercept year shape
19 beta(1, 1) zi
As I mentioned above, the initial model where I do not predict shape converges well:
Family: zero_inflated_negbinomial
Links: mu = log; shape = identity; zi = identity
Formula: cars95 ~ age * year + mort + bwev + prevsire + (age | id)
Data: c95data (Number of observations: 220)
Samples: 4 chains, each with iter = 15000; warmup = 5000; thin = 2;
total post-warmup samples = 20000
Group-Level Effects:
~id (Number of levels: 60)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 1.82 0.97 0.15 3.86 1.00 4798 9806
sd(age) 0.54 0.22 0.14 1.02 1.00 4365 6019
cor(Intercept,age) -0.63 0.37 -0.96 0.52 1.00 6200 6898
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -4.04 1.77 -7.93 -0.98 1.00 7970 13616
age 1.44 0.61 0.33 2.72 1.00 7962 14286
year 0.10 0.23 -0.35 0.58 1.00 14321 16159
mort1 0.58 0.90 -1.21 2.36 1.00 16313 17632
bwev1 -0.97 0.89 -2.81 0.70 1.00 16957 16902
prevsire1 0.07 0.50 -0.91 1.05 1.00 8872 14221
age:year -0.15 0.08 -0.30 -0.00 1.00 8792 14912
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 1.68 0.95 0.54 4.14 1.00 9684 14440
zi 0.12 0.09 0.01 0.34 1.00 14419 16585
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
I’ve conducted loo tests and a number of post predictive checks on this initial model. Many, but not all, check out. Problems that arise with the initial model are a higher p_loo count than the number of parameters in my model, and a somehwat poor-performing ppc_dens_overlay plot:
Loo output:
Computed from 20000 by 220 log-likelihood matrix
Estimate SE
elpd_loo -172.6 18.5
p_loo 34.2 6.5
looic 345.1 37.1
------
Monte Carlo SE of elpd_loo is 0.4.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 197 89.5% 13
(0.5, 0.7] (ok) 23 10.5% 121
(0.7, 1] (bad) 0 0.0% <NA>
(1, Inf) (very bad) 0 0.0% <NA>
All Pareto k estimates are ok (k < 0.7).
See help('pareto-k-diagnostic') for details.
Post-predictive checks:
PPC - bars:
PPC - ecdf overlay:
PPC- stat, mean
PPC - stat, mean and SD
PPC - ppc_check()
I have multiple questions.
- Why can I not assign a prior on a shape parameter in the model where I try to predict shape?
- What priors should I assign in the model where I try to predict shape by year to improve convergence? – and how should I specify these? (i.e., using ‘class=’ or ‘dpar=’ etc.) Note, I’m assuming I need stronger priors on the family specific parameters, but I don’t know what or how (i.e., the parameters where dpar=shape)
- In my initial model where I don’t predict shape, why is the p_loo from the loo output so high when the model appears to fit decently by the ppc’s?
- You’ll notice that my empirical response data does not exceed 8; is it problematic that the zero-inflated negative binomial distribution allows (occasional) responses that are unrealistically high for empirical data? If so, suggested alternatives or work-arounds?
I’ve included a reproducible example together with a link to data:
library(dplyr)
library(brms)
library(ggeffects)
library(GGally)
library(rstantools)
library(rstanarm)
library(loo)
library(bayesplot)
c95data <- read.csv("ARS95KAc.csv")
c95data$id <- as.factor(c95data$id)
c95data$year <- as.integer(c95data$year) #year (2-8) as integer NOT categorical
c95data$mort <- as.factor(c95data$mort)
c95data$prevsire <- as.factor(c95data$prevsire)
c95data$bwev <- as.factor(c95data$bwev)
#MODELS
mod1 <- brms(cars95 ~ age*year + mort + bwev + prevsire + (age|id) ,
data=c95data,
family=zero_inflated_negbinomial(),
prior= set_prior("gamma(1,1)", class="shape"),
iter=20000, warmup=5000, thin=2, chains=4,
control = list(adapt_delta = 0.99),
cores=6)
#THROWS ERROR ABOUT PRIOR
mod2 <- brms(bf(cars95 ~ age*year + mort + bwev + prevsire + (age|id) ,
shape ~ year),
data=c95data,
family=zero_inflated_negbinomial(),
prior= set_prior("gamma(1,1)", class="shape"),
iter=20000, warmup=5000, thin=2, chains=4,
control = list(adapt_delta = 0.99),
cores=6)
##DOES NOT CONVERGE
mod2_nopri <- brms(bf(cars95 ~ age*year + mort + bwev + prevsire + (age|id) ,
shape ~ year),
data=c95data,
family=zero_inflated_negbinomial(),
iter=20000, warmup=5000, thin=2, chains=4,
control = list(adapt_delta = 0.99),
cores=6)
#LOO
loo1 <- loo(mod1, reloo=T, cores=6)
##PPCs
#post predictive checks
yrep <- posterior_predict(mod1)
y <- c95data$cars95
bayesplot::ppc_stat(y, yrep[1:50, ], stat="mean")
brms::pp_check(mod1)
brms::pp_check(mod1, type="stat_2d", stat=c("mean", "sd"), nsamples=30)
brms::pp_check(mod1, type="ecdf_overlay", discrete=T, nsamples=30)
brms::pp_check(mod1, type="bars", nsamples=30)
DATA: https://drive.google.com/file/d/1T51Cp0m8qHtsf2fZcpqchMQJo6dIdQn1/view?usp=sharing
Thank you (again) in advance!
El