search

fishR-Core-Team / FSA

FSA (Fisheries Stock Assessment) package provides R functions to conduct typical introductory fisheries analyses.
https://fishr-core-team.github.io/FSA/
GNU General Public License v2.0
66 stars 22 forks source link

Inf in SE from mrOpen() #69

Closed droglenc closed 3 years ago

droglenc commented 3 years ago

Some of the SEs in mrOpen() will return as Inf. This appears to happen if any of the r values (an intermediate calculation of the future catch of fish released at time i) are 0. For example, the SE of M is computed with (line 272 in mrOpen() ....

df$M.se <- sqrt((df$M-df$m)*(df$M-df$m+df$R)*((1/df$r)-1/df$R))

and the SE of N is computed with (line 290 in mrOpen()) ...

df$N.se <- sqrt(df$N*(df$N-df$n)*((df$M-df$m+df$R)/df$M*((1/df$r)-(1/df$R))+((df$N-df$M)/(df$N*df$m))))

In both cases, for values of df$r that are 0 then the 1/df$r will return Inf.

An R default warning is already issued when this happens, but we should implement a more informative warning. The most logical place for this is inside of the iCalcr internal function (begins on line 242 in mrOpen()).

This code demonstrates the issue...

mb.top <- matrix(c(NA,NA,NA,NA,
                   1,NA,NA,NA,
                   0,5,NA,NA,
                   1,1,0,NA),nrow=4)

mb.bot <- matrix(c(0,3,3,3,
                   1,4,5,5,
                   5,0,5,5,
                   2,1,3,3),nrow=4)
rownames(mb.bot)<-c("m","u","n","R")

mr1 <- mrOpen(mb.top,mb.bot)
summary(mr1,verbose=TRUE)

Relatedly, note that there will be an issue with the SE of M if R is 0 (no marked fish are returned) or if (1/r-1/R)<0 (more marked fish are recaptured in the future than were returned). Probably should issue more informative warnings for these as well.

droglenc commented 3 years ago

This now "looks" like the following in v0.9.0.

> mb.top <- matrix(c(NA,NA,NA,NA,
+                    1,NA,NA,NA,
+                    0,5,NA,NA,
+                    1,1,0,NA),nrow=4)
> 
> mb.bot <- matrix(c(0,3,3,3,
+                    1,4,5,5,
+                    5,0,5,5,
+                    2,1,3,3),nrow=4)
> rownames(mb.bot)<-c("m","u","n","R")
> 
> mr1 <- mrOpen(mb.top,mb.bot)
Warning message:
Some SE for M and thus N will be infinity because some r=0; i.e., future recaptures for some times are 0. 
> summary(mr1,verbose=TRUE)
Observables:
  m n R  r  z
1 0 3 3  2 NA
2 1 5 5  6  1
3 5 5 5  0  2
4 2 3 3 NA NA

Estimates (phi.se includes sampling and individual variability):
     M M.se    N N.se   phi phi.se   B B.se
1   NA   NA   NA   NA 0.619  0.245  NA   NA
2  1.9  NaN  5.6  1.3 2.902    Inf 0.8  Inf
3 17.0  Inf 17.0  Inf    NA     NA  NA   NA
4   NA   NA   NA   NA    NA     NA  NA   NA