# [R] 2-parameter MLE problems

Ben Bolker bbolker at gmail.com
Wed Apr 13 02:45:19 CEST 2011

```Tyler Schartel <tes164 <at> msstate.edu> writes:

>
> Hi all,
>
> Sorry for the re-post, I sent my previous e-mail before it was complete. I
> am trying to model seroprevalence using the differential equation:
> dP/dt = beta*seronegative*.001*(seropositive)-0.35*(0.999)*(seropositive)-
>     r*seropositive.
> I would like to estimate my two parameters, beta and r, using maximum
> likelihood methods. I have included my code below:
>
> summary
>   Year   N SeroPos SeroNeg
> 1    1  75       1      74
> 2    2  12       3       9
> 3    3 139      11     128
> 4    4 178      22     156
> 5    5 203      18     185
> 6    6 244      37     207
> attach(summary)
> poisNLL=function(P){
> lambda=P*SeroNeg*0.001*SeroPos-0.35*0.999*SeroPos-P*SeroPos
> v=-sum(dpois(SeroNeg,lambda=lambda,log=TRUE))
> if (!is.finite(v)) v<- -2000000
> v
> }

Your basic problem here is that you have switched the order
of the parameters/neglected to tell R which was which.

opt1=optim(fn=poisNLL,par=c(10,.1),method='BFGS')

would have done what you were aiming for ...

but I think you've got bigger problems than that.
The analysis below shows pretty thoroughly that the answer
wants to converge on beta=52, r=0.  Are you sure your
equations make sense?

dat <- data.frame(year=1:6,N=c(75,12,139,178,203,244),
SeroPos=c(1,3,11,22,18,37))
dat <- transform(dat,SeroNeg=N-SeroPos)

calclambda <- function(beta,r,SeroPos,SeroNeg) {
beta*0.001*SeroPos*SeroNeg-0.35*0.999*SeroPos-r*SeroPos
}
poisNLL <- function(P,cdat=dat) {
lambda <- calclambda(beta=P,r=P,
SeroPos=cdat\$SeroPos,SeroNeg=cdat\$SeroNeg)
lambda <- pmax(lambda,0.001)
-sum(dpois(dat\$SeroNeg,lambda=lambda,log=TRUE))
}

poisNLL(c(beta=10,r=0.1),dat)
calclambda(beta=10,r=0.1,dat\$SeroPos,dat\$SeroNeg)
library(bbmle)
parnames(poisNLL) <- c("beta","r")
mle2(poisNLL,vecpar=TRUE,
start=list(beta=10,r=0.1),data=dat,
method="L-BFGS-B",lower=c(0,0))

with(dat,plot(year,SeroNeg,las=1,bty="l",ylim=c(0,300)))
lines(1:6,calclambda(beta=41,r=0,dat\$SeroPos,dat\$SeroNeg))

library(emdbook)
par(las=1)
cc <- curve3d(poisNLL(c(beta=x,r=y)),xlim=c(40,60),ylim=c(0,0.2),
xlab="beta",ylab="r",sys3d="image")

cc2 <- curve3d(poisNLL(c(beta=x,r=y)),xlim=c(50,55),ylim=c(0,0.02),
xlab="beta",ylab="r",sys3d="image")

> opt1=optim(poisNLL,start=c(10,.1),method='BFGS')
>
> I receive the following error from this code: "Error in optim(poisNLL, start
> = c(10, 0.1), method = "BFGS") :
>   cannot coerce type 'closure' to vector of type 'double'"
>
> Any assistance provided would be greatly appreciated!

```