R code for fitting a nonhomogeneous temporal Poisson process model using the spatstat package

In my previous blog post I showed how to fit power and loglinear nonhomogeneous Poisson process (NHPP) models to temporal data.

The following R code fits these same two models, but this time using the spatstat package.

NHPP spatstat

The data in Example 1 of the R code below was taken from Meeker and Escobar’s 1998 book Statistical methods for reliability data.
The data used for Example 2 can be found in: Guo, H., Mettas, A., Sarakakis, G., and Niu, P. (2010), Piecewise NHPP Models with Maximum Likelihood Estimation for Repairable Systems, Reliability and Maintainability Symposium (RAMS), 2010 Proceedings.

An excellent and comprehensive introduction to the spatstat package is given by Baddeley, Rubak, and Turner in their 2015 book Spatial point patterns: methodology and applications with R.

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(spatstat)


##EXAMPLE 1: fit a power and loglinear-NHPP model to temporal data

##data

#U.S.S. Halfbeak diesel engine data from Meeker & Escobar (1998), Example 16.12, pp. 415-416
halfb <- read.table("http://www.public.iastate.edu/~wqmeeker/anonymous/Stat533_data/Splida_text_data/halfbeak.txt", header=T)
head(halfb)

#plot of time against cumulative recurrences
ct <- halfb$Hours[2:72]
cf <- 1:length(ct)
plot(ct, cf, type="p", xlab="Time", ylab="Cumulative recurrences")



##set up data for the spatstat package

#observation window
Tstart <- 0 #lower bound of observation window
Tend <- halfb$Hours[73] #upper bound of observation window
v <- ppp(x=c(Tstart, Tend), y=c(0,0), xrange=c(Tstart, Tend), yrange=c(0,0))
edg <- matrix(c(1,2), ncol=2)
Tw <- linnet(v, edges=edg)

#insert event times into observation window
times <- halfb$Hours[2:72]
t <- list(x=times, y=rep(0,length(times)))
tmp <- lpp(t, Tw)

#plot the event times within the observation window
plot(tmp, main=NULL)



##fit NHPP model to temporal data

#power-NHPP model
#with the intensity at time x defined as alpha*x^k, which is
#identical to exp(theta0 + theta1*log(x)), where theta0=log(alpha)
#and theta=k
fit <- lppm(tmp~log(x), correction="none")
fit
vcov(fit) #variance covariance-matrix for the coefficients
logLik(fit) #log-likelihood
#transform the estimated coefficients such that the intensity at time x
#is defined as [beta/eta]*[x/eta]^[beta-1]
(beta <- coef(fit)[2]+1)
(eta <- (beta/exp(coef(fit)[1]))^(1/beta))
#note that these transformed coefficients are very similar to the
#estimates for beta and eta obtained in my previous post: http://blogs2.datall-analyse.nl/2016/02/17/rcode_temporal_nhpp/

#loglinear-NHPP model
#with the intensity at time x defined as exp(theta0 + theta1*x)
fit <- lppm(tmp~x, correction="none")
fit
vcov(fit) #variance covariance-matrix for the coefficients
logLik(fit) #log-likelihood
#note that the estimated coefficients and variance-covariance matrix
#are very similar to those obtained in my previous post: http://blogs2.datall-analyse.nl/2016/02/17/rcode_temporal_nhpp/





##EXAMPLE 2: fit a piecewise NHPP model to temporal data

##data

#data from Guo et al. (2010)
pw <- c(15.70, 29.39, 41.14, 56.47, 75.61, 98.83, 112.42, 125.61,
        129.39, 133.45, 138.94, 141.41, 143.67, 144.63, 144.95,
        145.16, 146.25, 146.70, 147.26, 148.15, 152.40)

#plot of time against cumulative recurrences
ct <- pw
cf <- 1:length(ct)
plot(ct, cf, type="p", xlab="Time", ylab="Cumulative recurrences")



##set up data for the spatstat package

#observation window
Tstart <- 0 #lower bound of observation window
Tend <- max(pw) #upper bound of observation window
v <- ppp(x=c(Tstart, Tend), y=c(0,0), xrange=c(Tstart, Tend), yrange=c(0,0))
edg <- matrix(c(1,2), ncol=2)
Tw <- linnet(v, edges=edg)

#event times
times <- pw

#insert event times into observation window
t <- list(x=times, y=rep(0,length(times)))
tmp <- lpp(t, Tw)

#plot the event times within the observation window
plot(tmp, legend=FALSE, main=NULL)

#Guo et al. suggested that the data should be modeled using a piecewise
#NHPP model, with a change point at t=125.5999

#create a function that represents such a change point at t=125.5999
cp <- function (x,y,seg,tp) {ifelse(x<=125.5999, 0, 1)}
changePoint <- linfun(cp, Tw)

#fit a piecewise NHPP model having a power law intensity function
fit <- lppm(tmp~log(x)*as.factor(changePoint(x,y)), correction="none")
fit
logLik(fit) #log-likelihood



##plot the fitted piecewise NHPP model

#coefficients of the piecewise NHPP model before the change point
gamma10 <- exp(coef(fit)[1])
gamma11 <- coef(fit)[2]
#predicted cumulative recurrences before the change point
ts1 <- seq(.1, 125.5999, length.out=30)
crb <- (gamma10/(gamma11+1))*ts1^(gamma11+1)

#coefficients of the piecewise NHPP model after the change point
gamma20 <- exp(coef(fit)[1]+coef(fit)[3])
gamma21 <- coef(fit)[2]+coef(fit)[4]
#predicted cumulative recurrences after the change point, and
#make sure that the predicted cumulative recurrences immediately
#before and after the change point are identical
ts2 <- seq(125.5999, 160, length.out=30)
cra <- ((gamma10/(gamma11+1))*125.5999^(gamma11+1)-
          (gamma20/(gamma21+1))*125.5999^(gamma21+1))+
  (gamma20/(gamma21+1))*ts2^(gamma21+1)

#plot predicted/observed cumulative recurrences
plot(ct, cf, type="p", pch=as.numeric(marks(tmp)), xlab="Time",
     ylab="Cumulative recurrences", main="piecewise Power-NHPP")
lines(ts1, crb, col="blue", lty=2) #predictions before the change point
lines(ts2, cra, col="blue", lty=2) #predictions after the change point
abline(v=125.5999, col="red", lty=2) #include change point



##compare fitted NHPP models with/without a change point
fit2 <- lppm(tmp~log(x), correction="none") #NHPP model without a change point
fit2
AIC(fit2) #without change point
AIC(fit) #with change point
anova(fit2, fit, test="LR") #likelihood ratio test





##EXAMPLE 3: data generated by a NHPP having some arbitrary intensity function
library(pracma)

##generate data

#function that specifies some arbitrary (loglinear) intensity function
Fintensity <- function (t) {exp(-1 + .003*t + .0015*t^2)}

#observation window
Tstart <- 0 #lower bound of the observation window
Tend <- 60 #lower bound of the observation window

#plot the specified intensity function
timesSeq <- seq(Tstart, Tend, length.out=Tend)
plot(timesSeq, Fintensity(timesSeq), type="l", xlab="Time", ylab="Intensity")

#function for generating temporal points with a NHPP having
#some arbitrary intensity function
rNHPP <- function (Tstart, Tend, intensity) {
  #-Tstart specifies the lower bound of the observation window
  #-Tend specifies the the upper bound of the observation window
  #-intensity is a function expressing the intensity function
  
  t <- Tstart
  i <- 2
  
  #function for computing the inverse of the intensity measure
  Finv <- function(Tnew, MU, U){integrate(intensity, 0, Tnew)$value - MU + U}
  
  while (t[length(t)] <= Tend) {
    Ti <- t[i-1]
    MU <- integrate(intensity, 0, Ti)$value
    U <- log(runif(1,0,1))
    Tnew <- fzero(Finv, (Tstart+Tend)/2, MU=MU, U=U)$x
    t <- c(t, Tnew)
    i <- i+1
  }
  #remove t=0 and the final event time, since the latter time > Tend
  t[-c(1,length(t))]
}

#generate temporal points
dataNHPP <- rNHPP(Tstart=Tstart, Tend=Tend, intensity=Fintensity)



##set up data for the spatstat package

#observation window
v <- ppp(x=c(Tstart, Tend), y=c(0,0), xrange=c(Tstart, Tend), yrange=c(0,0))
edg <- matrix(c(1,2), ncol=2)
Tw <- linnet(v, edges=edg)

#insert event times into observation window
times <- dataNHPP
t <- list(x=times, y=rep(0,length(times)))
tmp <- lpp(t, Tw)

#plot the event times within the observation window
plot(tmp, main=NULL)



##fit a loglinear NHPP model to temporal data

fit <- lppm(tmp~x + I(x^2), correction="none")
coef(fit) #parameter estimates
confint(fit) #confidence intervals for the parameter estimates
#remember that the true (population) values for the parameters were
#-1, .003, and .0015, respectively

#it is also possible to fit the above model using the ppm-function as follows
X <- ppp(x=dataNHPP, y=rep(1,length(dataNHPP)), xrange=c(Tstart, Tend), yrange=c(0,1))
fit2 <- ppm(quadscheme(X, ntile=c(200,1)) ~ polynom(x,2), correction="none")
coef(fit2) #parameter estimates
confint(fit2) #confidence intervals for the parameter estimates

#compare these parameter estimates with those from the lpp-function,
coef(fit)
confint(fit)
#notice that the estimates from ppm-function are reasonably close
#to the estimates from the lpp-function

#since the estimates from the ppm-function are reasonably close to those from
#the lpp-function, it seems reasonable to use some of the diagnostic tools of the
#ppm-function for model validation

#for instance, use the partial residuals plot for assessing the
#correctness of the functional form of the specified intensity function
tpar <- parres(fit2, "x")
plot(tpar, main="Time parameter")