same covariate on persistence and colonisation probabilities in dynamic models

[Marco Salvatori]

Dear all,

I am using dynamic occupancy and N-mixture models and I was wondering what is the effect of using the same set of covariates on both persistence and colonisation probabilities (or survival and recruitment in the case of N-mixture)? Is there any evidence that they can balance each other out resulting in a null net effect? Does anybody know which are the guidelines in such cases: should one avoid to use the same covariate on both probabilites or is it statistically reasonable if one doesn't have specific hypotheses regarding the effect of covariates on only one of the two probabilities (persistence and colonisation) ?

Thanks,

Marco

[Quresh Latif]

I am unaware of papers that formally test for the problem you describe, although I think there are plenty of papers that put the same covariate on both. I think it's worth thinking about what's going on under the hood to decide if there is a potential problem.

For occupancy, colonization estimates the probability that a previously unoccupied site becomes occupied, so the data informing this parameter are sites that are occupied in t-1 (or at least estimated to be unoccupied). Persistence is the probability that a previously occupied sites stay occupied, so the relevant data are sites that are occupied in t-1. The two parameters are therefore fundamentally complementary in terms of the data that inform them, so it seems unlikely to me that putting the same covariate on both would be problematic. In fact, without strong hypotheses about mechanisms, my inclination would be to put all covariates of interest for dynamics on both colonization and persistence.

For abundance, things might be a bit more muddled since any change in abundance from t-1 to t could reflect an infinite number of combinations of survival and recruitment events. My inclination would still be to put covariates of interest on both as I imagine declines in (estimated) abundance will tend to inform survival, whereas increases in abundance will tend to inform recruitment. An alternative would be to simply model population growth (i.e., the sum of survival and recruitment), in which case you would just put covariates on the one parameter.

[Matthijs Hollanders]

Hey Marco,

Good question. I went and simulated some data (attached). I did it pretty quickly so there might be an error but it doesn't seem like the covariates get recovered well when they're in opposite directions (that is, positive effect on colonisation and negative effect on emigration). 

One thing I've been working on with these models is a structural equation modeling (SEM) approach. The idea is that you postulate each site has a latent variable associated with it and this latent variable affects the colonisation and emigration probabilities through "factor loadings". Then you can model covariates on just that latent site-level trait, instead of both colonisation and emigration probabilities. It's a dimension-reducing technique that's also described in the new spAbundance paper.

Cheers,

Matt

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/48e7aee5-58d4-4500-84ed-ea7b95a9c6ffn%40googlegroups.com.

[Matthijs Hollanders]

Hey Quresh,

Depends on the parameterisation. If you do colonisation (gamma) and emigration (epsilon) probabilities:

matrix(c(1 - gamma, gamma,

                epsilon, 1 - epsilon),

             2, 2, byrow = T),

The matrix now is compatible with continuous time models too if you construct it as follows where they are now rates instead of probs:

matrix_exp(matrix(c(-gamma, gamma,

                                    epsilon, -epsilon),

                                 2, 2, byrow = T)). 

Now it’s more intuitive that it could actually be problematic putting the same covariates on both parameters. 

In my simulation there was no visual problem with sampling or posterior correlations. It might have just slipped under the radar. Another thing that seemingly slipped under the radar is that multistate models with state-specific transition and recapture probs aren’t identifiable (for example multistate CMR without secondary surveys). 

Matt

Dr. Matthijs Hollanders

Postdoctoral Research Fellow – Faculty of Science & Technology | University of Canberra

Statistical Consultant – Quantecol

--

*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/ad2f426b-997d-4619-b881-7f1cad0a054cn%40googlegroups.com.

[José Infante]

Hi Marco,

Let's see what the experts say, but colonization, persistence, survival, and recruitment are all estimated from different data. For example, in colonization probability, the "occurrences" are cases where an unoccupied site becomes occupied in the next season, and "absences" are those unoccupied sites that remain unoccupied. For persistence probability, the "occurrences" are cases where an occupied site remains occupied in the next season, and so on. Therefore, the estimation comes from separate response variables. If you hypothesize an association between both parameters, there shouldn't be an impediment to include a given covariate in both. You could run some simulations to be sure.

Best,

José 

--

*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/48e7aee5-58d4-4500-84ed-ea7b95a9c6ffn%40googlegroups.com.

[Jim Baldwin]

If it is reasonable to think that the same covariate can affect persistence and colonization probabilities, then that covariate should certainly be included in both parts of the model no matter what the "potential confounding" might be.  (Which is what Quresh Latif already mentioned.)

But I think your question is "How can one detect some sort of issue from using that covariate for estimating the persistence and the colonization probabilities?"  I would say one should (as with any other pair of covariates) look at the estimated parameter correlation matrix where values close to +1 or -1 should cause further investigation.

Jim

--

[Matthijs Hollanders]

I'm not sure about this (Google groups doesn't let you quote text?):

• colonization, persistence, survival, and recruitment are all estimated from different data. For example, in colonization probability, the "occurrences" are cases where an unoccupied site becomes occupied in the next season, and "absences" are those unoccupied sites that remain unoccupied. For persistence probability, the "occurrences" are cases where an occupied site remains occupied in the next season, and so on. Therefore, the estimation comes from separate response variables

This argument seems akin to saying beta is estimated differently depending on the value of y_i in y_i = alpha + beta * x_i. The probability of z is z ~ categorical(P[z[i, j - 1], 1:2]). Depending on the value of z[i, j - 1], either the colonisation or emigration probabilities (or whatever complement you happen to be choosing) are selected. I don't believe it is correct to see these as independent processes, as more fundamentally the probability of z[i, j] is given by the transition probability matrix which includes both parameters. This point is a bit more intuitive when you marginalise out z[i, j] using the forward algorithm, where you're essentially summing the probabilities associated with each possible state of z.

I also disagree with this:

• If it is reasonable to think that the same covariate can affect persistence and colonization probabilities, then that covariate should certainly be included in both parts of the model no matter what the "potential confounding" might be.

It's definitely a problem if the data generating process cannot be recovered in your model. 

I started writing out the SEM parameterisation and realised it was mathematically equivalent to the traditional approach in the case of one predictor. In both cases there's one covariate that influences both gamma and epsilon. Whether you estimate the coefficients directly (traditional approach) or construct a site-level trait completely determined by the covariate and estimate the factor loadings, there's nothing different happening really. However, this changes for multiple predictors. The number of parameters to estimate for the traditional occupancy is X * 2, where X is the number of covariates and 2 is because it's for both colonisation and emigration. With SEM and one latent trait per site, the numbers of parameters to estimate is X + 2, because you estimate one set of coefficients for each predictor in X and 2 factor loadings for the two rates The SEM approach makes a lot of sense with random effects, where instead of estimating site-level random effects for both colonisation and emigration, you just estimate one site-level trait which then maps to the colonisation and emigration rates using the two factor loadings. This random site effects can be then be affected the covariate.

I ended up running some simulation-based calibration (SBC). The idea is you simulate parameters from your priors, generate datasets from these parameters, and then estimate the parameters for each dataset. If everything is working as it should the ranks of the draws with respect to the input parameters should be uniform. I simulated 50 sites over 20 occasions with 3 secondary occasions per primary. I used normal(0, 1) priors for both logit intercepts and coefficients of colonisation and emigration probabilities, and beta(6, 6) prior for detection probability to avoid extremes. First I ran 500 simulations with intercepts-only, so no coefficients. This one already had a problem, it appears that emigration probabilities (epsilon) and recapture probabilities ( p)  are often overestimated (where rank is defined as the number of draws < simulated value) (see attached).

I'll keep looking into it, but I think it's a great question.

[Matthijs Hollanders]

Update: fixed an error in my code and ran SBC with the coefficients. These plots look like, so doesn't look like there's a problem with having the same covariate on both parameters in the ecological process!

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/e46c90ad-98bc-4815-86b5-93d1ce8811cen%40googlegroups.com.

[Daniel Linden]

The potential confounding comes from the fact that group membership for sites (z[i,t] = 0 or z[i,t] = 1) is partially latent when detection probability < 1.  If detection were perfect, then there would be no uncertainty in group membership and the logistic regressions would be estimated with separate groups of sites each time step.  Sample size issues aside, using simulations with p=1 and then p<1 would help illustrate this.

The SEM approach that Matt describes sounds promising.

[Qing Zhao]

Nice discussion.

Here, I ask a different question: what will happen if x only influences persistence (phi) but is wrongly put on both persistence and colonization (gamma)? I simulated data while x only influences phi, and gamma is influenced by the occupancy status of surrounding sites. Using a model that has x on both phi and gamma to analyze the data, based on 3 simulations, it seems the model tends to underestimate the effect of x on persistence.

Results and code are attached.

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/635556dd-6881-449c-a42f-efd0634e48b7n%40googlegroups.com.

[Matthijs Hollanders]

Hey,

I ran another SBC to check where beta[2] (the effect of the covariate on the emigration probability) was always 0 but was estimated in the model anyway. I did not detect any issues as per these rank histograms.

Cheers,

Matt

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/CAE-HdUFtR_pgKSfM4u-T-9CeGHxYLzcnaQ0sRDDdPchC2-9O9g%40mail.gmail.com.

[Quresh Latif]

Qing, can you explain the decay parameter in colonization? It seems like colonization is declining with distance from the source sampling unit in your simulation? Could the inclusion of that process explain why you are finding a problem but Matthis is not?

This seems related to another discussion we are having in our group, namely when is it advisable or valid to include the same covariate on both abundance/occupancy and detectability?

Both of these questions seem to fall into a broader topic - when is it useful/advisable (or not) to put the same covariate on two different sub-models within a hierarchical model. Without having implemented a substantive literature search, I am aware of one paper that looks into this topic (Monroe et al. 2019). Is anyone aware of any others off the top of your head? Does this seem like an unstudied or understudied question?

Quresh S. Latif
Biometrician

Bird Conservancy of the Rockies

230 Cherry St., Ste. 150, Fort Collins, CO 80521

970-482-1707 (ext. 15)

Connecting people, birds and land

www.birdconservancy.org

To view this discussion on the web visit https://link.edgepilot.com/s/90db0c71/2xSGu7mkFUyEBI465hMqGQ?u=https://groups.google.com/d/msgid/hmecology/48e7aee5-58d4-4500-84ed-ea7b95a9c6ffn%2540googlegroups.com.

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.

To view this discussion on the web visit https://link.edgepilot.com/s/ebca326e/Y4ocw0C6eUaMzZBven_LVA?u=https://groups.google.com/d/msgid/hmecology/635556dd-6881-449c-a42f-efd0634e48b7n%2540googlegroups.com.

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.

To view this discussion on the web visit https://link.edgepilot.com/s/c80bd588/Xq0EfQqOfk6E2UMgytn6Eg?u=https://groups.google.com/d/msgid/hmecology/CAE-HdUFtR_pgKSfM4u-T-9CeGHxYLzcnaQ0sRDDdPchC2-9O9g%2540mail.gmail.com.

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---

You received this message because you are subscribed to a topic in the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this topic, visit https://link.edgepilot.com/s/2ef53bd9/eQvwg4TO0kqwrojAGfBtgQ?u=https://groups.google.com/d/topic/hmecology/xnn3bP2tw5Y/unsubscribe.
To unsubscribe from this group and all its topics, send an email to hmecology+...@googlegroups.com.
To view this discussion on the web visit https://link.edgepilot.com/s/e3697d61/pIZevnUILE2F4iKsyoZuSw?u=https://groups.google.com/d/msgid/hmecology/CAKMRE2u4UgcsWgLY7z2GxhoAtmHOE8N2W2RwUz6mt4Xf4uYTbQ%2540mail.gmail.com.

[Matthijs Hollanders]

I don't have the reference off the top of my head but in one of Kery's books (maybe BPA?) I'm pretty sure in one of the first chapters introducing observation models they show the example of having a covariate that affects each process in different directions and it recovered it. The SBC that I've reported on in this thread also suggests that there doesn't seem to be any issues with at least having the same covariate on both parameters of the state process.

[Marc Kery]

Dear all,

yes, one of the big advantages of these 'explicit' hierarchical models is that we can put the same covariate in the state and the detection model. Many years ago I had been wondering about this and checked it out for the static Nmix model, which led to Kery, Auk, 2008. 

To be honest, I haven't been following this discussion about the dynocc model (sorry!). However, I think that in principle one should absolutely be able to put the same covariate into persistence and colonization. Ignoring the rescue effect, these are different processes !

Indeed, in a detection-naive analysis of single-visit data, one could use separate logistic regressions for each: for each annual interval {t, t+1}, we could first subset the data set to the cases where y[t] = 0 and then model y[t+1] in a logistic regression with colonization for the success probability, and then subset to the cases where y[t] =1 and doing the same, which yields estimates of persistence.

However, as we go towards very small sample sizes and parameter values close to 0, I'd expect things to fall apart. As always, there is no substitute to running a couple simple simulations to check your specific case, when you are in doubt.

Best regards  --- Marc

---

From: hmec...@googlegroups.com <hmec...@googlegroups.com> on behalf of Matthijs Hollanders <matthijs....@gmail.com>
Sent: Tuesday, July 2, 2024 09:59
To: Quresh Latif <quresh...@birdconservancy.org>
Cc: Qing Zhao <white...@gmail.com>; hmecology: Hierarchical Modeling in Ecology <hmec...@googlegroups.com>
Subject: Re: same covariate on persistence and colonisation probabilities in dynamic models

 

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/CAKMRE2tzvregHhjxxg%3D2hzFS%3D9N2NOcwRqF5Kb6m-GcHHg%3DFag%40mail.gmail.com.

[Jim Baldwin]

If I've read the responses correctly, no example of a problem with including the same covariate into persistence and colonization has been given.

While running simulations for any model of any complexity is good practice, I suspect that while those responding would typically do that, I wouldn't be surprised if the vast majority of folks using these models to make management decisions don't do that.

So for the situations where there might be some issue (as Marc states that inadequate sample size might be the major culprit), what summary statistic from the output should one use that would suggest there's an issue?  Or is defining potential issues what needs to happen first?

Jim

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/ZR0P278MB08695A80F668ABE0B18AF0BDEBDC2%40ZR0P278MB0869.CHEP278.PROD.OUTLOOK.COM.

[Qing Zhao]

Just based on a few additional tests, it seems there IS a problem (bias) with including the same x on both persistence and colonization when mean detection probability is low (0.2), but the bias seems to disappear when mean detection probability is moderate (0.5). Also, based on some quick search (two papers attached), it seems detection rarely goes above 0.5. I'm afraid people have to learn to run simulations for their specific situations, but overall, just blindly throwing the same covariate(s) on persistence and colonization doesn't seem to be a good idea to me.

Best,

Qing

[Qing Zhao]

And to clarify, this is not a spatial model; the analytical model matches the data-generating model.

[Marc Kery]

Dear Qing, all,

weeeelllll .... perhaps this is a case where there is a difference between "is" and "ought". I think that right here is the plan for a useful paper to be written, in which the conditions are explored under which the same covariate can be put on phi and on gamma, in terms of things such as sample sizes (number of sites, number of years, number of replicate occasions per year) and parameter values. Ideally, one could also obtain some simple diagnostic features of the data (?) that could guide one's model building when one does not want to run a simulation study tailored to one's study.

Best regards  --- Marc

---

From: Qing Zhao <white...@gmail.com>
Sent: Tuesday, July 2, 2024 17:23
To: Jim Baldwin <jbald...@gmail.com>
Cc: Marc Kery <marc...@vogelwarte.ch>; Matthijs Hollanders <matthijs....@gmail.com>; Quresh Latif <quresh...@birdconservancy.org>; hmecology: Hierarchical Modeling in Ecology <hmec...@googlegroups.com>

[Matthijs Hollanders]

I think another decent heuristic is to model site-level rates or probabilities are bivariate normal (log or logit depending). If the correlation is high between them, you might reduce the dimensions with a single latent trait per site with covariates and factor loadings (coefficients) for the two rates/probabilities. If not, you might proceed with independent sub models for each. 

Dr. Matthijs Hollanders

Postdoctoral Research Fellow – Faculty of Science & Technology | University of Canberra

Statistical Consultant – Quantecol

[John Clare]

Hi Marco (et al.),

For a dynamic occupancy model, the effect of x on colonization and persistence processes can be viewed (or specifically parameterized) as main effects and interaction effects depending on the value of a two-level factor (the previous occupancy state z_it-1, which happens to be latent). I don't see any conceptual/statistical problem with model specifications that range from just two main intercepts(phi~1, gamma~1), main effects and an additional additive term (gamma ~ X, phi~1, or an autologistic model),  interaction effects and an additional additive term (gamma ~ 1, phi ~ X), all the above (gamma ~ X, phi ~ X), or specific sub-cases (gamma ~x1, phi ~ x2) . Could  use any of those specifications in some standard regression analysis, and so I think the same applies here with all of the normal caveats.  (Some will fit better, some may pose additional challenges in general or very specific cases for reasons [data/parameter/other issues] previously mentioned).

With regard to balancing out, I think you mean that the partial derivative of psi_*t or psi_i* (itself derived from phi, gamma, etc.) with respect to x1 is 0?  I think this could happen, but also that it would be pretty unlikely. (It's easier to imagine when x1 is a factor; I'm struggling to imagine this when x1 is continuous, although I guess the mean derivative could be 0). 

John

---

From: hmec...@googlegroups.com <hmec...@googlegroups.com> on behalf of Marco Salvatori <marco.s...@unifi.it>
Sent: Thursday, June 27, 2024 4:45 AM
To: hmecology: Hierarchical Modeling in Ecology <hmec...@googlegroups.com>
Subject: same covariate on persistence and colonisation probabilities in dynamic models

 

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/48e7aee5-58d4-4500-84ed-ea7b95a9c6ffn%40googlegroups.com.

[John Clare]

(Whoops, gamma~x, phi~1 is not autologistic at all! general interaction idea is that f(x) could look one way or be a constant if previous z = one thing and look different if previous z = other thing).

On Jul 2, 2024, at 1:26 PM, John Clare <jcl...@wisc.edu> wrote:



To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/MW4PR06MB83153C1423EAC8F402B4F13E8FDC2%40MW4PR06MB8315.namprd06.prod.outlook.com.

[Marc Kery]

Dear all,

 

Here are results from a simple simulation, where I generated 10,000 data sets where a single site covariate affected both phi and gamma. The intercepts of phi and gamma were varying among the years according to a uniform(0.5, 1) and a uniform(0, 0.5), respectively, while p was randomly picked from a uniform(0.01, 0.99) and kept constant within each data set. Then, I fitted the following model in unmarked:

 

collect(psiformula= ~1, gammaformula =  ~ X, epsilonformula = ~ X, pformula = ~ 1, …)

 

That is, almost the data-generating model, but without the year effects in gamma and epsilon. Note also that the data simulation is in terms of persistence, while in the model we have as a parameter 1 minus persistence, i.e., extinction. This will switch the sign of the coefficient of X in this parameter.

 

Here is a graphical overview of the results: the overall sampling distribution of the coefficients of X in phi and eps are these (this is focused on values near 1, and there are some more extreme estimates on either side, sometimes far more extreme):

 

 

 

And here is a focused picture of the % relative bias of these coefficients, as a function of the value of detection probability p:

 

 

My conclusions: these parameters can be estimated OK, but there is a lot of variability when p is less than about 0.2. In addition, we note an ever so slight negative bias in the coefficients for both parameters, which I would guess is attributable to the model mis-specification in terms of the year variation in phi/eps and gamma that is present in the data, but not specified in the data analysis model.

 

R code attached for anyone who’s interested in checking out or perhaps building on this for more extensive sims.

 

Best regards  -- Marc

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/D301C054-1C8E-4687-A47D-A2617F8F61AA%40wisc.edu.