mebrooks / selfisher

selfisher - modeling the SELectivity of FISHERies gear in R
9 stars 4 forks source link

possible error in logit(phi) eqn. in selfisher.cpp #31

Closed wStockhausen closed 4 years ago

wStockhausen commented 4 years ago

In your TMB code, you have the following code (lines 524-527 in selfisher.cpp) for logit(phi):

} else { //covered codend or catch comparison for (int i = 0; i < r.size(); i++) { logit_phi(i) = log(q(i)) + logit_inverse_linkfun(etar(i), etad(i), link); //log(q) + logit(r) }

which I think is wrong--I think the term involving etar(i)... should evaluate to log(r), not logit(r) [see detail provided below]. Thus, the term should be "log(inverse_linkfun(etar(i), etad(i), link))", or just "log(r(i))" since you already computed r(i)=inverse_linkfun(etar(i), etad(i), link) back on line 514. However, it's also possible I'm completely off track!

Results from a similar comparison between AFSC and BSFRF side-by-side catches for snow crab by Somerton et al (2013: CJFAS 70:1699-1708) were based on fitting the equation logit(phi) = log(r(z)) + log(q) [not logit(r(z)) + log(q)] which was developed from phi = (S_apr(z))/[S_apr(z)+S_b(1-p)] where z is size, r(z) is the fraction of individuals entering the AFSC gear that were retained, S_a is the sampling fraction for the AFSC gear, S_b is the sampling fraction for the BSFRF gear, and p is the assumed fractional split for the AFSC gear of the total number of crab entering the combined gear. p for each haul was assumed to be the fraction of the total area swept by the AFSC gear (i.e., p=A_a/[A_a+A_b] with A_a = area swept by AFSC gear, A_b=area swept by BSFRF gear). Setting q=A_bS_b/(A_a*S_a) in the equation for phi and taking the logit of both sides gets you to the first equation.

Buck Stockhausen (william.stockhausen@noaa.gov)

mebrooks commented 4 years ago

This model should be fit with psplit=TRUE, i.e. using the "trouser trawl or alternate haul" model rather than the code quoted above for covered codend or catch comparison.

If you fix p=0.5, then the "trouser trawl or alternate haul" model (logit_phi = log(q) + log(p) + log(r) - log(Type(1.0)-p)) simplifies to logit_phi = log(r) + log(q). To do that, use p=~0 in the selfisher call as in the example below

dat <- transform(haddock, tot=nfine+nwide, prop=nwide/(nfine+nwide))
m0 <- selfisher(prop~Lengths, pformula=~0, psplit=TRUE, total=tot, dat)

If you wanted to fix p at another value based on prior knowledge of area swept, then it will be more difficult, but it is possible.

mebrooks commented 4 years ago

Also, I welcome any suggested improvements for the documentation. It's possible that it's written in terms that are more familiar to European researchers, but we aim to be global.

mebrooks commented 4 years ago

I'm going to close this for now since it isn't an error. Feel free to suggest documentation updates or request that we make it easier to fix p at values other than 0.5.

wStockhausen commented 4 years ago

Hi Molly,

Thanks! Having worked on this for a week, I now understand your original comments above regarding the use of psplit and setting formula to "~0" in my case. Took awhile since "p" was not obviously part of the formula I was fitting--it was buried in the ratio of areas swept by the two side-by-side gears.

Actually, though, it's quite easy to incorporate fixed values for "p" other than 0.5---that's essentially what I'm doing in my case and why "p" isn't in the formula I'm fitting.You just incorporate the fixed value for "p" into the "qratio": new qratio = your qratio/(p/(1-p)) and again set pformula = ~0 [or new qratio = your qratio*(p/(1-p)), depending on which way you define the ratio of areas swept). And you can do this on a per-haul basis, which is what I need to do, because qratio is defined on a per haul basis.

Cheers,

Buck


On Wed, Jun 10, 2020 at 7:58 AM Mollie Brooks notifications@github.com wrote:

I'm going to close this for now since it isn't an error. Feel free to suggest documentation updates.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/mebrooks/selfisher/issues/31#issuecomment-642065429, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABMAD5AZEUCZMYQL4VSNR43RV6NQPANCNFSM4NQDRFUQ .

mebrooks commented 4 years ago

That's a cool trick! I'll have to keep it in mind in case it comes up again. I'm glad it's documented here on GitHub now. Thanks!