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 that applied to my data and in the console is as follows:> `desman.y` <- structure(c(3L,4L,3L,2L,1L), .Names = c("1", "2", "3", "4", > "5")) > > RN<-function(y=desman.y,J=5,nsites=39,Nmax=100){+ #### + lik<-function(parms){ + r<-expit(parms[1]) + lambda<-exp(parms[2]) + pvec<-1-(1-r)^(0:Nmax) + gN<-dpois(0:Nmax,lambda) + gN<-gN/sum(gN) + lik<-rep(NA,nsites) + for(i in 1:nsites){ + lik[i]<-sum(dbinom(y[i],J,pvec)*gN) + } + -1*sum(log(lik)) + } + #### + + tmp<-nlm(lik,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) + } but when I try to see the results (i.e., the list in the function), I've the follow errors:> summary(RN)Error in object[[i]] : object of type 'closure' is not subsettable or obtains something like:> print(RN)function(y=desman.y,J=5,nsites=39,Nmax=100){ #### lik<-function(parms){ r<-expit(parms[1]) lambda<-exp(parms[2]) pvec<-1-(1-r)^(0:Nmax) gN<-dpois(0:Nmax,lambda) gN<-gN/sum(gN) lik<-rep(NA,nsites) for(i in 1:nsites){ lik[i]<-sum(dbinom(y[i],J,pvec)*gN) } -1*sum(log(lik)) } #### tmp<-nlm(lik,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) }>this second one seems like the analyses are not done at all. Any suggestion? Many thanks, Pablo -- View this message in context: http://r.789695.n4.nabble.com/problem-running-a-function-tp3390117p3390117.html Sent from the R help mailing list archive at Nabble.com.
Try RN() - Phil Spector Statistical Computing Facility Department of Statistics UC Berkeley spector at stat.berkeley.edu On Sat, 19 Mar 2011, garciap wrote:> 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 > > that applied to my data and in the console is as follows: >> `desman.y` <- structure(c(3L,4L,3L,2L,1L), .Names = c("1", "2", "3", "4", >> "5")) >> >> RN<-function(y=desman.y,J=5,nsites=39,Nmax=100){ > + #### > + lik<-function(parms){ > + r<-expit(parms[1]) > + lambda<-exp(parms[2]) > + pvec<-1-(1-r)^(0:Nmax) > + gN<-dpois(0:Nmax,lambda) > + gN<-gN/sum(gN) > + lik<-rep(NA,nsites) > + for(i in 1:nsites){ > + lik[i]<-sum(dbinom(y[i],J,pvec)*gN) > + } > + -1*sum(log(lik)) > + } > + #### > + > + tmp<-nlm(lik,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) > + } > > but when I try to see the results (i.e., the list in the function), I've the > follow errors: > >> summary(RN) > Error in object[[i]] : object of type 'closure' is not subsettable > > or obtains something like: > >> print(RN) > function(y=desman.y,J=5,nsites=39,Nmax=100){ > #### > lik<-function(parms){ > r<-expit(parms[1]) > lambda<-exp(parms[2]) > pvec<-1-(1-r)^(0:Nmax) > gN<-dpois(0:Nmax,lambda) > gN<-gN/sum(gN) > lik<-rep(NA,nsites) > for(i in 1:nsites){ > lik[i]<-sum(dbinom(y[i],J,pvec)*gN) > } > -1*sum(log(lik)) > } > #### > > tmp<-nlm(lik,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) > } >> > > this second one seems like the analyses are not done at all. > > Any suggestion? Many thanks, > > Pablo > > -- > View this message in context: http://r.789695.n4.nabble.com/problem-running-a-function-tp3390117p3390117.html > Sent from the R help mailing list archive at Nabble.com. > > ______________________________________________ > R-help at 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. >
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)