3
$\begingroup$

I feel like this is an extremely basic question, but have not been able to dig up an answer.

We're enumerating migrating fish, and there are some counting periods (hours) in which we were not able to count during the full period. In these cases, we've estimated the total passage according to the form $$\hat{m}=\frac{n}{p}$$ in which n represents the number of fish actually counted, and p represents the proportion of the period we were able to count.

Is there a closed-form estimate of $\hat{var}(\hat{m})$?

Feels like this should be trivial, but I'm having a hard time even wrapping my brain around which term(s) to treat as a random variable.


FURTHER EXPLORATION, thanks to @RoryT's answer...

I first jumped at the Poisson process idea like a salmon after a salmon fly. We make assumptions all the time in analyses... no big deal, right? Plus, any Francophone would love the idea of using the Poisson to model fish. #BadStatsPun

But then I started thinking about how to go about violating the Poisson process assumptions, and an obvious choice is time dependency. It seems logical that migration would follow trends or pulses on some scale, rather than truly at random.

So first, a hacked-together R function to simulate random draws from a density function defined by a portion of a sinusoidal curve...

rsin <- function(n, scale=1, min=0, max=1) {
  # n = size of output vector
  # scale = contribution (on a scale of 0 to 1) of sinusoidal process to density
  # min & max = support
  # note: x is scaled to number of full periods

  x_unif <- runif(5*n, min, max)
  d_curve <- function(x) (sin(x*2*pi)/2 + 0.5)*scale + (1-scale)

  p_accept <- d_curve(x_unif)
  x_accept <- rbinom(5*n, 1, p_accept)
  x <- x_unif[x_accept==1]
  x <- x[1:n]
}

... then a wad of duct tape to simulate total runs defined by this density, and counts from some random temporal "window" of the total run.

windowsim <- function(nsim, N, p, timebounds, scale) {
  n <- rep(NA, nsim)
  plength <- p*(timebounds[2]-timebounds[1])

  for(i in 1:nsim) {
    fishtimes <- rsin(N, min=timebounds[1], max=timebounds[2], scale=scale)
    countstart <- runif(1,timebounds[1],timebounds[2]-plength)
    n[i] <- sum(fishtimes>countstart & fishtimes<(countstart+plength))
  }
  Nhat <- n/p
  parmfrow <- par("mfrow")
  par(mfrow=c(2,1))
  hist(rsin(10000, min=timebounds[1], max=timebounds[2], scale=scale), main="time density", xlab="")
  lines(c(timebounds[1],timebounds[1]+plength),rep(par("usr")[4],2),lwd=3,lend=1)
  hist(Nhat)
  abline(v=N, lwd=2, lty=2)
  par(mfrow=parmfrow)

  return(list(mean_n_sim=mean(n), n_Pois=N*p, sd_N_sim=sd(Nhat), sd_N_Pois=sqrt(mean(n)/p*(1/p-1))))
}

First, setting scale=0 gives the sinusoidal component zero influence, and forces true temporal randomness. In this case, the Poisson process approximation performs just as advertised!

windowsim(nsim = 10000, N = 1000, p = 0.6, timebounds = c(-0.5, 
    0.5), scale = 0)
## $mean_n_sim
## [1] 599.8594
## 
## $n_Pois
## [1] 600
## 
## $sd_N_sim
## [1] 25.97738
## 
## $sd_N_Pois
## [1] 25.81686

Including a trending time component, even a small one, seems to add a surprising amount of variance.

windowsim(nsim = 10000, N = 1000, p = 0.6, timebounds = c(-0.25, 
    0.25), scale = 0.5)
## $mean_n_sim
## [1] 600.4284
## 
## $n_Pois
## [1] 600
## 
## $sd_N_sim
## [1] 103.1553
## 
## $sd_N_Pois
## [1] 25.82911

The greater the time dependence, the stronger the effect.

windowsim(nsim = 10000, N = 1000, p = 0.6, timebounds = c(-0.5, 
    0.5), scale = 1)
## $mean_n_sim
## [1] 601.3204
## 
## $n_Pois
## [1] 600
## 
## $sd_N_sim
## [1] 314.2516
## 
## $sd_N_Pois
## [1] 25.84828
$\endgroup$
2
$\begingroup$

One way would be to treat the hours observed as a sample from the total hours over the period you're interested in - each hour has a fixed but unknown level of migration, and by selecting an hour to observe you're learning that level of migration. If the hours you chose to observe are independent of the process that generates the level of migration, then you could treat the observed hours as a simple random sample without replacement, and use the form of the variance for this. The SRSWOR gives the same estimate for $m$ as in your question, and gives a variance:

$$ \left(\frac{1}{p}-1\right) \frac{H}{h-1} \sum_{i=1}^h\left(f_i^2 - \frac{1}{h} \left(\sum_{i=1}^h f_i\right)^2\right)$$

Where $H$ is the total number of hours, $h$ is the number of hours observed, and $f_i$ is the number of fish observed in the $i$th hour observed.

Another way is that you could assume the fish migration process follows a particular distribution. For example, you could model this as a Poisson Process, where the expected waiting time between each fish is the same. In this case your estimate of the total number of fish is again the formula you've presented in the question, but this time the variance of this estimate is given by:

$$\frac{n}{p}\left(\frac{1}{p}-1\right)$$

This formula is simpler because we've made stronger assumptions about the nature of the random process that is generating observations. In fact, although I am by no means a subject matter expert, I would guess that observations of fish aren't well explained by a Poisson Process.

I suppose the point of this answer is to say: variance of this estimate depends on the assumptions that underlie your model. I'd recommend bringing a statistician into your analysis if you possibly can.

$\endgroup$
3
  • $\begingroup$ Thanks! Appreciate the answer, it was really helpful - though it looks like I have more homework to do. Check out my edits - I got curious about the Poisson process idea and tried a sim. $\endgroup$ – Matt Tyers Jun 9 '17 at 20:22
  • $\begingroup$ Hey Rory - could I trouble you for a reference or a derivation? I can't give you any more points, but I'll buy you a beer next time you're in Alaska ;) $\endgroup$ – Matt Tyers Jun 13 '17 at 22:00
  • $\begingroup$ Hi Matt, the first formula is from survey methodology (it's my area) - the text this from is Sarndel, Swenson and Wretman - "Model-Assisted Survey Sampling", but I read a For the second formula I did a quick derivation myself, using the Law of Total Variance, conditioning on the Poisson rate parameter being known. I double checked, but I can't be sure I didn't make an error. In any case, I reckon you should still consult a statistician about it. $\endgroup$ – RoryT Jun 19 '17 at 6:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.