question about multi species hybrid mixture models

[Brendan Alting]

Hey Jeffrey, thanks very much for getting back to me the other day. I had another couple of questions if that's ok. I've been running these models with a few covariates for detection and abundance, and for four species, 3 mammals and one reptile. All of the parameters for the covariates converged with rhat <1.05, but the intercepts for all of the species really struggled to converge, sometimes up to 3 rhat and low ess. I'm not entirely sure why this is happening, or if you have any suggestions. It could potentially be due to the sparsity of data, but i think there are enough detections for it to be OK. I was running with n.batch = 2000 and batch.length = 50, so i think the iterations are enough. 

I expect that the community mean will be quite different to the values for the individual species, so I thought that may have been causing an issue. I've experimented with looser priors, but so far no luck. I also removed the reptile as they will have higher abundance than the other 3, but having the same issue with just the 3 mammals. 

I am mostly interested in the responses to covariates, and not specifically in abundance. do you think it's justifiable to state the issues with the convergence of intercept parameters, and focus on the covariate influences? 

This is an example of the code with priors I've used (but I've tested a lot of combinations). n.batch = 1000 (results similar for 2000). Any advice or help would be very much appreciated, or if you can see if theres anything wrong with my prior specifications. I have a combination of different camera trap types, which have clear influences on the different species detections. I've tried poisson models too, but it provides the same results. 

Thanks in advance, any help would be super appreciated. 

beta_inits <- matrix(c(
    4.0, -0.12, -0.65,
    3.8, 0.48, 0.06,
    3.9, 0.31, 0.66
  ), nrow = 3, byrow = TRUE)

alpha_inits <- matrix(c(
    -4.5, 0.0,  1.0,  0, 0,    
    -4.8, 0.0,  1.2,  0, 0,  
    -2, 0.0, -1.2, 0, 0    
  ), nrow = 3, byrow = TRUE)



n.factors <- 1




lambda <- matrix(0, nrow = n_species, ncol = n.factors)
diag(lambda) <- 1
lambda[lower.tri(lambda)] <- rnorm(sum(lower.tri(lambda)), 0, 0.1)


init_list <- list(
  beta = beta_inits,
  alpha = alpha_inits,
  phi = runif(n.factors, 0.1, 5),
  lambda = lambda,
  w = matrix(rnorm(n.factors * n_sites, 0, 0.1), nrow = n.factors),
  fix = FALSE
)

priors <- list(
  beta.comm.normal = list(
    mean = c(3.5, 0, 0),        
    var  = c(2, 1, 1)),
  alpha.comm.normal = list(
    mean = c(-3.5, 0, 0, 0,0),
    var  = c(2, 2, 2, 2, 2)),
  tau.sq.beta.ig = list(
    shape = rep(1.1, 3),
    scale = rep(3, 3)),
  tau.sq.alpha.ig = list(
    shape = rep(1.5, 5),
    scale = rep(3, 5)),
  kappa.unif = list(1, 10),
  phi.unif = list(rep(0.01, n.factors),              
                  rep(10, n.factors)))

spabund_model_hab_det <- sfMsNMix(
  abund.formula = abund.formula,
  det.formula = det.formula,
  data = final_data,
  n.factors = n.factors,
  n.batch = 1000,    
  batch.length = 50,
  priors = priors,
  inits = init_list,
  NNGP = TRUE,
  n.neighbors = 15,
  cov.model = "exponential",
  family = "NB",
  n.chains = 3,
  n.report = 50,
  verbose = TRUE
).

summary(spabund_model_hab_det) 

Samples per Chain: 50000 Burn-in: 5000 Thinning Rate: 1 Number of Chains: 3 Total Posterior Samples: 135000 Run Time (min): 12.912 ---------------------------------------- Community Level ---------------------------------------- Abundance Means (log scale): Mean SD 2.5% 50% 97.5% Rhat ESS (Intercept) 3.8610 0.8421 2.1582 3.8713 5.4817 1.0995 1323 scaledcampdist 0.1273 0.6209 -1.1254 0.1303 1.3508 1.0005 59051 scaledurbdist 0.0062 0.6381 -1.2657 0.0060 1.2746 1.0013 61280 Abundance Variances (log scale): Mean SD 2.5% 50% 97.5% Rhat ESS (Intercept) 3.0127 3.9569 0.6106 2.0156 11.3469 1.0114 35233 scaledcampdist 2.3431 2.8154 0.5296 1.6207 8.5337 1.0005 122633 scaledurbdist 2.6766 3.5414 0.6014 1.8447 9.5842 1.0016 124290 Detection Means (logit scale): Mean SD 2.5% 50% 97.5% Rhat ESS (Intercept) -5.1494 0.9315 -6.8377 -5.2009 -3.1468 1.0510 2286 Habitat2Forest -0.2073 0.6901 -1.5665 -0.2149 1.1843 1.0002 9574 CameraTypeUnsealedroad 0.8587 0.9577 -1.1453 0.9027 2.6381 1.0086 12227 CameraType4wdtrack 0.6701 0.9933 -1.3972 0.7078 2.5375 1.0013 38681 CameraTypeWalktrail 1.3729 0.9174 -0.5958 1.4290 3.0244 1.0355 4145 Detection Variances (logit scale): Mean SD 2.5% 50% 97.5% Rhat ESS (Intercept) 3.3744 4.2798 0.6626 2.2748 12.8542 1.0403 17836 Habitat2Forest 1.9721 2.0948 0.4942 1.4238 6.7681 1.0011 99762 CameraTypeUnsealedroad 5.0715 6.3933 0.9254 3.5177 18.2892 1.0049 11574 CameraType4wdtrack 6.6314 7.3523 1.3907 4.7475 23.0014 1.0059 18420 CameraTypeWalktrail 3.5346 4.2529 0.6848 2.4176 13.0408 1.0035 11195 ---------------------------------------- Species Level ---------------------------------------- Abundance (log scale): Mean SD 2.5% 50% 97.5% Rhat ESS (Intercept)-dingo 4.6598 0.5122 3.5445 4.7496 5.4586 2.0913 83 (Intercept)-fox 4.1393 0.6375 3.0178 4.0876 5.4424 1.4469 57 (Intercept)-quoll 3.1840 0.8350 1.5981 3.2085 4.6767 1.2133 51 scaledcampdist-dingo -0.1731 0.2138 -0.5902 -0.1744 0.2446 1.0636 939 scaledcampdist-fox 0.5340 0.2352 0.0678 0.5347 0.9923 1.0147 776 scaledcampdist-quoll 0.2559 0.2263 -0.2011 0.2574 0.6959 1.0081 1496 scaledurbdist-dingo -0.7256 0.2072 -1.1351 -0.7261 -0.3213 1.0586 1103 scaledurbdist-fox -0.0570 0.2456 -0.5289 -0.0591 0.4303 1.0100 869 scaledurbdist-quoll 0.8315 0.2603 0.3352 0.8250 1.3515 1.0603 1190 Detection (logit scale): Mean SD 2.5% 50% 97.5% Rhat ESS (Intercept)-dingo -6.3620 0.5525 -7.2638 -6.4268 -5.2351 1.9285 52 (Intercept)-fox -6.7480 0.7182 -8.1419 -6.7303 -5.3960 1.4901 89 (Intercept)-quoll -4.5452 0.8315 -5.9567 -4.5770 -2.8321 1.3274 40 Habitat2Forest-dingo -0.4754 0.3775 -1.2050 -0.4767 0.2907 1.0794 159 Habitat2Forest-fox -0.3785 0.4283 -1.2469 -0.3655 0.4196 1.0749 280 Habitat2Forest-quoll 0.0526 0.4495 -0.9339 0.0793 0.8732 1.0162 240 CameraTypeUnsealedroad-dingo 2.5900 0.6228 1.3039 2.6311 3.7234 1.0477 97 CameraTypeUnsealedroad-fox 2.7048 0.6668 1.2837 2.7289 3.9570 1.0202 224 CameraTypeUnsealedroad-quoll -1.0965 1.0548 -3.2029 -1.0942 0.9555 1.0672 743 CameraType4wdtrack-dingo 3.0899 0.4002 2.3241 3.0946 3.8939 1.0524 151 CameraType4wdtrack-fox 2.5746 0.5065 1.5912 2.5746 3.5818 1.0265 346 CameraType4wdtrack-quoll -2.0286 0.8967 -3.8716 -2.0081 -0.2940 1.0447 1134 CameraTypeWalktrail-dingo 2.7107 0.6813 1.1309 2.7578 3.9646 1.3770 125 CameraTypeWalktrail-fox 2.8637 0.6262 1.5990 2.8704 4.0420 1.1054 219 CameraTypeWalktrail-quoll 0.4653 0.9241 -1.3139 0.4587 2.2841 1.0752 393 ---------------------------------------- Spatial Covariance ---------------------------------------- Mean SD 2.5% 50% 97.5% Rhat ESS phi-1 5.029 2.8656 0.2983 5.0375 9.7442 1.0002 29802 ---------------------------------------- NB overdispersion ---------------------------------------- Mean SD 2.5% 50% 97.5% Rhat ESS kappa-dingo 3.2368 2.3508 1.0519 2.2327 9.2780 1.0281 464 kappa-fox 2.1079 1.8858 1.0091 1.3144 8.4564 1.1563 329 kappa-quoll 1.6823 1.4304 1.0059 1.1921 6.9910 1.0155 506

[Brendan Alting]

Wow that is not clear at all! Apologies, screenshots of model output: 

[Jeffrey Doser]

Hi Brendan,

N-mixtures can be notoriously difficult to get to converge, particularly when used with something like camera traps where detection rates are generally quite low. In your case, there could be a few things going on. Here are some suggestions:

• Fitting a multi-species model with just four species is difficult. These models estimate species-specific parameters as random effects, which means that you are trying to estimate random effects with four levels. While in theory this is possible, it makes estimating the variance of the random effects difficult. You might try to fit a single species model to some of your species and see if you can get that to converge.
  
• I can't tell the number of spatial locations you have, but based on the amount of run time it takes for the model and the amount of uncertainty in the spatial decay parameter that it is not extremely high (e.g., less than 50?). If so, that can also contribute to making these models hard to estimate, so it could be worth trying to estimate a non-spatial model as well. However, I don't think this is the main cause of the issue here.
  
• You are only specifying a burn-in of 5000 samples out of the total 50,000 samples that you fit, which may not be enough. You should also take a look at the traceplots to assess what the MCMC chains look like. This might make it clear that you need a larger burn-in.

Hope that helps!

Jeff

--
You received this message because you are subscribed to the Google Groups "spOccupancy and spAbundance users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to spocc-spabund-u...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/spocc-spabund-users/bde12d96-bad4-4742-899d-9703bcd71c91n%40googlegroups.com.

--

Jeffrey W. Doser, Ph.D.

Assistant Professor

Department of Forestry and Environmental Resources

North Carolina State University

Statistical Ecology and Forest Science Lab

Pronouns: he/him/his

[Brendan Alting]

Hey Jeff,

Thanks very much for the advice, I'll try to run the species independently and see if that helps. You're right, there's about 61 cameras we had deployed. The traceplots showed it wasn't a burn-in issue, it was more that they continually bounced around. 

Thanks again, 

Brendan