[R] 2-parameter MLE problems

[Ben Bolker]

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=read.delim('summary.txt',header=T)
> 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[1]*SeroNeg*0.001*SeroPos-0.35*0.999*SeroPos-P[2]*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[1],r=P[2],
                       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")
contour(cc$x,cc$y,cc$z,add=TRUE)

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


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