Uncertainty in fixed parameter values between 0-1

[mcmelnychuk]

Hello all,

I'm wondering if it's possible to incorporate uncertainty into a fixed parameter value. I’m especially wondering about the final, joint “beta” parameter in a CJS model where Phi to the last occasion and p at the last occasion are confounded. If the final p is fixed at some value 0<p<1 with zero variance (to estimate the final Phi, accepting the risk of bias in the final Phi from an incorrect fixed value of the final p), then along with the untangled final Phi, the SE of the final Phi is estimable and returned.

If uncertainty were admitted for that final p, however, I would expect the uncertainty of the final Phi to increase as well. I'm just wondering if it's possible to somehow compensate for this underestimation of uncertainty in the final Phi when the final p is fixed without error.

Thank you for your thoughts,
Mike

[dhewitt]

Mike,

This is curious. Can you explain more about your motivations for fixing the final p and its variance? I.e., what external information are you relying on?

[sbonner]

Hi Mike,

As you point out, the final capture and survival parameters in the CJS are confounded so it's not possible to estimate one or the other separately. What is estimable is their product, phi[k-1]*p[k] (where k is the number of capture occasions). Fixing p on the final occasion is just a trick to reduce the number of parameters by 1 so that all of the remaining parameters are estimable.

But, don't be fooled into thinking that by fixing the capture probability you can suddenly estimate the final survival parameter. By fixing p[k]=1 (the usual choice) the final survival probability output by MARK is really an estimate of the product p[k]*phi[k-1], and the standard error is an estimate of the variability in this product. If you set p[k] equal to another value, say c, then the final survival probability is an estimate of c*p[k]*phi[k-1] and the reported standard error is an estimate of the variability in this product.

In short -- you can't interpret the estimate of phi[k-1] as an estimate of survival. I'm sure that Evan has a better explanation of this in the Gentle Introduction to MARK.

Hope that helps.

Cheers,

Simon

[mcmelnychuk]

Thank you for your replies, Simon and David. I realize that what I’m trying to do is atypical and prone to biases.

This is curious. Can you explain more about your motivations for fixing the final p and its variance? I.e., what external information are you relying on?

I’m using “spatial” CJS models with tagged salmon smolts detected at stations along a migratory route, like the studies in the “blue book” (Burnham et al., 1987). Segments of the migration occur between detection stations, and fish populations typically pass several stations along a migration route, so instead of “time” we have a form like Phi(segment)p(station) for a CJS model. The final few stations are in the ocean, and because of feasibility and funding issues there were few of them. Biologically, I’m interested in making survival inferences to the final detection station even if survival estimates in the final segment, Phi[k-1], are slightly incorrect as a result of incorrect assumptions of a fixed p[k] at the final station. (I’m actually interested in survivorship, the product of segment-specific Phi estimates. Low Phi estimates in the earlier segments generally result in fairly low survivorship leading up to the final segment. Consequently, even if Phi[k-1] is incorrect, the overall survivorship estimates are relatively insensitive to fixed values of p[k].)

Estimates of p at ocean stations before the final station are all fairly similar. Detection probability is determined more by acoustic tags and the geometry of receiver stations than by fish behaviour, and these seem to be fairly constant across several stations (around p=0.9 for one particular tag type and 0.75 for another). Since the final detection station shares similar receiver geometry and environmental conditions as the other ocean lines, I’ve estimated p[k-1, k-2...] in separate (though not quite external) analyses and assumed this value (as fixed) for p[k], with slight adjustments for variation in receiver geometry among stations.

Because of the low variability in p[k-1, k-2...] and the relative insensitivity of survivorship to fixed values of p[k], I’m willing to take the risk of fixing p[k] to estimate Phi[k-1]. I realize that incorrect assumptions of fixed p[k] would lead to incorrect estimates of Phi[k-1]. I suppose I’m trying to estimate Phi[k-1] conditioned on the assumed, fixed value of p[k] at 0<p[k]<1.

By fixing p[k]=1 (the usual choice) the final survival probability output by MARK is really an estimate of the product p[k]*phi[k-1], and the standard error is an estimate of the variability in this product. If you set p[k] equal to another value, say c, then the final survival probability is an estimate of c*p[k]*phi[k-1]

I’ve run two simple CJS models, identical except for in one model, p[k] is fixed at c where 0<c<1, and in the other it is not fixed. In the real estimates generated, the product p[k]*Phi[k-1] is identical between models. So if p[k] is fixed at 0<c<1 then I believe the estimate for Phi[k] becomes what the “beta parameter” product would have been without fixing p, divided by c (rather than multiplied).
I.e., c*Phi[k-1_pFixed] = p[k_pNotFixed]*Phi[k-1_pNotFixed] = B
Can the B parameter be untangled this way and the reported Phi[k-1_pFixed] be considered an estimate of survival, conditioned on the fixed value of p[k]?


Now, to the question of variances. In those separate analyses I mentioned for estimating p[k-1, k-2...], I also estimated their variances with CJS models, and wish to somehow incorporate these into the fixed values of p[k]. Without fixing p[k], the SEs reported in the real estimates for p[k] and Phi[k-1] are non-sensical as we all know. However, if p[k] is fixed, the reported SEs for Phi[k-1] no longer seem unreasonable – they are slightly larger on average than those of Phi[k-2], for example. As Simon indicated, this SE is actually an estimate of the variability of p[k]*Phi[k-1] (divided by?) c. I think another interpretation could be that this is the SE for Phi[k-1] if p[k] were known with zero error. (Am I off track on that?) The reported SEs for the fixed p[k] are, of course, 0 since they are fixed values. I’m wondering if it is possible to fix the SE of p[k] at the same time as fixing p[k] (or else to compensate for what I believe is an underestimation of the SE of Phi[k-1] when fixing p[k] at 0<p[k]<1). I assume that this would increase the SE of the reported Phi[k-1] to better represent the actual uncertainty. My logic in this assumption is that if the SE of p[k]>0, there would be more possible ways or combinations of fitting the product p[k]*Phi[k-1] to the data than in the case where the SE of p[k] is “fixed” at 0 by default.

Any thoughts would be appreciated. Thank you,
Mike

[sbonner]

Hi Mike,

Sorry I didn't respond earlier -- your reply snuck past me.

I’ve run two simple CJS models, identical except for in one model, p[k] is fixed at c where 0<c<1, and in the other it is not fixed. In the real estimates generated, the product p[k]*Phi[k-1] is identical between models. So if p[k] is fixed at 0<c<1 then I believe the estimate for Phi[k] becomes what the “beta parameter” product would have been without fixing p, divided by c (rather than multiplied).
I.e., c*Phi[k-1_pFixed] = p[k_pNotFixed]*Phi[k-1_pNotFixed] = B
Can the B parameter be untangled this way and the reported Phi[k-1_pFixed] be considered an estimate of survival, conditioned on the fixed value of p[k]?

I'm not sure I follow exactly, but I think that this is right. Suppose you set p[k]=c. Then \hat{\phi[k-1]} will actually be an estimate of \phi[k-1]*p[k]/c. If you happened to guess the correct value of p[k] (i.e., p[k]=c for the true population) then you're home free -- p[k] cancels with c and you are actually estimating \phi[k-1]. As you noticed, no matter what value of c you choose to fix p[k] to, the product \hat{\phi[k-1]}/c will be the same, but \hat{\phi[k-1]} will only be an estimate of \phi[k-1] if you happened to hit on the correct value of p[k].

I think another interpretation could be that this is the SE for Phi[k-1] if p[k] were known with zero error. (Am I off track on that?)

Hmmm... again, I think that makes sense. If you are willing to believe that you have chosen the correct value of p[k] then SE(\hat{\phi[k-1]}) will be the correct standard error for \hat{\phi[k-1]}, which in turn will be a proper estimate of \phi[k-1]. If not, then SE(\hat{\phi[k-1]}) will still be the standard error of \hat{\phi[k-1]} -- but it's not clear what value the latter is estimating. I understand your desire to inflate the variance, but I can't see how this makes sense, at least in a frequentist framework....

Since the final detection station shares similar receiver geometry and environmental conditions as the other ocean lines, I’ve estimated p[k-1, k-2...] in separate (though not quite external) analyses and assumed this value (as fixed) for p[k], with slight adjustments for variation in receiver geometry among stations.

It really sounds to me like you want to be a Bayesian here. You have information about p[k] from other receivers with similar structure and you would like to incorporate this information into the estimation of p[k]. The ideal way to handle this is through Bayesian hierarchical modelling. I realize that this may seem a little flippant or even overkill (your hammer does quite work, so try my sonic screwdriver) but this seems to me to be the most sensible way to proceed. A hierarchical model would allow you to put an exact structure to the idea that there is variability in p[k] but that some information is available from other stations. Unfortunately, such models are (currently?) not part of MARK, but a hierarchical CJS can easily be implemented in the BUGS language. You could have a look at the examples in:

Gimenez, O., S. Bonner, R. King, R. A. Parker, S.P. Brooks, L. E. Jamieson, V. Grosbois, B. J. T. Morgan, and L. Thomas (2008). WinBUGS for population ecologists: Bayesian modeling using Markov Chain Monte Carlo methods. Environmental and Ecological Statistics. In press.

which you can find on Olivier Gimenez's website here:
http://www.cefe.cnrs.fr/biom/Permanents ... blications.

Hope this helps.

Cheers,

Simon

[mcmelnychuk]

Thank you, that helps to clear things up.

One approach I’ve used instead of being able to allow the SE of p[k] to be >0 when p[k] is fixed between 0-1 is to also fix p[k] at lower and higher values to provide a bounded range of uncertainty for Phi[k-1]. In the separate analysis where I estimated p[k-1...], I also estimated 95% c.l., and applied these as the low and high fixed values of p[k], which produce high and low estimates of Phi[k-1], respectively. So even if p[k] is fixed without error, there can at least be bounded estimates of Phi[k-1] that hopefully make up for the underestimation of the SE of Phi[k-1]. This seems logical enough to me, although there are so many steps in this estimation of uncertainty that, even if it’s not severely flawed, it’s difficult to explain or justify. Would have liked to KISS.

I agree that a Bayesian approach would make sense with respect to borrowing information for p[k]. The models that I’m mainly dealing with are two large ones, each with 20-30 release groups that have a variable number of detection occasions. Some detection stations are shared among release groups while others are not. Although I consider these groups independently in terms of Phi, I combine them in the same model so groups share p parameters at a given station where it makes sense, and also to maintain certain additive effects on p consistent across groups. This borrowing of information, I think, at least takes the flavour of a hierarchical model (for a different reason than the one you suggested) even though it’s implemented in a frequentist framework. (?) There are other reasons why I would ideally like to specify the likelihoods manually to provide additional flexibility not found in standard models, but for these large models I suspect it would be intractable so haven’t gone far down that road. The Gimenez et al. method and code you referred to (thank you for that) of reading in an m-array and using built-in WinBugs distribution functions seems like an efficient way of coding. Do you (or anyone) have a sense of whether running a CJS-type analysis where p is modelled as a hyperparameter across many release groups would be relatively straightforward to code and computationally feasible?