The following R code demonstrates how to fit a multitype NonHomogeneous Poisson Process (NHPP) model to temporal data. Multitype Poisson processes belong to the family of marked point processes.

The R code may be used to fit 1) a multitype NHPP model with a loglinear intensity function, with the intensity at time t defined by , or 2) a multitype NHPP model with a power law intensity function, with the intensity at time t defined by .

Furthermore, the R code also shows how to fit a multitype Homogeneous Poisson Process (HPP) model to temporal data.

See this blog post for fitting a multitype NHPP model 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)
#The R functions for fitting a multitype temporal NHPP model
#can be downloaded from:
source("http://www.datall-analyse.nl/R/NHPP_multitype_temporal.R")
#Have a look at the function multiNHPP and notice that at the beginning
#of the script you will find a description of the function's arguments.
##EXAMPLE 1: "parallel systems"
#Data for three different types of temporal points are generated.
#The three different types of points are labeled "marks1", "marks2", and
#"marks3". The intensity of all three types is described by a loglinear
#function, with the intensity at time t defined by exp(gamma0 + gamma1*t).
#Furthermore, the observation windows for the three different types
#all start at t=0 and end at t=50.
#Finally, note that this data also can be regarded as three parallel
#operating systems. These three parallel systems all start operating at t=0,
#and stop operating at t=50. The three parallel systems are labeled "marks1",
#"marks2", and "marks3", respectively.
##Generate data
#Sequence of event times generated by a loglinear-NHPP for marks1.
#The observation window for marks1 starts at t=0 and ends at t=30.
t1 <- rloglinNHPP(Tstart=0, Tend=30, gamma0=-1.5, gamma1=0.10)
nr1 <- 1:length(t1) #cumulative recurrences
#Sequence of event times generated by a loglinear-NHPP for marks2.
#The observation window for marks2 also starts at t=0 and ends at t=30.
t2 <- rloglinNHPP(Tstart=0, Tend=30, gamma0=-3, gamma1=0.15)
nr2 <- 1:length(t2) #cumulative recurrences
#Sequence of event times generated by a loglinear-NHPP for marks3.
#Again, the observation window starts at t=0 and ends at t=30.
t3 <- rloglinNHPP(Tstart=0, Tend=30, gamma0=-2.1, gamma1=0.12)
nr3 <- 1:length(t3) #cumulative recurrences
#combine the generated data
tm <- data.frame(t=c(t1,t2,t3), nr=c(nr1,nr2,nr3),
marks= c(rep("marks1", length(t1)),
rep("marks2", length(t2)),
rep("marks3", length(t3))))
#plot the generated data
plot(tm$t, tm$nr, type="p", xlim=c(0, 30), pch=as.numeric(tm$marks),
col=as.numeric(tm$marks), xlab="Time", ylab="Cumulative recurrences")
legend("topleft", pch=as.numeric(unique(tm$marks)), col=as.numeric(unique(tm$marks)),
legend=as.character(unique(tm$marks)), bty="n", y.intersp=1.5)
##Fit a multitype temporal NHPP model
#First, check the order of the type labels:
#The vector elements of Tstart and Tend have to be ordered according to
#the order the type labels
levels(tm$marks)
#Thus, the first element of Tstart should contain the start time of the
#observation window for "marks1", the second element the start time for "marks2",
#and the third element the start time for "marks3".
#Similarly, the first element of Tend should contain the end time of the
#observation window for "marks1", the second element the end time for
#"marks2", and the third element the end time for "marks3".
#Fit an interaction model with a loglinear intensity function:
#The interaction model assumes that the three marks (=types) all have different
#values for the gamma0 and gamma1 parameters of the loglinear intensity function.
#Note that all models estimated by the multiNHPP function are fitted
#by using Maximum Likelihood Estimation (MLE).
ppgLoglin <- multiNHPP(Times=tm$t, Tstart=rep(0, 3), Tend=rep(30,3),
marks=tm$marks, model="interaction",
intensity="loglinear", nsegments=rep(6,3))
ppgLoglin[1:3]
#Note that the parameter values for gamma0 and gamma1 of "marks1" are labeled
#gamma0.marks1 and gamma1.marks1, respectively, in the output.
#The parameter values for the other two marks are labeled in a similar way.
##Variance covariance matrix for the model parameters
round(solve(-1*ppgLoglin$hessian), 6)
#Note that the zero's in the matrix reflect that the temporal points (=event times)
#of each of the three types (=marks) are generated by three
#independent temporal point processes
##Compute confidence intervals for the gamma0 and gamma1 parameter estimates
#95% confidence intervals (normal-approximation) for gamma0 and gamma1
#of the fitted loglinear-NHPP model
#compute standard errors
SEparms <- sqrt(diag(solve(-1*ppgLoglin$hessian)))
#marks1
ppgLoglin$parameters[1] + c(1,-1)*qnorm(.05/2)*SEparms[1] #95% ci for gamma0
ppgLoglin$parameters[4] + c(1,-1)*qnorm(.05/2)*SEparms[4] #95% ci for gamma1
#compare these with the population values for gamma0 and gamma1 of marks1
#remember, the population values for gamma0 and gamma1 were -1.5 and 0.1, respectively
#marks2
ppgLoglin$parameters[2] + c(1,-1)*qnorm(.05/2)*SEparms[2] #95% ci for gamma0
ppgLoglin$parameters[5] + c(1,-1)*qnorm(.05/2)*SEparms[5] #95% ci for gamma1
#population values for gamma0 and gamma1 of marks2 were -3 and 0.15, respectively
#marks3
ppgLoglin$parameters[3] + c(1,-1)*qnorm(.05/2)*SEparms[3] #95% ci for gamma0
ppgLoglin$parameters[6] + c(1,-1)*qnorm(.05/2)*SEparms[6] #95% ci for gamma1
#population values for gamma0 and gamma1 of marks3 were -2.1 and 0.12, respectively
##Diagnostics
#Use the fitted model to transform, for this "parallel systems configuration",
#the NHPP-times into corresponding Homogeneous Poisson Process (HPP) times,
#where the HPP is a unit-rate process (i.e., with intensity equal to 1).
hppInt <- tHPPparallel(ppgLoglin)
#For a unit-rate HPP, the number of cumulative recurrences is expected
#to be 1 at t=1, the number of cumulative recurrences at t=2 is expected
#to be 2, et cetera.
#In other words, when plotting the HPP-times against the cumulative recurrences
#for a unit-rate HPP, points are expected to lie close to the line y=x.
plot(hppInt$eventIndex, hppInt$hppTimes, type="l", col="blue",
xlab="Cumulative recurrences", ylab="HPP-time")
abline(a=0, b=1, col="gray", lty=2) #line y=x
#For a HPP, the interrecurrence times should be independent.
acf(diff(hppInt$hppTimes), main="Interrecurrence times HPP")
##Predictions: cumulative recurrences
#separate fits for the three types (i.e., "marks1", "marks2", and "marks3")
plot(tm$t, tm$nr, type="p", xlim=c(0, 30), pch=as.numeric(tm$marks),
col=as.numeric(tm$marks), xlab="Time", ylab="Cumulative recurrences")
parms <- ppgLoglin$parameters #extract parameters of fitted NHPP-model
tis <- seq(0, 30, length.out=100) #generate sequence of times
lines(tis, (exp(parms[1])*(exp(parms[4]*tis)-1))/parms[4], col="black") #marks1
lines(tis, (exp(parms[2])*(exp(parms[5]*tis)-1))/parms[5], col="red") #marks2
lines(tis, (exp(parms[3])*(exp(parms[6]*tis)-1))/parms[6], col="green") #marks3
legend("topleft", pch=as.numeric(unique(tm$marks)),
col=as.numeric(unique(tm$marks)), legend=as.character(unique(tm$marks)),
bty="n", y.intersp=1.5)
#Remember that the data from the three types may be seen as a
#"parallel systems configuration".
#As a consequence, it is possible to compute the cumulative recurrences
#of the three systems together by using the superposition principle.
sptimes <- sort(tm$t) #recurrence times of the three superimposed systems
#observed cumulative recurrences of the three superimposed systems
sprec <- 1:length(sptimes)
#generate a sequence of times within the observation window of the three systems
tis <- seq(0.1, 30, length.out=100)
#compute predicted cumulative recurrences (and their standard errors)
#for the three superimposed systems by using the superposition principle:
sp <- "(exp(gamma0.marks1)*(exp(gamma1.marks1*ts)-1))/gamma1.marks1+
(exp(gamma0.marks2)*(exp(gamma1.marks2*ts)-1))/gamma1.marks2+
(exp(gamma0.marks3)*(exp(gamma1.marks3*ts)-1))/gamma1.marks3"
SEs <- sapply(1:length(tis),
function (i) deltaMethod(ppgLoglin$parameters,
g=sp,
vcov.=solve(-1*ppgLoglin$hessian),
constants=c(ts=tis[i])))
#95% confidence interval for the predicted cumulative 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 recurrences
plot(sptimes, sprec, col=gray(level=.5), xlab="Time",
ylab="Cumulative recurrences", ylim=c(min(llci), max(luci)),
main="Superimposed systems")
#add predicted cumulative recurrences to the plot
lines(tis, unlist(SEs[1,]), col="red", lty=2, lwd=2)
#add the 95% confidence intervals
lines(tis, llci, col="blue", lty=2)
lines(tis, luci, col="blue", lty=2)
##Fit alternative models to the temporal data
#1) Fit a simple model
#The simple model assumes that the three types (=marks) all have different
#values for the gamma0 parameter of the loglinear intensity function. However,
#in contrast to the interaction model, the simple model assumes that the gamma1
#parameter is identical for all three types.
(modsimple <- multiNHPP(Times=tm$t, Tstart=rep(0, 3), Tend=rep(30,3),
marks=tm$marks, model="simple",
intensity="loglinear", nsegments=rep(6,3)))[1:3]
#diagnostics
hppSim <- tHPPparallel(modsimple)
plot(hppSim$eventIndex, hppSim$hppTimes, type="l", col="blue",
xlab="Cumulative recurrences", ylab="HPP-time")
abline(a=0, b=1, col="gray", lty=2) #line y=x
acf(diff(hppSim$hppTimes), main="Interrecurrence times HPP")
#2) Fit a homogenous model
#The homogeneous model assumes that the three types (=marks) all have different
#values for the gamma0 parameter of the loglinear intensity function.
#However, this model assumes that the intensity is constant (homogeneous) over
#time. In other words, the homogeneous model assumes that the gamma1
#parameter is equal to zero.
#The intensity for each type is thus defined as: exp(gamma0).
(modhom <- multiNHPP(Times=tm$t, Tstart=rep(0, 3), Tend=rep(30,3),
marks=tm$marks, model="homogeneous"))[1:3]
exp(modhom$parameters) #intensity estimate for the three types
#diagnostics
hppHom <- tHPPparallel(modhom)
#notice the bad fit of the homogeneous model
plot(hppHom$eventIndex, hppHom$hppTimes, type="l", col="blue",
xlab="Cumulative recurrences", ylab="HPP-time")
abline(a=0, b=1, col="gray", lty=2) #line y=x
acf(diff(hppHom$hppTimes), main="Interrecurrence times HPP")
#EXAMPLE 2: "two systems in series with a change point"
#Data for two different types of temporal points are generated.
#The two types of points are labeled "marks1" and "marks2".
#The intensity of both types is described by a loglinear
#function, with the intensity at time t defined by exp(gamma0 + gamma1*t).
#The observation window for "marks1" starts at t=0 and ends at t=25.
#The observation window for "marks2" starts at t=25 and ends at t=50.
#Note that this data also can be regarded as two systems operating in series.
#System 1 (labeled as "marks1") is the only system operating up to the
#change point at t=25. However, at the change point, System 2 (="marks2")
#takes over and from then on will be the only operating system.
##Generate data
#Sequence of event times generated by a loglinear-NHPP for marks1.
#The observation window for marks1 starts at t=0 and ends at t=25.
t1 <- rloglinNHPP(Tstart=0, Tend=25, gamma0=-1.5, gamma1=0.10)
#Sequence of event times generated by a loglinear-NHPP for marks2.
#The observation window for marks2 starts at t=25 and ends at t=50.
t2 <- rloglinNHPP(Tstart=25, Tend=50, gamma0=-6.5, gamma1=0.15)
#combine the generated data
tm <- data.frame(t=c(t1,t2),
marks= c(rep("marks1", length(t1)), rep("marks2", length(t2))))
#plot generated data
plot(tm$t, 1:nrow(tm), type="p", pch=as.numeric(tm$marks), xlim=c(0, 50),
xlab="Time", ylab="Cumulative recurrences")
legend("topleft", pch=as.numeric(unique(tm$marks)),
legend=as.character(unique(tm$marks)), bty="n", y.intersp=1.5)
##Fit a multitype temporal NHPP model
#Check the order of the type labels:
#The vector elements of Tstart and Tend have to be ordered according to
#the order the type labels
levels(tm$marks)
#Thus, the first element of Tstart should contain the start time of the
#observation window for "marks1", while the second element the start time
#for "marks2". Thus, Tstart=c(0, 25).
#Similarly, the first element of Tend should contain the end time of the
#observation window for "marks1", and the second element the end time for
#"marks2". Hence, Tend=c(25, 50).
#Fit an interaction model with a loglinear intensity function
ppgLoglin <- multiNHPP(Times=tm$t, Tstart=c(0, 25), Tend=c(25, 50),
marks=tm$marks, model="interaction",
intensity="loglinear", nsegments=rep(5,2))
ppgLoglin[1:3]
##Compute confidence intervals for the gamma0 and gamma1 parameter estimates
#95% confidence intervals (normal-approximation) for gamma0 and gamma1
#of the fitted loglinear-NHPP model
#compute standard errors
SEparms <- sqrt(diag(solve(-1*ppgLoglin$hessian)))
#system operating before the change point (labeled as "marks1")
ppgLoglin$parameters[1] + c(1,-1)*qnorm(.05/2)*SEparms[1] #95% ci for gamma0
ppgLoglin$parameters[3] + c(1,-1)*qnorm(.05/2)*SEparms[3] #95% ci for gamma1
#compare these with the population values for gamma0 and gamma1 of "marks1"
#remember, the population values for gamma0 and gamma1 were -1.5 and 0.1, respectively
#system operating after the change point (labeled as "marks2")
ppgLoglin$parameters[2] + c(1,-1)*qnorm(.05/2)*SEparms[2] #95% ci for gamma0
ppgLoglin$parameters[4] + c(1,-1)*qnorm(.05/2)*SEparms[4] #95% ci for gamma1
#population values for gamma0 and gamma1 of System 2 were -6.5 and 0.15, respectively
#Diagnostics
#Use the fitted model to transform, for this "series systems configuration",
#the NHPP-times into corresponding Homogeneous Poisson Process (HPP) times,
#where the HPP is a unit-rate process (i.e., with intensity equal to 1).
hppInt <- tHPPseries(ppgLoglin)
#Are the points close to the line y=x?
plot(hppInt$eventIndex, hppInt$hppTimes, type="l", col="blue",
xlab="Cumulative recurrences", ylab="HPP-time")
abline(a=0, b=1, col="gray", lty=2) #line y=x
#For a HPP, the interrecurrence times should be independent.
acf(diff(hppInt$hppTimes), main="Interrecurrence times HPP")
##Predictions: cumulative recurrences
#observed cumulative recurrences of the two systems together
recurrencesT <- 1:nrow(tm)
#generate sequence of times within the observation window of System 1
tis1 <- seq(.1, 25, length.out=50)
#generate sequence of times within the observation window of System 2
tis2 <- seq(25, 50, length.out=50)
#compute predicted cumulative recurrences (and their standard errors)
#System 1 ("marks1")
sp1 <- "(exp(gamma0.marks1)*(exp(gamma1.marks1*ts)-1))/gamma1.marks1"
SEs1 <- sapply(1:length(tis1),
function (i) deltaMethod(ppgLoglin$parameters,
g=sp1,
vcov.=solve(-1*ppgLoglin$hessian),
constants=c(ts=tis1[i])))
#95% confidence interval for the predicted cumulative recurrences
predrec1 <- unlist(SEs1[1,])
llci1 <- unlist(SEs1[1,])*exp((qnorm(.025)*unlist(SEs1[2,]))/unlist(SEs1[1,]))
luci1 <- unlist(SEs1[1,])*exp((qnorm(.975)*unlist(SEs1[2,]))/unlist(SEs1[1,]))
#System 2 ("marks2")
#Make sure that the cumulative recurrences at the ending time of System 1 and
#the starting time of System 2 (that is, at the change point) are identical.
sp2 <- "((exp(gamma0.marks1)*(exp(gamma1.marks1*tcp)-1))/gamma1.marks1-
(exp(gamma0.marks2)*(exp(gamma1.marks2*tcp)-1))/gamma1.marks2)+
(exp(gamma0.marks2)*(exp(gamma1.marks2*ts2)-1))/gamma1.marks2"
tcp <- tis1[length(tis1)] #position of change point
SEs2 <- sapply(1:length(tis2),
function (i) deltaMethod(ppgLoglin$parameters,
g=sp2,
vcov.=solve(-1*ppgLoglin$hessian),
constants=c(tcp=tcp,
ts2=tis2[i])))
#95% confidence interval for the predicted cumulative recurrences
predrec2 <- unlist(SEs2[1,])
llci2 <- unlist(SEs2[1,])*exp((qnorm(.025)*unlist(SEs2[2,]))/unlist(SEs2[1,]))
luci2 <- unlist(SEs2[1,])*exp((qnorm(.975)*unlist(SEs2[2,]))/unlist(SEs2[1,]))
#combine data
predrec <- c(predrec1, predrec2)
lci <- c(llci1, llci2)
uci <- c(luci1, luci2)
#plot the observed cumulative recurrences
plot(tm$t, recurrencesT, pch=as.numeric(tm$marks), col=gray(level=.2),
xlab="Time", ylab="Cumulative recurrences", ylim=c(min(lci), max(uci)))
#add predicted cumulative recurrences to the plot
lines(c(tis1,tis2), predrec, col="red", lty=2)
#add the 95% confindence intervals
lines(c(tis1,tis2), lci, col="blue", lty=2)
lines(c(tis1,tis2), uci, col="blue", lty=2)
#add change point
abline(v=25, lty=2, col="gray")
#EXAMPLE 3: "two systems in series with a change point"
#Similar to Example 2, but this time the intensity of both types
#(i.e., "marks1" and "marks2") is described by a power law intensity
#function, with the intensity at time t defined as gamma0*t^gamma1.
##Generate data
#Sequence of event times generated by a power-NHPP for marks1.
#The observation window for marks1 starts at t=0 and ends at t=30.
t1 <- rpowerNHPP(Tstart=0, Tend=30, gamma0=0.3, gamma1=0.7)
#Sequence of event times generated by a power-NHPP for marks2.
#The observation window for marks2 starts at t=30 and ends at t=60.
t2 <- rpowerNHPP(Tstart=30, Tend=60, gamma0=0.008, gamma1=1.4)
#combine data
tm <- data.frame(t=c(t1,t2),
marks= c(rep("marks1", length(t1)), rep("marks2", length(t2))))
#plot generated data
plot(tm$t, 1:nrow(tm), type="p", pch=as.numeric(tm$marks), xlim=c(0, 60),
xlab="Time", ylab="Cumulative recurrences")
legend("topleft", pch=as.numeric(unique(tm$marks)),
legend=as.character(unique(tm$marks)), bty="n", y.intersp=1.5)
##Fit a multitype temporal NHPP model
#Check the order of the type labels:
#The vector elements of Tstart and Tend have to be ordered according to
#the order the type labels
levels(tm$marks)
#Thus, Tstart=c(0, 30), and Tend=c(30, 60).
#Fit an interaction model with a power law intensity function
ppgPower <- multiNHPP(Times=tm$t, Tstart=c(0, 30), Tend=c(30, 60),
marks=tm$marks, model="interaction",
intensity="power", nsegments=rep(7,2))
ppgPower[1:3]
#Note that gamma0 in the power law model is modeled as log(gamma0).
#Thus, the value for gamma0 of "marks1" is actually exp(gamma0.marks1), and the
#the value for gamma0 of "marks2" is exp(gamma0.marks2).
##Compute confidence intervals for the gamma0 and gamma1 parameter estimates
#95% confidence intervals (normal-approximation) for gamma0 and gamma1
#of the fitted power law NHPP model
#compute standard errors
SEparms <- sqrt(diag(solve(-1*ppgPower$hessian)))
#system before the change point (labeled as "marks1")
exp(ppgPower$parameters[1] + c(1,-1)*qnorm(.05/2)*SEparms[1]) #95% ci for gamma0
ppgPower$parameters[3] + c(1,-1)*qnorm(.05/2)*SEparms[3] #95% ci for gamma1
#compare these with the population values for gamma0 and gamma1 of "marks1"
#remember, the population values for gamma0 and gamma1 were 0.3 and 0.7, respectively
#system after the change point (labeled as "marks2")
round(exp(ppgPower$parameters[2] + c(1,-1)*qnorm(.05/2)*SEparms[2]), 5) #95% ci for gamma0
ppgPower$parameters[4] + c(1,-1)*qnorm(.05/2)*SEparms[4] #95% ci for gamma1
#population values for gamma0 and gamma1 of System 2 were 0.008 and 1.4, respectively
##Diagnostics
#Use the fitted model to transform, for this "series systems configuration",
#the NHPP-times into corresponding Homogeneous Poisson Process (HPP) times,
#where the HPP is a unit-rate process (i.e., with intensity equal to 1).
hppInt <- tHPPseries(ppgPower)
#Are the points close to the line y=x?
plot(hppInt$eventIndex, hppInt$hppTimes, type="l", col="blue",
xlab="Cumulative recurrences", ylab="HPP-time")
abline(a=0, b=1, col="gray", lty=2) #line y=x
#For a HPP, the interrecurrence times should be independent.
acf(diff(hppInt$hppTimes), main="Interrecurrence times HPP")
##Predictions: cumulative recurrences
#observed cumulative recurrences of the two systems together
recurrencesT <- 1:nrow(tm)
#generate sequence of times within the observation window of System 1
tis1 <- seq(.1, 30, length.out=50)
#generate sequence of times within the observation window of System 2
tis2 <- seq(30, 60, length.out=50)
#compute predicted cumulative recurrences (and their standard errors)
#System 1 ("marks1")
sp1 <- "((exp(gamma0.marks1)/(gamma1.marks1+1))*ts^(gamma1.marks1+1))"
SEs1 <- sapply(1:length(tis1),
function (i) deltaMethod(ppgPower$parameters,
g=sp1,
vcov.=solve(-1*ppgPower$hessian),
constants=c(ts=tis1[i])))
#95% confidence interval for the predicted cumulative recurrences
predrec1 <- unlist(SEs1[1,])
llci1 <- unlist(SEs1[1,])*exp((qnorm(.025)*unlist(SEs1[2,]))/unlist(SEs1[1,]))
luci1 <- unlist(SEs1[1,])*exp((qnorm(.975)*unlist(SEs1[2,]))/unlist(SEs1[1,]))
#System 2 ("marks2")
#Make sure that the cumulative recurrences at the ending time of System 1 and
#the starting time of System 2 (that is, at the change point) are identical.
sp2 <- "(((exp(gamma0.marks1)/(gamma1.marks1+1))*tcp^(gamma1.marks1+1))-
((exp(gamma0.marks2)/(gamma1.marks2+1))*tcp^(gamma1.marks2+1)))+
((exp(gamma0.marks2)/(gamma1.marks2+1))*ts2^(gamma1.marks2+1))"
tcp <- tis1[length(tis1)]
SEs2 <- sapply(1:length(tis2),
function (i) deltaMethod(ppgPower$parameters,
g=sp2,
vcov.=solve(-1*ppgPower$hessian),
constants=c(tcp=tcp,
ts2=tis2[i])))
#95% confidence interval for the predicted cumulative recurrences
predrec2 <- unlist(SEs2[1,])
llci2 <- unlist(SEs2[1,])*exp((qnorm(.025)*unlist(SEs2[2,]))/unlist(SEs2[1,]))
luci2 <- unlist(SEs2[1,])*exp((qnorm(.975)*unlist(SEs2[2,]))/unlist(SEs2[1,]))
#combine data
predrec <- c(predrec1, predrec2)
lci <- c(llci1, llci2)
uci <- c(luci1, luci2)
#plot the observed cumulative recurrences
plot(tm$t, recurrencesT, pch=as.numeric(tm$marks), col=gray(level=.2),
xlab="Time", ylab="Cumulative recurrences", ylim=c(min(lci), max(uci)))
#add predicted cumulative recurrences to the plot
lines(c(tis1, tis2), predrec, col="red", lty=2)
#add the 95% confidence intervals
lines(c(tis1, tis2), lci, col="blue", lty=2)
lines(c(tis1, tis2), uci, col="blue", lty=2)
#add change point
abline(v=30, lty=2, col="gray")
```