MCMC random effects model to resolve boundary parameters?

[chris_oosthuizen]

Hi all,

I used MARK to estimate survival probability of rodents caught in live-traps. I have 12 capture occasions and 300 individuals (of which 80 were captured and released at the first occasion).

I started by fitting CJS and time-constant models with trap-dependence effects for p. The model where survival was fully time-dependent was preferred by AIC over the constant model (delta AIC = 3.2) and other covariate models. However, 3 of the 11 survival probabilities of the time dependent model was estimated at the boundary (phi = 1). I was not too happy with this result, and considered two options to improve upon the analysis: model averaging, and random effect models.

I decided to use the MCMC procedures in MARK to fit models with random time effects (1 hyperdistribution - i.e., to estimate mu and sigma for survival). To compare results, I also ran “method of moments” (MoM) random effect models as an alternative approach.

Some of my results are as follows:
# Maximum likelihood (MLE) constant model
Phi (MLE) = 0.69
p (MLE) = 0.32 (caught on previous occasion)
p* (MLE) = 0.16 (not caught on previous occasion)

# MCMC
phi (mu) on the logit scale = -0.06, which translates to a mean survival probability of 0.48 (in R: mu.probability = plogis(mu.logit))
p (mcmc) = 0.32 (caught on previous occasion)
p* (mcmc) = 0.17 (not caught on previous occasion)

# Method of moments random effect
phi (mu) = 0.73.

In the Method of Moments model, 4 of the shrinkage estimates for phi is on the boundary (=1). So this does not appear to solve that issue. None of the annual survival estimates in the MCMC approach are at the boundary.

I have some reservations about these results, where the mean survival rate is similar for the MLE and MoM approaches, but comparatively low for the MCMC analysis (at least according to my calculations).

White et al (2009) suggests that MCMC work well for estimating mu, but that performance may be poor for sigma with ≤ 10 occasions. I have 12 occasions, which is borderline, but I would expect mu to be more similar between the three analyses?

Any thoughts or suggestions will be much appreciated. I can also provide the RMark code and MARK import file (for MCMC) in a private message should anyone be interested to replicate the results (which should not take long to do).

[cooch]

Hi all,

I used MARK to estimate survival probability of rodents caught in live-traps. I have 12 capture occasions and 300 individuals (of which 80 were captured and released at the first occasion).

So, what is it you're interested in - mean survival over the study, or something else? I'm guessing mean, because that seems to be the focus of much of what comes next...

I started by fitting CJS and time-constant models with trap-dependence effects for p. The model where survival was fully time-dependent was preferred by AIC over the constant model (delta AIC = 3.2) and other covariate models. However, 3 of the 11 survival probabilities of the time dependent model was estimated at the boundary (phi = 1). I was not too happy with this result, and considered two options to improve upon the analysis: model averaging, and random effect models.

You should probably consult Appendix F - 'data cloning', to shed some light on the 'boundary estimation' problem. If the estimates are for 'interior' parameters (meaning, not at the beginning or end of a time-series) then you need to suss out if the boundary estimates represent parameters that really are at the boundary, or if its just a data insufficiency problem ('external nonidentifiability').

I decided to use the MCMC procedures in MARK to fit models with random time effects (1 hyperdistribution - i.e., to estimate mu and sigma for survival). To compare results, I also ran “method of moments” (MoM) random effect models as an alternative approach.

<snipped...>

# MCMC
phi (mu) on the logit scale = -0.06, which translates to a mean survival probability of 0.48 (in R: mu.probability = plogis(mu.logit))
p (mcmc) = 0.32 (caught on previous occasion)
p* (mcmc) = 0.17 (not caught on previous occasion)

Whenever anyone runs MCMC, and tries to talk about a 'point estimate' (say, of mu, or sigma, or anything else), the usual questioon I ask to start the chat is 'what moment(s) from the posterior are you looking at?'. Unless the posterior is unimodal, abnd symmetrical, then the mean, median, and mode will diverge -- often dramatically so. I'd be willing to bet (sight unseen) that the -0.06 was based on one moment, while the othe other two (MARK reports 3: mean, median, mode) are quite different. For [0,1] bounded parameters, especially for any parameter nearer 0 or 1 than 0.5, assymetrical posteriors are typical, and the issue of what to base inference on is challenging. Most 'real Bayesians' will say that focus should be on the credible interval, and not the moment from the posterior. Mevin Hooten and coleagues have written a fair but about 'what do you do if you want to talk about point estimates from a Bayesian analysis'. You might seek out his papers.

In the meantime, does the back-transformed HPD (careful how you do it) bound the values from the MoM and 'dot' model approach? If so, then you may be worrying about nothing. Also, what do the median and mode look like? Have you actually plotted the posterior density for mu? Is it multi-modal, or strongly skewed? If either, then why presume that a single moment is even remotely meaningful? This issue is discussed at somne length in Appendix E.

In the Method of Moments model, 4 of the shrinkage estimates for phi is on the boundary (=1). So this does not appear to solve that issue. None of the annual survival estimates in the MCMC approach are at the boundary.

Whether or not it is an issue depends on whether or not the parameters are really 1.

I have some reservations about these results, where the mean survival rate is similar for the MLE and MoM approaches, but comparatively low for the MCMC analysis (at least according to my calculations).

See above.

White et al (2009) suggests that MCMC work well for estimating mu, but that performance may be poor for sigma with ≤ 10 occasions. I have 12 occasions, which is borderline, but I would expect mu to be more similar between the three analyses?

Again, see above. Your MCMC analysis is based largely on you trying to make inference from a single moment from the posterior distribution. In my experience, often a bad idea.

[chris_oosthuizen]

Thanks for the useful reply Evan. It unfortunately took me a while to get a chance to look at this again.

So, what is it you're interested in - mean survival over the study, or something else? I'm guessing mean, because that seems to be the focus of much of what comes next...

Seasonal survival (winter/summer), and survival correlating with environmental (e.g., temperature) or individual covariates (e.g., age and body mass) are all of interest to the principle investigators (who I am assisting with the analysis). However, none of these covariate models were preferred over the fully time dependent or constant models in AIC model selection.

I thought a random effect model is a useful ‘intermediate’ between the constant model and fully time-dependent model as survival are (apparently) time-varying, but because of data insufficiency (see below) some parameters of the fully time-dependent model is poorly estimated. I was also hoping that the ML random effects or MCMC approaches will ‘shrink’ the boundary estimates toward the model-specified mean. Alas, I see this does not always work this way (see Appendix F.6. Limitations & other thoughts and approaches).

Finally, having a mean survival probability is probably useful as it can be used to approximate mean longevity (or is it lifespan?) for the species? (which is unknown).

You should probably consult Appendix F - 'data cloning', to shed some light on the 'boundary estimation' problem.

Thanks for the good recommendation to look at Appendix F. In this case, the 'boundary estimation' problem is (mostly) for 'interior' parameters (parameters 1, 6 and 10). I computed profile likelihood confidence intervals for all the survival parameters. The profile likelihood CI did not change/improve the CI of parameters 1, 6 or 10:

1: Phi = 0.999 (SE = 0.004) Profile CI = (0.00 - 1.00)
6: Phi = 1.00 (SE = 0) Profile CI = (1.00 - 1.00)
10: Phi = 1.00 (SE = 0) Profile CI = (1.00 - 1.00)

Cloning the data from the profile likelihood estimates also did not change the 0 standard error for parameters 1, 6 and 10. Given that the profile likelihood CI are poorly estimated, even with simulated annealing, I conclude these parameters are extrinsically non-identifiable because of sparse data.

Whenever anyone runs MCMC, and tries to talk about a 'point estimate' (say, of mu, or sigma, or anything else), the usual question I ask to start the chat is 'what moment(s) from the posterior are you looking at?'. Unless the posterior is unimodal, and symmetrical, then the mean, median, and mode will diverge -- often dramatically so. I'd be willing to bet (sight unseen) that the -0.06 was based on one moment, while the other other two (MARK reports 3: mean, median, mode) are quite different.

The posterior distribution for mu is unimodel and quite symmetrical. The mean for mu was -0.0602, the median -0.0221, and the mode 0.1671. Back-transformed these correspond to 0.484, 0.494, 0.541, respectively.

For [0,1] bounded parameters, especially for any parameter nearer 0 or 1 than 0.5, assymetrical posteriors are typical, and the issue of what to base inference on is challenging.

Yes indeed: the posterior distributions for the boundary parameters in particular were not unimodal or symmetric. Instead, parameter 1 was multi-modal, and parameters 6 and 10 were strongly skewed. These parameters are clearly difficult to estimate.

Have you actually plotted the posterior density for mu?

I plotted the posterior density for mu on the logit scale - it appeared to be unimodel and symmetrical. At this stage I am not sure how to get a posterior density plot for mu on the real scale?

In the meantime, does the back-transformed HPD (careful how you do it) bound the values from the MoM and 'dot' model approach?

The the 95% HPD for mu on the logit scale is [-1.809; 1.265]. Back transformed, this corresponds to [0.14 - 0.78] if I am correct? (I simply back-transformed the two HPD from the logit scale to the real scale). In this case the HPD for mu will (just) bound the values from the MoM and 'dot' model approach.

back-transformed HPD (careful how you do it)

I think this warning refers to sigma only - or should I also be careful how I back-transform the mu HPD?
(the reference to the Book for back transformation is: Section E.1.2. Estimating the hyperparameters μ and σ should this help anyone else)

[cooch]

Seasonal survival (winter/summer), and survival correlating with environmental (e.g., temperature) or individual covariates (e.g., age and body mass) are all of interest to the principle investigators (who I am assisting with the analysis). However, none of these covariate models were preferred over the fully time dependent or constant models in AIC model selection.

Which should already raise red flags about the quality of the data -- I'm guessing low encounter probabilities?

I thought a random effect model is a useful ‘intermediate’ between the constant model and fully time-dependent model as survival are (apparently) time-varying,

Correct -- it is, in effect, intermediate (and the shrinkage estimates are basically ~equivalent to what model averaged estimates would look like over a saturated set of models between dot -> time).

but because of data insufficiency (see below) some parameters of the fully time-dependent model is poorly estimated.

As suspected (above).

I was also hoping that the ML random effects or MCMC approaches will ‘shrink’ the boundary estimates toward the model-specified mean. Alas, I see this does not always work this way (see Appendix F.6. Limitations & other thoughts and approaches).

Also correct (based on experience to date...there may be particular circumstances where it 'works' better, but I've not had such luck).

Finally, having a mean survival probability is probably useful as it can be used to approximate mean longevity (or is it lifespan?) for the species? (which is unknown).

Well, maybe. If you assume survival is a constant, and doesn't vary with age (both pretty tenuous assumptions), then there ways to calculate 'mean time to event' (being the death event), but given the assumptions, I personally wouldn't hold much faith in this approach.

<following suggestion to look at 'data cloning' approaches...appendix F>

In this case, the 'boundary estimation' problem is (mostly) for 'interior' parameters (parameters 1, 6 and 10). I computed profile likelihood confidence intervals for all the survival parameters. The profile likelihood CI did not change/improve the CI of parameters 1, 6 or 10:

1: Phi = 0.999 (SE = 0.004) Profile CI = (0.00 - 1.00)
6: Phi = 1.00 (SE = 0) Profile CI = (1.00 - 1.00)
10: Phi = 1.00 (SE = 0) Profile CI = (1.00 - 1.00)

Cloning the data from the profile likelihood estimates also did not change the 0 standard error for parameters 1, 6 and 10. Given that the profile likelihood CI are poorly estimated, even with simulated annealing, I conclude these parameters are extrinsically non-identifiable because of sparse data.

I would concur with your conclusions, based on your solid summary.

<reference to symmetry of posterior>

The posterior distribution for mu is unimodel and quite symmetrical. The mean for mu was -0.0602, the median -0.0221, and the mode 0.1671. Back-transformed these correspond to 0.484, 0.494, 0.541, respectively.

Symmetrical? Perhaps...

Yes indeed: the posterior distributions for the boundary parameters in particular were not unimodal or symmetric. Instead, parameter 1 was multi-modal, and parameters 6 and 10 were strongly skewed. These parameters are clearly difficult to estimate.

Basically, difficult -> not possible.

I plotted the posterior density for mu on the logit scale - it appeared to be unimodel and symmetrical. At this stage I am not sure how to get a posterior density plot for mu on the real scale?

Take the chain on the logit scale, back-transform the whole thing to the real scale, and go from there. That's basically all you need to do.

The the 95% HPD for mu on the logit scale is [-1.809; 1.265]. Back transformed, this corresponds to [0.14 - 0.78] if I am correct? (I simply back-transformed the two HPD from the logit scale to the real scale). In this case the HPD for mu will (just) bound the values from the MoM and 'dot' model approach.

Simplest approach is to take the chain, back-transform to real, and construct 95% HPD on the real.

Generally, if you have sparse data, which you acknowledge, and which is consistent with multiple boundary issues, etc, there are limits to what you can do (cue reference to Darryl MacKenzies reminding us that: 'these methods are statistical, not magical'). You have sparse data. You have a bunch of estimates near the boundary. None of the covariates or anything 'biologically interesting' have any model-based support. Your choice then might be to (i) re-calculate mu and sigma excuding the years with boundary estimate problems, and/or (ii) generate a pretty picture with the estimates, indicating in some fashion which estimates might be speculative/outright wrong, and let the 'reader' decide for themselves what to make of the estimates that are well-estimated. There is always option (iii) regret stepping up to help colleagues without looking at their data first, beg forgiveness for not being able to rescue them, and consume beer.

[chris_oosthuizen]

should already raise red flags about the quality of the data -- I'm guessing low encounter probabilities?

Yes. In this case capture probabilities are roughly 0.22 for trap "unaware" individuals, and 0.36 for trap "aware" individuals. I have worked with worse (e.g., dolphin photo-ID recapture rates...) but such probabilities make me miss my usual data set with near perfect detection!

If you assume survival is a constant, and doesn't vary with age (both pretty tenuous assumptions), then there ways to calculate 'mean time to event' (being the death event), but given the assumptions, I personally wouldn't hold much faith in this approach.

Very true!

Your choice then might be to...

The approach I thought to take was to acknowledge the estimation trouble with the time-dependent model, and to simply obtain model averaged estimates (i.e., drop random effect models (ML or MCMC) from analysis altogether). I am partly motivated in this choice by what I read in Francis and Saurola (2009) (given below) and because of the evident limitations of the data.

Francis and Saurola compared survival estimates from Bayesian Hierarchical and Maximum-Likelihood Methods. They also had a number of boundary parameters, and concluded that "In these cases, because the standard error was estimated at zero, the method of moments estimator did not perform well".

Further, they concluded: "In contrast, the hierarchical models did not have problems with boundary estimates, because the prior distributions for families of parameters provided additional information, in an empirical Bayes framework for estimating parameters with few data."

I also found that the MCMCM point estimates (and 95% HPD) were shifted away from the boundary. But given the distribution of the posteriors (see earlier message) I have little confidence in some of these results.

[Francis CM, Saurola P (2009) Estimating demographic parameters from complex data sets: a comparison of Bayesian hierarchical and maximum-likelihood methods for estimating survival probabilities of tawny owls, Strix aluco in Finland. In: Thomson DL, Cooch EG, Conroy MJ, Patil GP (eds) Modeling demographic processes in marked populations. Environmental and Ecological Statistics Series Vol 3. Springer, New York, NY, pp 617−637.]

Currently, when model averaging, the time-dependent model is "dropped from the model averaging because one or more beta variances are not positive" (RMark warning). I have not tried to fix any of the boundary estimates, which may allow this model to be used in model averaging. Though I would not know what value to fix the specific survival estimates to. Since the data describes survival of a small rodent, the "estimated" phi = 1 boundary parameters is unrealistic (even over the relatively short intervals between capture occasions).

There is always option (iii) regret stepping up to help colleagues without looking at their data first, beg forgiveness for not being able to rescue them, and consume beer.

One can only hope that the colleagues will be paying for the beer! But on the positive side, this was my first foray into Method of Moments and MCMC analysis in MARK, so I did learn many things (also about RMark!) which will be of use in future.