## implement richards-incidence (="revised superlogistic")
## with analytic gradients
## from Junling's code:
model_richardson <- function(times, theta, N)
{
x0 = theta[1]
lambda = theta[2]
K = theta[3] * N
alpha = theta[4]
return(K/(1+((K/x0)^alpha-1)*exp(-lambda*alpha*times))^(1/alpha))
}
## equivalent model, in terms of sigma and as a symbolic expression
Rcum <- expression((sigma*N)/(1+(((sigma*N)/x0)^alpha-1)*exp(-lambda*alpha*times))^(1/alpha))
pnames <- c("x0","lambda","sigma","alpha")
## function to compute gradient (and value), derived by R
Rderiv <- deriv(Rcum,pnames, function.arg=c(pnames,"N","times"))
## equivalent (using Rcum): return incidence (incid=TRUE) or cumulative incidence (incid=FALSE)
calc_mean <- function(p,times,N,incid=TRUE) {
## this is more 'magic' than I would like it to be ...
## have to create an environment and populate it with the contents of p (and N and times),
## then evaluate the expression in this environment
pp <- c(as.list(p),list(times=times,N=N))
## e0 <- new.env()
## mapply(assign,names(pp),pp,MoreArgs=list(envir=e0))
cumvals <- eval(Rcum,envir=pp)
if (incid) diff(cumvals) else cumvals
}
## Poisson likelihood function
likfun <- function(p,dat,times,N,incid=TRUE) {
-sum(dpois(dat,calc_mean(p,times,N,incid=incid),log=TRUE))
}
## deriv of P(x,lambda) = -sum(dpois(x,lambda,log=TRUE)) wrt lambda == sum(1-lambda/x) = N - lambda/(sum(x))
## deriv of P(x,lambda) wrt p = dP/d(lambda) * d(lambda)/dp
## compute gradient vector
gradlikfun <- function(p,dat,times,N,incid=TRUE) {
gcall <- do.call(Rderiv,c(as.list(p),list(times=times,N=N))) ## values + gradient matrix
lambda <- gcall
attr(lambda,"gradient") <- NULL
if (incid) lambda <- diff(lambda)
gmat <- attr(gcall,"gradient") ## extract gradient
if (incid) gmat <- apply(gmat,2,diff) ## differences
totderiv <- sweep(gmat,MARGIN=1,(1-dat/lambda),"*") ## apply chain rule (multiply columns of gmat by dP/dlambda)
colSums(totderiv) ## deriv of summed likelihood = sum of derivs of likelihod
}
N <- 1000
p0 <- c(x0=0.1,lambda=1,sigma=0.5,alpha=0.5)
t0 <- 1:10
## deterministic versions of data (cumulative and incidence)
dcdat <- model_richardson(t0,p0,N)
ddat <- diff(dcdat)
plot(t0,dcdat)
plot(t0[-1],ddat)
set.seed(1001)
ddat <- rpois(length(ddat),ddat)
likfun(p0,ddat,t0,N)
gradlikfun(p0,ddat,t0,N)
library(numDeriv)
grad(likfun,p0,dat=ddat,times=t0,N=N) ## finite differences
## matches!
library(bbmle)
parnames(likfun) <- names(p0)
m1 <- mle2(likfun,start=p0,gr=gradlikfun,data=list(times=t0,N=N,dat=ddat),
vecpar=TRUE)
plot(t0[-1],ddat)
lines(t0[-1],calc_mean(coef(m1),times=t0,N=N))
if (FALSE) {
## too slow ..
pp0 <- profile(m1)
pp0C <- profile(m1,continuation="naive")
}
pp1 <- profile(m1,which="lambda")
m0 <- mle2(likfun,start=p0,data=list(times=t0,N=N,dat=ddat),
vecpar=TRUE)
pp0 <- profile(m0,which="lambda")
par(mfrow=c(1,2))
plot(pp1,show.points=TRUE)
plot(pp0,show.points=TRUE)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.