cschwarz-stat-sfu-ca / BTSPAS

Bayesian Time Stratified Petersen Analysis System
1 stars 3 forks source link

Prior for betas for logit P should be passable from call to _fit routines #27

Closed cschwarz-stat-sfu-ca closed 3 years ago

cschwarz-stat-sfu-ca commented 3 years ago

Tyler Pilger, PhD Fisheries Biologist

I have been using the BTSPAS package in r to estimate annual juvenile salmon abundance from rotary screw trap data. The mark-recapture methods include collecting fish from a single trap, marking them, then transporting and releasing them upstream of the trap to estimate trap efficiency (p). There are circumstances that prevent weekly releases of marked fish which has led to some sparse empirical estimates of p within each year. For an average 24 week season, there are anywhere from four to 13 releases in a given season. However, I have 14 years of mark-recapture releases and these estimates fit somewhat to linear regression where logit(p) ~ intercept + slope*log(flow).

I have followed the instructions in the BTSPAS vignette for including flow as a covariate where I created a two column matrix of 1s and weekly mean flow and passed this to the TSPDE_fit function. When I plot the estimated logit(p)'s against log(flow) it suggests that there is no relationship with flow because there is high variability and low number of empirical p estimates in a given year.

I tried to use the same covariate approach but instead of a matrix of 1s and mean flow, I used the fitted intercept and slope*log(mean flow) where slope and intercept were from all years, but this did not affect the estimated logit(p)s. The resulting plot looked the same.

Is there a way to use the predicted relationship of logit(p)~intercept + slope*log(flow) across all years to help inform the estimate of p within a year and how would I implement it?

cschwarz-stat-sfu-ca commented 3 years ago

This may be possible using a "trick". For example, if you look at the actual code for the diagonal case found in the TimeStratPetersenDiagError.R code, you will see the priors for the intercept and the coefficients associated with the other covariates ...

Capture probabilities covariates

beta.logitP[1] ~ dnorm(mu_xiP,tau_xiP) # first term is usually an intercept for(i in 2:NlogitP.cov){ beta.logitP[i] ~ dnorm(0, .01) # rest of beta terms are normal 0 and a large variance }

... At the moment the first logit parameter has a prior (on the logit scale) based on mean mu_xiP and 1/var of tau_xiP. The rest of the logitP are hard coded to have a mean of 0 and a 1/var of .01.

The mu_xiP and tau_xiP are parameters in the call to the TimeStratPetersenDiagError_fit so you should be able to set these in the function call, but you cannot change the priors on the other covariates. Unfortunately, I see that I don't pass these parameters internally, so I would need to fix that first before doing the trick below...(rats...but a simple fix)

The trick involves switching the covariate matrix to put the log(flow) first and the intercept second. Then call TimeStratPetersenDiagError_fit and set the mu_xiP and tau_xiP parameter to give a prior for the first coefficient (which is now the slope for log(flow)). where mu_xiP will be the prior slope you want (on the logit scale) with tau_xiP=1/var controls how sure you are of the prior, e.g. try different tau_xiP's.

If you would not mind me sending a sample dataset, I can also try it out here. I can then post the updated code on the GitHub site that you could download until the changes get propagated to CRAN....

If I'm modifying the code, I likely should also modify the code to allow the user to set the priors for all the covariates for logitP and then yu don't have to resort to the trick above (which won't work at the moment) . It appears to be a relatively simple fix (as he says naively) that I can try and implement for all the models in BTSPAS.

cschwarz-stat-sfu-ca commented 3 years ago

You can now add prior to the beta coefficients for logitP vs covariates using the prior.beta.logitP.mean and prior.beta.logitP.sd. Refer to the vignette on the diagonal case. The other routines work similarly.

Carl Schwarz