R Under development (unstable) (2023-12-07 r85661 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ## 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)
[1] 22.3544
> gradlikfun(p0,ddat,t0,N)
      x0   lambda    sigma    alpha 
15.42028 30.95135 19.33690 30.04404 
> 
> library(numDeriv)
> grad(likfun,p0,dat=ddat,times=t0,N=N)  ## finite differences
[1] 15.42028 30.95135 19.33690 30.04404
> ## matches!
> 
> library(bbmle)
Loading required package: stats4
> parnames(likfun) <- names(p0)
> 
> 
> m1 <- mle2(likfun,start=p0,gr=gradlikfun,data=list(times=t0,N=N,dat=ddat),
+            vecpar=TRUE)
Warning messages:
1: In dpois(dat, calc_mean(p, times, N, incid = incid), log = TRUE) :
  NaNs produced
2: In dpois(dat, calc_mean(p, times, N, incid = incid), log = TRUE) :
  NaNs produced
3: In dpois(dat, calc_mean(p, times, N, incid = incid), log = TRUE) :
  NaNs produced
> 
> 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")
There were 50 or more warnings (use warnings() to see the first 50)
> 
> m0 <- mle2(likfun,start=p0,data=list(times=t0,N=N,dat=ddat),
+            vecpar=TRUE)
Warning messages:
1: In dpois(dat, calc_mean(p, times, N, incid = incid), log = TRUE) :
  NaNs produced
2: In dpois(dat, calc_mean(p, times, N, incid = incid), log = TRUE) :
  NaNs produced
3: In dpois(dat, calc_mean(p, times, N, incid = incid), log = TRUE) :
  NaNs produced
> 
> pp0 <- profile(m0,which="lambda")
There were 50 or more warnings (use warnings() to see the first 50)
> par(mfrow=c(1,2))
> plot(pp1,show.points=TRUE)
> plot(pp0,show.points=TRUE)
Warning message:
In .local(x, ...) :
  non-monotonic profile: reverting to linear interpolation.  Consider setting std.err manually
> 
> proc.time()
   user  system elapsed 
   8.20    0.21    8.42