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!