rm(list=ls(all=TRUE))

# the working directory has to be correctly set
# otherwise R cannot load the next two files
source("mfreg.txt")
load("cutepi.rdata")

# to analyze new data, they must be correctly supplied:
# follow-up times must be stored in the vector tau
# total number of recurrences observed in each subject must be stored in the vector nr
# recurrence times must be stored in the matrix tempi; each row of this matrix
# refers to a subject; the number of columns equals the maximum number of recurrences
# (8 in our case). Therefore most of the cells of matrix tempi are possibly empty
# (filled with zeroes).




# Script 1
allpts <- cmfplot(tau,tempi,nr,out=TRUE,plot=FALSE)
plot(c(0,2180),c(0,0.3),type="n",xlab="Days",ylab="Cumulative Mean Number of Contacts")
x <- allpts[,1]
y <- allpts[,2]; lines(x,y,lty=1,lwd=2)
y <- allpts[,4]; lines(x,y,lty=3,lwd=1.5)
y <- allpts[,5]; lines(x,y,lty=3,lwd=1.5)




# Script 2
ok <- which(Xcov[,2]==0)
mal <- cmfplot(tau[ok],tempi[ok,],nr[ok],out=TRUE,plot=FALSE)
ok <- which(Xcov[,2]==1)
fem <- cmfplot(tau[ok],tempi[ok,],nr[ok],out=TRUE,plot=FALSE)
plot(c(0,2180),c(0,0.32),type="n",xlab="Days",ylab="Cumulative Mean Number of Tumours",axes=FALSE)
axis(2)
axis(1,at=c(0,seq(365,365*6,by=365)))
box()
lines(fem[,1],fem[,2],type="l",lwd=2)
lines(mal[,1],mal[,2],type="l",lty=3,lwd=2)
legend(0,0.30,c("females","males"),lty=c(1,2),lwd=2)




# Script 3
xcov <- matrix(Xcov[,1],ncol=1)
fit <- mfreg(tau,tempi,nr,xcov)
fit$estimates
xcov <- matrix(Xcov[,2],ncol=1)
fit <- mfreg(tau,tempi,nr,xcov)
fit$estimates
xcov <- Xcov[,1:2]
fit <- mfreg(tau,tempi,nr,xcov)
fit$estimates
xcov <- Xcov[,1:3]
fit <- mfreg(tau,tempi,nr,xcov)
fit$estimates




# Script 4
ok <- which(Xcov[,2]==0)
mal <- cmfplot(tau[ok],tempi[ok,],nr[ok],out=TRUE,plot=FALSE)
ok <- which(Xcov[,2]==1)
fem <- cmfplot(tau[ok],tempi[ok,],nr[ok],out=TRUE,plot=FALSE)
plot(c(0,2180),c(0,0.32),type="n",xlab="Days",ylab="Cumulative Mean Number of Tumours",axes=FALSE)
axis(2)
axis(1,at=c(0,seq(365,365*6,by=365)))
box()
lines(fem[,1],fem[,2],type="l",lwd=2)
lines(mal[,1],mal[,2],type="l",lty=3,lwd=2)
legend(0,0.30,c("females","males"),lty=c(1,2),lwd=2)
xcov <- matrix(Xcov[,2],ncol=1)
fit <- mfreg(tau,tempi,nr,xcov)
x <- fit$times; y <- cumsum(fit$basemf)
lines(x,y,lty=3)
x <- fit$times; y <- cumsum(fit$basemf*exp(fit$est[1,1]))
lines(x,y,lty=3)