[R] Plotting a curve for a Holling Type III Functional Response
Ben Bolker
bbolker at gmail.com
Mon Apr 8 15:22:24 CEST 2013
Student <mpchilds <at> middlebury.edu> writes:
>
> Hey,
>
> So I have a scatter plot and I am trying to plot a curve to fit the data
> based on a Holling Type III functional response. My function is this:
>
It's hard to see how 'size=DBH' could make sense; 'size' is
an overdispersion parameter ... I guess it's *possible*, but
I can't think of circumstances where it would actually be
what you wanted to do
nll2<-function(a,b,logk,DBH, cones) {
conefun<-(a*DBH^2)/(b^2+DBH^2)
nlls <- dnbinom(x=cones ,size=exp(logk), mu=conefun,log=TRUE)
-sum(nlls)
}
You didn't give a reproducible example, but I'm going to guess
that your data are coming from here:
library(emdbook)
dd <- with(FirDBHFec_sum,data.frame(DBH,cones=fecundity))
plot(cones~DBH,data=dd)
curve(100*x^2/(9^2+x^2),add=TRUE,col=2)
## OK, these are reasonable numbers (i.e., the curve looks reasonable)
library(bbmle)
with(dd,nll2(150,9,5,DBH,cones)) ## answer is NA: why?
conefun <- with(dd,(100*DBH^2)/(9^2+DBH^2))
## there are NA values in the data!
dd <- na.omit(dd)
with(dd,nll2(150,9,1,DBH,cones))
## now it's infinite ...
## debugging inside nll2 and inspecting:
## elements 33 36 40 41 83 126 are -Inf
## these are half-integer values, which have NB probability equal to zero
dd <- dd[dd$cones %% 1 == 0 , ]
mfit <- mle2(nll2,start=list(a=150,b=9,logk=0),data=dd,
control=list(maxit=1000))
with(as.list(coef(mfit)),
curve(a*x^2/(b^2+x^2),add=TRUE,col=4))
This whole thing would be slightly more compact as follows:
mfit2 <- mle2(cones~dnbinom(mu=a*DBH^2/(b^2+DBH^2),size=exp(logk)),
start=list(a=150,b=9,logk=0),data=dd,
control=list(maxit=1000))
newdat <- data.frame(DBH=seq(3,16,length=101))
lines(newdat$DBH,predict(mfit2,newdata=newdat),col="purple")
> and my plot is this:
>
> plot (DBH,cones)
>
> DBH is on the x-axis and cones is on the y-axis. How do I get the curve for
> my function onto the scatterplot?
>
More information about the R-help
mailing list