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)
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
