MCMCglmm中过离散与Hurdle Poisson模型先验设定合理性咨询

Hey there, let's break down your MCMCglmm prior challenges with those count models—kudos for tackling such a large nested longitudinal dataset, that's no small feat!

First: Is Fixing Both Residual Variances to 1 Appropriate? #

Short answer: No, this isn't ideal, and it is an overly informative prior that can skew your results. Here's why:

• For your overdispersed Poisson model: Standard Poisson assumes variance = mean, but overdispersed data has variance > mean. Fixing the residual variance to 1 forces the model to ignore this overdispersion, leading to underestimated standard errors and unreliable inference.
• For your hurdle Poisson (hupoisson) model: This family has two components: a logit model for zero/non-zero status, and a truncated Poisson model for non-zero counts. The logit component’s theoretical residual variance is π²/3 ≈ 3.29—fixing this to 1 misrepresents the inherent variability in the binary process, which will bias your estimates for the hurdle part of the model.

Alternative Prior & Model Adjustment Strategies #

Let’s walk through actionable fixes to avoid relying on fix=1:

1. Strengthen R-Structure Priors (Instead of Fixing) #

The "Mixed model equations singular" error often stems from overly vague priors (your original nu=0.002 is extremely non-informative). Try these adjusted priors:

For Overdispersed Poisson:

Give the residual variance a weak but informative prior instead of fixing it:

prior.poi <- list(
  R = list(V = diag(1)*1, nu = 1, n=0),  # Weak half-Cauchy prior for overdispersion
  G = list(
    G1=list(V = diag(3)*0.02, nu =4),
    G2=list(V=diag(3)*0.02, nu=4)
  )
)

• nu=1 uses a half-Cauchy distribution (common for variance parameters), which puts most weight on small values but allows for larger overdispersion estimates if the data supports it.
• V=diag(1)*1 sets a reasonable starting scale for the residual variance.

For Hurdle Poisson:

Respect the logit component’s theoretical variance while letting the Poisson part estimate overdispersion:

prior.hp <- list(
  R = list(
    V = diag(c(3.29, 1)),  # First value = logit's theoretical variance, second = Poisson scale
    nu = c(3, 1),           # Tighter prior for logit, weaker for Poisson
    n=0
  ),
  G = list(
    G1=list(V = diag(6)*0.02, nu =7),
    G2=list(V=diag(6)*0.02, nu=7)
  )
)

This gives the logit process a prior aligned with its inherent variability, while still letting the data drive the Poisson overdispersion estimate.

2. Simplify Your Model Gradually #

Your current model is extremely complex: multiple fixed effects, interaction terms, and 3x3 unstructured covariance matrices for both community and county random effects. This can overwhelm the MCMC sampler and lead to singular equations. Try:

• Starting with a basic model: Fixed effects + only community random intercepts (no slopes or county effects). Get this running without fix=1 first.
• Add random slopes (for t, then x, then the x:t interaction) one at a time.
• Finally, layer in the county-level random effects.
  This helps isolate which component is causing the singularity.

3. Tune MCMC Sampling Parameters #

Your current sampler setup (burnin=15000, nitt=65000, thin=50) gives only 1000 posterior samples—too few for a complex model. Try:

model.hp <- MCMCglmm(
  # ... your model formula ...
  burnin = 20000, nitt = 100000, thin = 50,  # Gives 1600 samples
  verbose = TRUE,  # Monitor sampler stability
  pr = T, pl = T
)

More samples give the sampler more time to converge, which can resolve numerical instability causing the singularity error.

4. Try Alternative Families for Count Data #

If hurdle Poisson continues to fight you, consider:

• Negative binomial (nbinom): Directly models overdispersion without needing a separate hurdle, and is often more stable than overdispersed Poisson.
• Zero-inflated Poisson (zipoisson): Another approach to excess zeros that might have more intuitive prior handling in MCMCglmm.

5. Check for Collinearity & Data Issues #

• Compute variance inflation factors (VIF) for your fixed effects. Values >10 indicate severe collinearity (e.g., between t and time-varying dummy variables, or county-level aggregates and community-level variables), which can cause singular equations.
• Verify your nested structure: Do any counties have only 1-2 communities? Small groups can make random effect estimation unstable—you might need to aggregate or exclude those counties temporarily.

内容的提问来源于stack exchange，提问作者Ling