A Question

Recently, a fishR user asked me the following question:

I have recently used the FSA package to calculate Zippin depletion population estimates resulting from 3 pass removals on smaller trout streams. I have currently been doing this the “long way”:

##GRYSUS2012
GRYSUS2012 <- c(31,23,9)
a <- removal(GRYSUS2012)
summary(a)
confint(a)

##GRYSUS2013
GRYSUS2013 <- c(60,12,9)
b <- removal(GRYSUS2013)
summary(b)
confint(b)

##GRYSDS2012
GRYSDS2012<- c(45,12,8)
c <- removal(GRYSDS2012)
summary(c)
confint(c)

However, I was wondering if there was a way I could format my data (in tabular format) so that I may run zippin estimates on large numbers of sites at once, using a table similar to the one shown below (Note: I deleted this table but it looks like the data frame in the code below)? Any help you can provide is much appreciated. Thanks!

An Answer

There is a similar example at the very bottom of the help page for removal(), though that example is more complicated. Here is an answer to this question.

Begin by loading the FSA package …

library(FSA)

… and reading in an external tab-delimited text file that simply contains the table that the questioner provided …

( df <- read.table("shanktest.txt",header=TRUE) )

##         Site Pass_1 Pass_2 Pass_3
## 1 GRYSUS2012     31     23      9
## 2 GRYSUS2013     60     12      9
## 3 GRYSDS2012     45     12      8

Using apply() once can send each row of the data frame to removal(). The tricks here are that the first column of df must be excluded (hence the -1 in the square brackets) and removal() should return just the parameter estimates (i.e., the use of just.ests=TRUE; see the help for removal()).

( res <- apply(df[,-1],MARGIN=1,FUN=removal,just.ests=TRUE) )

##           [,1]     [,2]    [,3]
## No    76.00000 83.00000 68.0000
## p      0.44056  0.69231  0.6373
## No.se  9.05884  2.09625  2.6939
## p.se   0.09387  0.05683  0.0696

The resulting matrix, however, should be transposed so that the parameter estimates form the columns and the original groupings will form the rows. In addition, the matrix should be converted to a data frame so that the group names can be appended.

( res <- data.frame(t(res)) )

##   No      p No.se    p.se
## 1 76 0.4406 9.059 0.09387
## 2 83 0.6923 2.096 0.05683
## 3 68 0.6373 2.694 0.06960

The original groupings (from the first column of df) can then be appended to make the results more clear.

( res <- cbind(site=df[,1],res))

##         site No      p No.se    p.se
## 1 GRYSUS2012 76 0.4406 9.059 0.09387
## 2 GRYSUS2013 83 0.6923 2.096 0.05683
## 3 GRYSDS2012 68 0.6373 2.694 0.06960

Finally, the approximate 95% CIs can be appended.

res$No.LCI <- res$No-1.96*res$No.se
res$No.UCI <- res$No+1.96*res$No.se

The resulting data frame contains the required information.

res

##         site No      p No.se    p.se No.LCI No.UCI
## 1 GRYSUS2012 76 0.4406 9.059 0.09387  58.24  93.76
## 2 GRYSUS2013 83 0.6923 2.096 0.05683  78.89  87.11
## 3 GRYSDS2012 68 0.6373 2.694 0.06960  62.72  73.28

A Modified Answer

I suppose the questioner could wrap all of this into a function to make the call very easy (though this is not a very general function).

multRemoval <- function(df,colwgrp,type) {
  # df is dataframe with the pass data in columns, groups in rows
  # colwgrp is the column numbers with the group names
  # type is the type of calculation to use (see removal())
  tmp <- apply(df[,-colwgrp],MARGIN=1,FUN=removal,type=type,just.ests=TRUE)
  tmp <- data.frame(t(tmp))
  tmp <- cbind(df[,colwgrp],tmp)
  names(tmp)[1] <- names(df)[colwgrp]
  tmp$No.LCI <- tmp$No-1.96*tmp$No.se
  tmp$No.UCI <- tmp$No+1.96*tmp$No.se
  tmp
}

multRemoval(df,1,"Zippin")

##         Site No      p No.se    p.se No.LCI No.UCI
## 1 GRYSUS2012 76 0.4406 9.059 0.09387  58.24  93.76
## 2 GRYSUS2013 83 0.6923 2.096 0.05683  78.89  87.11
## 3 GRYSDS2012 68 0.6373 2.694 0.06960  62.72  73.28