From: Ben Bolker <bbolker_at_gmail.com>

Date: Mon, 04 Apr 2011 14:42:25 +0000

legend("topright",paste("p",c(6,4,2),sep="="),col=c(1,2,4),lty=1)

645,0,

675,0),

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. Received on Mon 04 Apr 2011 - 14:48:36 GMT

Date: Mon, 04 Apr 2011 14:42:25 +0000

<bialozyt <at> biologie.uni-marburg.de> writes:

*>
*

> Dear Ben,

*>
**> you answerd to Nancy Shackelford about Clarks 2Dt function.
**> Since the thread ended just after your reply,
**> I would like to ask, if you have an idea how to use this function in R
**>
*

Dear Ronald,

I got started on your problem, but I didn't finish it. I got a plausible answer to start with, but when checking the answer I ran into some trouble. Unfortunately, fitting these functions is a bit harder than one might expect ... it takes quite a bit of fussing to get a good, reliable answer.

My partly-worked solution is below.

> I defined it the following way:

You were multiplying instead of dividing by the second term (I changed it by raising the term to a negative power instead.

The lesson here: *always* do some sanity checks (graphical or otherwise) of your functions. I actually did the whole fit before I tried to plot the curves and found that they were increasing rather than decreasing ...

## fixed

clark2Dt <- function(x , p, u=1) {

(p/(pi*u))/(1+(x^2/u))^(p+1)

}

It might be preferable to define this in terms of s=sqrt(u) instead (then s would be a scale parameter with the same units as x, more easily interpretable ...

Sanity checks:

par(las=1,bty="l") ## personal preferences

curve(clark2Dt(x,p=6),from=0,to=5) curve(clark2Dt(x,p=4),col=2,add=TRUE) curve(clark2Dt(x,p=2),col=4,add=TRUE)

legend("topright",paste("p",c(6,4,2),sep="="),col=c(1,2,4),lty=1)

Grab data (in the future, if possible, please use dput(), which puts your data in the most convenient form, or write out a statement like this to define your data ...)

X <- as.data.frame(matrix(

c(15,12, 45,13, 75,10, 105,8, 135,16, 165,5, 195,15, 225,8, 255,9, 285,12, 315,5, 345,4, 375,1, 405,1, 435,1, 465,0, 495,1, 525,2, 555,0, 585,0, 615,0,

645,0,

675,0),

ncol=2,byrow=TRUE, dimnames=list(NULL,c("dist","count"))))

## assume these are traps/samples with unit size ## (if not, it will get absorbed into the "fecundity" constant

library(bbmle)

m1 <- mle2(count~dnbinom(mu=f*clark2Dt(dist,p,u),size=k),

data=X,start=list(f=20,u=10,p=5,k=2), lower=rep(0.002,4),method="L-BFGS-B")

## we get a plausible-looking fit ...

with(X,plot(count~dist,pch=16,las=1,bty="l"))
newdat <- data.frame(dist=1:700) ## overkill but harmless
lines(newdat$dist,predict(m1,newdata=newdat))

## but the coefficients look funny, especially f

coef(m1)

## tried resetting parscale but it's bogus (gets stuck at a worse likelihood) m2 <- mle2(count~dnbinom(mu=f*clark2Dt(dist,p,u),size=k),

data=X,start=list(f=20,u=10,p=5,k=2), control=list(parscale=abs(coef(m1))), lower=rep(0.002,4),method="L-BFGS-B") m3 <- mle2(count~dnbinom(mu=exp(logf)*clark2Dt(dist,exp(logp),exp(logu)), size=exp(logk)), data=X,start=list(logf=log(20),logu=log(10),logp=log(5), logk=log(2)), method="Nelder-Mead")

exp(coef(m3))

coef(m1)

summary(m1)

## hmm. Redefine in terms of s instead of u and (more importantly) ## with f = seed density at r=0 rather the

cov2cor(vcov(m1)) ## shows that f and u are horribly correlated

newclark2Dt <- function(x , p, s=1, eps=1e-70) {
d <- (1+(x/s)^2)

r <- 1/d^(p+1)

if (any(!is.finite(r))) browser()

r

}

dnbinom_pen <- function(x,mu,size,pen=1000,log=TRUE) {
mu <- rep(mu,length.out=length(x))

logval <- ifelse(mu==0 && x==0,pen*x^2,dnbinom(x,mu=mu,size=size,log=TRUE))
if (log) logval else exp(logval)

}

## needed for predict()

snbinom_pen <- snbinom

m4 <- mle2(count~dnbinom(mu=f*newclark2Dt(dist,p,s),size=k),

data=X,start=list(f=20,s=10,p=5,k=2), lower=rep(0.002,4),method="L-BFGS-B") m5 <- mle2(count~dnbinom_pen(mu=f*newclark2Dt(dist,1/(pinv),s),size=exp(logk)), data=X,start=list(f=15,s=10,pinv=100,logk=1),trace=TRUE, ## control=list(parscale=c(200,0.002,1.66,3600)), lower=rep(0.002,4),method="L-BFGS-B")

with(X,plot(count~dist,pch=16,las=1,bty="l")) newdat <- data.frame(dist=1:700) ## overkill but harmless lines(newdat$dist,predict(m1,newdata=newdat)) lines(newdat$dist,predict(m5,newdata=newdat),col=2)

> but I am not able to fit anything.

*> Do you have an idea?
**> I guess there is something wrong in my formula for Clarks 2Dt
**>
**> Thank you for reading
**>
**> Ciao
**> Ronald Bialozyt
**>
**>
*

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. Received on Mon 04 Apr 2011 - 14:48:36 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Mon 04 Apr 2011 - 15:10:28 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*