R code for fitting a loglinear-NHPP model to grouped temporal data

The following R code fits a NonHomogeneous Poisson Process (NHPP) model to grouped temporal data. Grouped temporal data occurs when only the number of recurrences within a given time interval are known.
The R code fits a NHPP model with a loglinear intensity function, where the intensity at time t is given by .

See this post for fitting a NHPP model with a power law intensity function to grouped temporal data.

Suggestions and/or questions? Please contact Stefan Gelissen (email: info at datall-analyse.nl).
See this page for an overview of all of Stefan’s R code blog posts.

R code

library(car)

##data
#generate temporal data for a nonhomogeneous Poisson process
#having a loglinear intensity function
#(with the intensity at time t defined as: exp[gamma0 + gamma1*t] )

#function for computing the intensity measure of a loglinear-NHPP
intensityMeasure <- function (gamma0, gamma1, t) {
  exp(gamma0)*(exp(gamma1*t)-1)/gamma1}

#function for computing the inverse of the specified intensity measure
intensityMeasureInv <- function(gamma0, gamma1, s) {
  (1/gamma1)*log(gamma1*exp(-gamma0)*s + 1)}

Tend <- 26 #upper bound of observation window
t <- 0 #start of observation window
i <- 2 #initialize counter

#population values for gamma0 and gamma1 of the intensity measure
gamma0 <- -1.5
gamma1 <- 0.15

#sequence of event times generated by the loglinear-NHPP
while (t[length(t)] <= Tend) {
  Ti <- t[i-1]
  s <- intensityMeasure(gamma0, gamma1, Ti) - log(runif(1,0,1))
  Tnew <- intensityMeasureInv(gamma0, gamma1, s)
  t <- c(t, Tnew)
  i <- i+1
}
#remove t=0 and the final event time, since the latter time > Tend
t <- t[-c(1,length(t))]

#plot generated data
plot(t, 1:length(t), type="p", xlim=c(0,Tend),
     xlab="Time", ylab="Cumulative recurrences")



#create grouped temporal data

#construct time intervals
int <- seq(0, 25, 5)
#lower and upper bound of each time interval
tL <- int #lower bounds
tU <- c(int[-1], max(t)) #upper bounds

#plot the generated data with the boundaries of the time intervals included
plot(t, 1:length(t), type="p", xlim=c(0,Tend),
     xlab="Time", ylab="Cumulative recurrences")
abline(v=tU, lty=2, col="blue")

#assign observations to the time intervals
assignT <- findInterval(t, int)

#number of recurrences within each time interval
recurrences <- sapply(1:length(tU), function (i) sum(assignT==i))





##fit loglinear-NHPP model to grouped temporal data

#function for computing the log-likelihood of a grouped loglinear-NHPP model
llikNHPPg <- function (theta, LowerTimes, UpperTimes, Ns) {
  gamma0 <- theta[1]
  gamma1 <- theta[2]
  mu <- (exp(gamma0+gamma1*UpperTimes)-exp(gamma0+gamma1*LowerTimes))/gamma1
  sum(Ns*log(mu) - mu)
}

#function for fitting a grouped loglinear-NHPP model
groupedNHPP <- function (LowerTimes, UpperTimes, Ns, starts=NULL) {
  #-LowerTimes is a vector with the lower values of the time intervals
  #-UpperTimes is a vector containing the upper values of the time intervals
  #-Ns is a vector with the total number of recurrences within each time interval
  #-starts is a vector with the starting values for gamma0 and gamma1, respectively
  
  #obtain starting values for gamma0 and gamma1
  if (is.null(starts)) {
    iwd <- UpperTimes-LowerTimes #width of time interval
    midp <- (UpperTimes-LowerTimes)/2 #midpoint of time interval
    ti <- LowerTimes+midp
    modstart <- glm(Ns~ti, family=quasi(link="log", variance="mu"), offset=log(iwd))
    starts <- as.vector(coef(modstart)) #starting values
  }
  
  #optimization of log-likelihood
  gloglin <- optim(starts, llikNHPPg,
                   method="BFGS", control=list(fnscale=-1),
                   LowerTimes=LowerTimes, UpperTimes=UpperTimes, Ns=Ns,
                   hessian=TRUE)
  parameters <- gloglin$par
  names(parameters) <- c("gamma0", "gamma1")
  loglik <- gloglin$value
  hessian <- gloglin$hessian
  
  #return values
  list(parameters=parameters, loglik=loglik, AIC=-2*loglik+2*length(parameters),
       hessian=hessian, LowerTimes=LowerTimes, UpperTimes=UpperTimes,
       Ns=Ns, starts=starts)
}



#fit loglinear-NHPP to grouped temporal data
(ppgLoglin <- groupedNHPP(LowerTimes=tL, UpperTimes=tU, Ns=recurrences))[1:3]

#95% confidence intervals (normal-approximation) for gamma0 and gamma1
#of the fitted loglinear-NHPP model
SEparms <- sqrt(diag(solve(-1*ppgLoglin$hessian)))
ppgLoglin$parameters[1] + c(1,-1)*qnorm(.05/2)*SEparms[1] #95% ci for gamma0
ppgLoglin$parameters[2] + c(1,-1)*qnorm(.05/2)*SEparms[2] #95% ci for gamma1
#compare these with the population values for gamma0 and gamma1
#remember, the population values for gamma0 and gamma1 were -1.5 and 0.15, respectively



#prediction of the number of recurrences:
#observed/predicted cumulative number of recurrences at the end of a time interval

upperBoundInterval <- tU #upper bounds of time intervals
#compute predicted number of recurrences (and their standard errors)
#at the end of each time interval
SEs <- sapply(1:length(upperBoundInterval),
              function (i) deltaMethod(ppgLoglin$parameters,
                                       g="(exp(gamma0)*(exp(gamma1*ts)-1))/gamma1",
                                       vcov.=solve(-1*ppgLoglin$hessian),
                                       constants=c(ts=upperBoundInterval[i])))
#95% confidence interval for the predicted number of recurrences
llci <- unlist(SEs[1,])*exp((qnorm(.025)*unlist(SEs[2,]))/unlist(SEs[1,]))
luci <- unlist(SEs[1,])*exp((qnorm(.975)*unlist(SEs[2,]))/unlist(SEs[1,]))
#plot observed cumulative number of recurrences at the end of each time interval
plot(upperBoundInterval, cumsum(recurrences),
     xlab="Upper bound time interval", ylab="Cumulative recurrences",
     ylim=c(min(llci), max(luci)))
#add predicted cumulative number of recurrences to the plot
points(upperBoundInterval, unlist(SEs[1,]), col="red")
#add error bars representing the 95% confindence intervals
arrows(upperBoundInterval, llci, upperBoundInterval , luci,
       code=3, length=0.1, angle=90, col="red")
legend("topleft", pch=1, col=c("black", "red"),
       legend=c("Observed", "Predicted"))





#finally, explore the log-likelihood function in the neighborhood
#of the starting values that were used for fitting the loglinear-NHPP model
#more specifically, locate the positions of both the starting values and the
#MLEs for gamma0 and gamma1 in a contour plot of log-likelihood values

#construct grid and calculate log-likelihood values
gamma0values <- seq(-3.5,-.5,length.out=75)
gamma1values <- seq(.05,.25,length.out=75)
lbgrid <- expand.grid(gamma0=gamma0values, gamma1=gamma1values)
lbgrid$ll <- apply(lbgrid, 1, llikNHPPg, LowerTimes=tL,
                   UpperTimes=tU, Ns=recurrences)
ll.matrix <- matrix(lbgrid$ll, nrow=length(gamma0values))

#contour plot of log-likelihood values
contour(gamma0values, gamma1values, ll.matrix,
        levels=seq(max(ll.matrix)+55, min(ll.matrix), length.out=20),
        xlab=bquote(gamma[.(0)]), ylab=bquote(gamma[.(1)]))
#position of starting values for gamma0 and gamma1 (blue cross)
points(x=ppgLoglin$starts[1], y=ppgLoglin$starts[2], col="blue", pch=4)
#position of MLEs for gamma0 and gamma1 (red dot)
points(x=ppgLoglin$parameters[1], y=ppgLoglin$parameters[2], col="red", pch=20)