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.

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")
```