library(survival)



# Example 1

n <- 10
t <- round(rexp(n,rate=0.2),2)
c <- round(rexp(n,rate=0.1),2)

d <- t < c
x <- ifelse(d,t,c)

ord <- order(x)
cat(paste(x[ord],ifelse(d,"","+")[ord],sep=""),sep=",") ;  cat("\n")

t.proc <- seq(10^(-6),20,length=1000)
alpha.t <-   rep(0.2,1000)

x.seg <- c(0,sort(x[d]))
n.seg <- length(x.seg)-1


y.proc      <- rep(NA,1000)
n.proc      <- rep(NA,1000)
lambda.proc <- rep(NA,1000)
Lambda.proc  <- rep(0,1000)

for(i in 1:1000)
{
        y.proc[i]      <- sum( t.proc[i] <= x)
        n.proc[i]      <- sum( x <= t.proc[i] & d)
        lambda.proc[i] <- y.proc[i] * alpha.t[i]
        if(i>1)
                Lambda.proc[i] <- Lambda.proc[i-1] + lambda.proc[i] * (t.proc[i]-t.proc[i-1])
       
}

ylm <- range(c(Lambda.proc,n.proc))
xlm <- c(0,max(x)+4)

par(mfcol=c(2,1))

plot(t.proc,Lambda.proc,type="l",xlim=xlm,ylim=ylm,xlab="Time",ylab="Number of obs. - intensity",
main="Event counting process and Intensity process")
segments(x.seg,0:n.seg,c(x.seg[-1],40),0:n.seg)

plot(t.proc,n.proc - Lambda.proc,type="l",xlim=xlm,xlab="Time",ylab="Error",main = "M process")





#____________________________________________________________________________
#
#	Input counting process simulation study - RUN FOLLOWING COMMANDS:
#____________________________________________________________________________


count.proc.sim <- 
function(n=50,plt=c(T,T,F,F))
{
#
t <- rexp(n,rate=0.2)
c <- rexp(n,rate=0.1)
#
d <- t < c
x <- ifelse(d,t,c)
#
t.proc <- seq(0.00001,20,length=1000)
alpha.t <-   rep(0.2,1000)
#
y.proc      <- rep(NA,1000)
n.proc      <- rep(NA,1000)
lambda.proc <- rep(NA,1000)
gamma.proc  <- rep(0,1000)

m.proc      <- rep(0,1000)
j.proc      <- rep(1,1000)
h.proc      <- rep(0,1000)
z.proc      <- rep(0,1000)
a.hat       <- rep(0,1000)
a.star      <- rep(0,1000)
var.proc    <- rep(0,1000)
#
#


for(i in 1:1000)
{
        y.proc[i]      <- sum( t.proc[i] <= x)
        n.proc[i]      <- sum( x <= t.proc[i] & d)
        lambda.proc[i] <- y.proc[i] * alpha.t[i]
        if(i>1)
                gamma.proc[i] <- gamma.proc[i-1] + lambda.proc[i] * (t.proc[i]-t.proc[i-1])

        m.proc[i] <- n.proc[i] - gamma.proc[i]

        if(y.proc[i] == 0)
                 j.proc[i] <- 0
        if(y.proc[i] >  0)
                 h.proc[i] <- 1 / y.proc[i]

        if(i>1)
        {
                a.star[i]   <- a.star[i-1] + j.proc[i] * alpha.t[i] * (t.proc[i]-t.proc[i-1])
                a.hat[i]    <- a.hat[i-1]  + h.proc[i] * (n.proc[i]-n.proc[i-1])
                z.proc[i]   <- z.proc[i-1] + h.proc[i] * (m.proc[i]-m.proc[i-1])
                var.proc[i] <- var.proc[i-1] + y.proc[i]^(-2) * (n.proc[i]-n.proc[i-1])
        }       
}
#
#
num.plt <- sum(plt)
if(num.plt>0)
{
par(mfcol=c(1,num.plt))
if(plt[1])
plot(survfit(Surv(x,d),type="kaplan-meier",se.fit=F))
#
if(plt[2])
{
ylm <- range(c(gamma.proc,n.proc))
plot(t.proc,gamma.proc,type="l",ylim=ylm,xlab="Time",ylab="Number of obs. - intensity")
lines(t.proc,n.proc)
}
#
if(plt[3])
{
plot(t.proc,m.proc,type="l",xlab="Time",ylab="Martingale process")
abline(0,0)
}
#
if(plt[4])
{
plot(sqrt(var.proc),z.proc,type="l",xlab="Time",
ylab="Z process")
abline(0,0)
}
}
#
#
return(t.proc,gamma.proc,n.proc,m.proc,a.star,a.hat,z.proc,var.proc)
}



#  Example 2


a <- count.proc.sim(plt=c(T,T,T,F),n=20)



#		Plot 8 replications of N and Gamma processes

n.mat <- array(dim=c(1000,8))
gamma.mat <- array(dim=c(1000,8))
a.star.mat <- array(dim=c(1000,8))
a.hat.mat <- array(dim=c(1000,8))
z.mat <- array(dim=c(1000,8))
var.mat <- array(dim=c(1000,8))

for(i in 1:8)
{
a <- count.proc.sim(plt=c(F,F,F,F),n=50)
n.mat[,i]      <- a$n.proc
gamma.mat[,i]  <- a$gamma.proc
a.star.mat[,i] <- a$a.star 
a.hat.mat[,i]  <- a$a.hat 
z.mat[,i]      <- a$z.proc
var.mat[,i]    <- a$var.proc
}

t.proc <- a$t.proc
n <- 50



#	Plot 8 replications of N-A cumulative hasard estimate

par(mfcol=c(1,1))
a.star <- apply(a.star.mat,1,max)
ylm <- range(a.hat.mat)
plot(t.proc,a.hat.mat[,1],ylim=ylm,xlab="Time",ylab="Predicted cummulative risk",main="N-A estimate",type="l")
for(i in 2:8) lines(t.proc,a.hat.mat[,i])
lines(t.proc,a.star,col=2)



#  1000 replcations


NA.est.10 <- rep(NA,1000)
NA.var.10 <- rep(NA,1000)

for(i in 1:1000)
{
a <- count.proc.sim(plt=c(F,F,F,F),n=800)
NA.est.10[i]  <- (a$a.hat[500] + a$a.hat[501])/2
NA.var.10[i]  <- (a$var.proc[500] + a$var.proc[501])/2
}

par(mfcol=c(1,2))
hist(NA.est.10,main=paste("NA-estimates - sd =",round(sqrt(var(NA.est.10)),2)))
boxplot(sqrt(NA.var.10),main="sqrt -- of variance estimates")