[R] problem running a function

[Ben Bolker]

garciap <garciap <at> usal.es> writes:

> 
> Dear people,
> 
> I'm trying to do some analysis of a data using the models by Royle & Donazio
> in their fantastic book, particular the following function:
> http://www.mbr-pwrc.usgs.gov/pubanalysis/roylebook/panel4pt1.fn
> 

[snip]

 I'm guessing you're fairly new to R, because you seem confused about
how functions work.  You don't need to redefine the function with
your particular data: instead, you 'call' the function with your
data inserted in the call.

  Below I've taken the definition of RN from the URL you provided;
I change the 'expit' function to 'plogis', which is an
equivalent function already built into R; and tried to run the RN()
function on your data.  I ran into some trouble running the function
before I realized that you had specified 'nsites' different from
the length of your y vector; I don't know under what circumstances
that would ever make sense, since the function is trying to 
calculate values from each element of y between 1 and nsites ...
Then the function did run, although it gave lots of warnings
and couldn't compute the standard error for lambda.  It also gave
a very small value for r and a very large value for lambda; I suspect
you don't have enough data to do a very good job estimating in this
case ...

RN<-function(y,J=11,nsites=length(y),Nmax=100,old=TRUE){
####
  lik<-function(parms){
    r<-plogis(parms[1])
    lambda<-exp(parms[2])
    pvec<-1-(1-r)^(0:Nmax)
    gN<-dpois(0:Nmax,lambda)
    gN<-gN/sum(gN)
    if (old) {
      likvec <- numeric(nsites)
      for(i in 1:nsites){
        likvec[i]<-sum(dbinom(y[i],J,pvec)*gN)
      }
    } else {
      likvec <- rowSums(sweep(outer(y,pvec,dbinom,size=J),2,
                           gN,"*"))
    }
    res <- -sum(log(likvec))
    res
  }
####
  tmp<-nlm(f=lik,
           p=c(0,0),
           hessian=TRUE)
  ests<-tmp$estimate
  aic<-tmp$minimum*2 + 2*length(ests)
  se<- sqrt(diag(solve(tmp$hessian)))
  list(ests=ests,se=se,aic=aic)
}

desman <- RN(y=c(3,4,3,2,1),J=5)