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!
______________________________________________
R-help@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.
