[R] Plotting a curve for a Holling Type III Functional Response

[Ben Bolker]

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?
>