msmdata<-function(mydata,ncov,ncovdata,nstate){
	msmdata<-matrix(data=NA,nrow=nrow(mydata)+max(mydata[,1]),ncol=(4+ncovdata))
	
	i<-0
	q<-2
	msmdata[1,1]<-1
	msmdata[1,2]<-0
	msmdata[1,3]<-1
	msmdata[1,4]<-mydata[1,5]
	msmdata[2,1]<-1
	msmdata[2,2]<-mydata[1,3]
	msmdata[2,4]<-mydata[1,5]
	if (ncov!=0) {for (k in 1:ncovdata) {	msmdata[1,4+k]<-mydata[1,5+k]
			   		msmdata[2,4+k]<-mydata[1,5+k]}}
				     		    
	ifelse (mydata[1,4]==1,msmdata[2,3]<-nstate,ifelse(mydata[1,3]==mydata[2,2],msmdata[2,3]<-mydata[2,5]+1,msmdata[2,3]<-1))
	if (mydata[2,1]==1) {	msmdata[3,1]<-1
				msmdata[3,2]<-mydata[2,3]
				ifelse (mydata[2,4]==1,msmdata[3,3]<-nstate,msmdata[3,3]<-mydata[2,5]+1)
				if (ncov!=0) {for (j in 1:(ncovdata+1)) {msmdata[3,3+j]<-mydata[2,4+j]} }
				q<-q+1}	

	for (k in q:(nrow(mydata)-1)) { 	msmdata[k+i+1,1]<-mydata[k,1]   #ptnum
					if (ncov!=0) {for (j in 1:(ncovdata+1)) {msmdata[k+i+1,3+j]<-mydata[k,4+j]} } #covariaveis
					msmdata[k+i+1,2]<-mydata[k,2]	#dtime
					if (msmdata[k+i+1,2]==0) msmdata[k+i+1,3]<-1	#state
	 				if (msmdata[k+i+1,2]!=0 & mydata[k,1]==mydata[k-1,1]) msmdata[k+i+1,3]<-mydata[k,5]+1	#state
					if (mydata[k,1]!=mydata[k+1,1]) {	#linha extra
						i<-i+1
						msmdata[k+i+1,1]<-mydata[k,1]
						if (ncov!=0) {for (j in 1:(ncovdata+1)) {msmdata[k+i+1,3+j]<-mydata[k,4+j]} }
			   			msmdata[k+i+1,2]<-mydata[k,3]
 						ifelse (mydata[k,4]==1,msmdata[k+i+1,3]<-nstate,msmdata[k+i+1,3]<-msmdata[k+i,3])}}

	msmdata<-msmdata[1:(nrow(msmdata)-1),]
	return(msmdata)
}		

nelsonomp<-function(msmdata,nstate){
nmsm<-nrow(msmdata)
times<-sort(unique(msmdata[,2]))
time<-data.frame(times)
ntimes<-dim(time)[1]+1
nomp<-matrix(data=NA,nrow=ntimes,ncol=(4*nstate-3))
for (k in 1:(ntimes-1)) nomp[k,1]<-times[k]
y<-vector(length=(nstate-1))
for (k in 1:(nstate-1)) y[k]<-0
aux<-vector(length=(nstate-1))
aux2<-vector(length=(nstate-1))
for (k in (3*nstate-1):(4*nstate-4)) nomp[1,k]<-0
#primeira linha da matriz
for (k in 1:(4*nstate-4)) nomp[1,k]<-0
nomp[1,(2*nstate-1)]<-max(msmdata[,1])
nomp[1,(3*nstate-2)]<-1
nomp[1,(4*nstate-3)]<-1
#ultima linha
nomp[ntimes,1]<-nomp[ntimes-1,1]+1

for (k in 2:(nstate-1)) nomp[2,2*nstate-2+k]<-0
nomp[2,2*nstate-1]<-max(msmdata[,1])

for (j in 2:(ntimes-1)) {
	for (t in 1:(nstate-1)) {	aux[t]<-0
				aux2[t]<-0}
	for (l in 1:(nmsm-1)) {	 if (msmdata[l,2]==nomp[j,1]) { 	for (k in 2:nstate) { if (msmdata[l-1,3]==1 & msmdata[l,3]==k) { 	aux[k-1]<-aux[k-1]+1
																y[1]<-y[1]+1
																y[k]<-y[k]+1}}}}
								
	for (k in 2:nstate) nomp[j,k]<-aux[k-1]
	
	for (l in 1:(nmsm-1)) { if (msmdata[l,2]==nomp[j,1]) {for (k in 2:(nstate-1)) { if (msmdata[l-1,3]==k & msmdata[l,3]==nstate) { 	aux2[k-1]<-aux2[k-1]+1
																y[k]<-y[k]-1}}}}
	for (k in 2:(nstate-1)) nomp[j,k+nstate-1]<-aux2[k-1]

	for (k in 2:(nstate-1)) nomp[j+1,2*nstate-2+k]<-y[k]
	nomp[j+1,2*nstate-1]<-max(msmdata[,1])-y[1]
	
	#definir probabilidades p110t
	sum<-sum(nomp[j,2:nstate])
	
	if (nomp[j,(2*nstate-1)]==0) nomp[j,(3*nstate-2)]<-0
	else nomp[j,(3*nstate-2)]<-nomp[j-1,(3*nstate-2)]*(1-sum/nomp[j,(2*nstate-1)])  #p110t    

	#definir probabilidades p1j0t com j=2,...,nstate-1
	for (k in 2:(nstate-1)) {	nomp[1,k+3*nstate-3]<-0
				if (nomp[j,k+2*nstate-2]!=0 & nomp[j,2*nstate-1]!=0) nomp[j,k+3*nstate-3]<-nomp[j-1,k+3*nstate-3]*(1-nomp[j,k+nstate-1]/nomp[j,k+2*nstate-2])+nomp[j-1,3*nstate-2]*nomp[j,k]/nomp[j,2*nstate-1]
				if (nomp[j,k+2*nstate-2]==0 & nomp[j,2*nstate-1]==0) nomp[j,k+3*nstate-3]<-nomp[j-1,k+3*nstate-3]
				if (nomp[j,k+2*nstate-2]==0) nomp[j,k+3*nstate-3]<-nomp[j-1,3*nstate-2]*nomp[j,k]/nomp[j,2*nstate-1]+nomp[j-1,k+3*nstate-3]
				if (nomp[j,2*nstate-1]==0) nomp[j,k+3*nstate-3]<-nomp[j-1,k+3*nstate-3]*(1-nomp[j,k+nstate-1]/nomp[j,k+2*nstate-2])} 	
			
	#definir sobrevivência
	nomp[j,4*nstate-3]<-sum(nomp[j,(3*nstate-2):(4*nstate-4)])}

	for (k in 2:(4*nstate-3)) nomp[ntimes,k]<-nomp[ntimes-1,k]

sobrevi<-matrix(ncol=(nstate+2),nrow=(2*nrow(nomp)-1))
for (k in 1:(nrow(nomp)-1)) { 	sobrevi[2*k-1,1]<-nomp[k,1]
			  	sobrevi[2*k,1]<-nomp[k+1,1]-0.001
				for (j in (3*nstate-2):(4*nstate-4)) { sobrevi[2*k-1,j-3*nstate+5]<-nomp[k,j]
								 sobrevi[2*k,j-3*nstate+5]<-nomp[k,j]}
			  	sobrevi[2*k-1,2]<-nomp[k,ncol(nomp)]
			  	sobrevi[2*k,2]<-nomp[k,ncol(nomp)]
				sobrevi[2*k-1,nstate+2]<-1-sum(nomp[k,(3*nstate-2):(4*nstate-4)])
			  	sobrevi[2*k,nstate+2]<-1-sum(nomp[k,(3*nstate-2):(4*nstate-4)])}
sobrevi[2*nrow(nomp)-1,1]<-nomp[nrow(nomp),1]
for (j in (3*nstate-2):(4*nstate-4)) { sobrevi[2*nrow(nomp)-1,j-3*nstate+5]<-nomp[nrow(nomp),j]}
sobrevi[2*nrow(nomp)-1,nstate+2]<-1-sum(nomp[nrow(nomp),(3*nstate-2):(4*nstate-4)])
sobrevi[2*nrow(nomp)-1,2]<-nomp[nrow(nomp),ncol(nomp)]
#plot(sobrevi[,1],sobrevi[,2],type="l",xlab="time",ylab="survival probability")
return(sobrevi)
}

#Auxiliary function for the calculation of the new matrix
newmatrix<-function(q,dim,dim1,nhaz,haz)
{
dim2<-dim+dim1
time<-vector(length=dim2)
for (k in 1:dim1) time[k]<-nhaz[k,1]
for (k in (dim1+1):dim2) time[k]<-haz[k-dim1,1]
times<-sort(unique(time))
dimtimes<-dim.data.frame(times)[2]
mymat<-matrix(data=0,ncol=q,nrow=dimtimes)
mymat[,1]<-times
for (k in 1:dimtimes) {
	for (r in 2:(q-1)) {	for (j in 1:dim1) {if (nhaz[j,1]==mymat[k,1]) mymat[k,r]<-nhaz[j,r] }}
	for (j in 1:(dim2-dim1)) {if (haz[j,1]==mymat[k,1]) mymat[k,q]<-haz[j,2]}}
nhaz<-mymat
return(nhaz)
}

# Function to calculate survival probability in CMM or CSMM
survival<-function(nhaz,nstate) 
{
haz<-matrix(ncol=ncol(nhaz),nrow=(nrow(nhaz)+1))
for (k in 1:ncol(nhaz)) {	haz[1,k]<-0
			for (j in 1:nrow(nhaz)) haz[j+1,k]<-nhaz[j,k]}
dim<-dim(haz)[1]
col<-dim(haz)[2]
pmatrix<-matrix(data=NA,ncol=(3*nstate-3),nrow=dim)
pmatrix[,1]<-haz[,1]
#p11(0,t|z) is column 2
pmatrix[1,2]<-1
for (j in 2:dim) {  aux<-1
		for (k in 2:nstate) { pmatrix[j,2]<- aux-haz[j,k]
				   aux<-pmatrix[j,2]}
				   pmatrix[j,2]<-pmatrix[j-1,2]*aux}
#p12(0,t|z) for column 3,..., p1a(0,t|z) for column a+1
for (j in 2:(nstate-1)) { 	pmatrix[1,j+1]<-0
			for (k in 2: dim) {	pmatrix[k,j+1]<-pmatrix[k-1,2]*haz[k,j]+pmatrix[k-1,j+1]*(1-haz[k,nstate+j-1])}}
#survival for column nstate+1
pmatrix[1,nstate+1]<-1
for (k in 2:dim) {	aux<-0
		for (j in 2:nstate) { pmatrix[k,nstate+1]<-aux+pmatrix[k,j]
				  aux<-pmatrix[k,nstate+1]}}
#pjj for columns from nstate+2 to 2nstate-1
for (k in (nstate+2):(2*nstate-1)) pmatrix[1,k]<-1
for (k in 2:dim) { for (j in (nstate+2):(2*nstate-1)) pmatrix[k,j]<-pmatrix[k-1,j]*(1-haz[k,j-1])}
#pjN for columns from 2nstate to 3nstate-3
for (k in (2*nstate):(3*nstate-3)) pmatrix[1,k]<-0
for (k in 2:dim) { for (j in (2*nstate):(3*nstate-3)) pmatrix[k,j]<-pmatrix[k-1,j]+pmatrix[k-1,j-nstate+2]*haz[k,j-nstate+1]}

sobrevi<-matrix(ncol=2*nstate,nrow=(2*nrow(pmatrix)-1))
for (k in 1:(nrow(pmatrix)-1)) { 	sobrevi[2*k-1,1]<-pmatrix[k,1]
			  	sobrevi[2*k,1]<-pmatrix[k+1,1]-0.001
			  	sobrevi[2*k-1,2]<-pmatrix[k,nstate+1]
			  	sobrevi[2*k,2]<-pmatrix[k,nstate+1]
				for (j in 2:nstate) {	sobrevi[2*k-1,j+1]<-pmatrix[k,j]
			  				sobrevi[2*k,j+1]<-pmatrix[k,j]}
				sobrevi[2*k-1,nstate+2]<-1-sum(sobrevi[2*k-1,3:(nstate+1)])
				sobrevi[2*k,nstate+2]<-1-sum(sobrevi[2*k-1,3:(nstate+1)])
				for (j in (nstate+3):(2*nstate)) { 	sobrevi[2*k-1,j]<-pmatrix[k,j+nstate-3]
								sobrevi[2*k,j]<-pmatrix[k,j+nstate-3]}}
sobrevi[2*nrow(pmatrix)-1,1]<-pmatrix[nrow(pmatrix),1]
sobrevi[2*nrow(pmatrix)-1,2]<-pmatrix[nrow(pmatrix),nstate+1]
for (j in 2:nstate) sobrevi[2*nrow(pmatrix)-1,j+1]<-pmatrix[nrow(pmatrix),j]
sobrevi[2*nrow(pmatrix)-1,nstate+2]<-1-sum(sobrevi[2*nrow(pmatrix)-1,3:(nstate+1)])
for (j in (nstate+3):(2*nstate)) sobrevi[2*nrow(pmatrix)-1,j]<-pmatrix[nrow(pmatrix),j+nstate-3]
return(sobrevi)
}

# Function to obtain the crudeinits for the multi-state homogeneous/non-homogeneous model
inits.tdc.surv<-function (state, time, subject, qmatrix, data = NULL, fromto = FALSE, 
    fromstate = NULL, tostate = NULL, timelag = NULL, check = FALSE) 
{
    if (!is.null(data)) data <- as.data.frame(data)
    if (fromto) {
        if (missing(fromstate) || missing(tostate) || missing(timelag)) 
            stop("fromstate, tostate and timelag must be specified if fromto is TRUE")
        state <- eval(substitute(fromstate), data, parent.frame())
        tostate <- eval(substitute(tostate), data, parent.frame())
        timelag <- eval(substitute(timelag), data, parent.frame())
        tottime <- tapply(timelag, state, sum)
        ntrans <- table(state, tostate)
    }
    else {
        if (missing(state) || missing(time) || missing(subject)) 
            stop("state, time, and subject must be specified if fromto is FALSE")
        state <- eval(substitute(state), data, parent.frame())
        n <- length(state)
        subject <- if (missing(subject)) 
            rep(1, n)
        else eval(substitute(subject), data, parent.frame())
        nextsubj <- c(subject[2:n], NA)
        lastsubj <- (subject != nextsubj)
        time <- eval(substitute(time), data, parent.frame())
        timecontrib <- ifelse(lastsubj, NA, c(time[2:n], 0) - time)
        tottime <- tapply(timecontrib[!lastsubj], state[!lastsubj], sum)
        ntrans <- statetable.msm(state, subject, data = NULL)
    }
    qmatrix <- msm.check.qmatrix(qmatrix)
    msm.check.state(nrow(qmatrix), state, fromto, tostate)
    if (check) 
        msm.check.model(state, subject, qmatrix, fromto, tostate)
    nst <- nrow(qmatrix)
    estmat <- matrix(0, nst, nst)
    rownames(estmat) <- colnames(estmat) <- paste(1:nst)
    tab <- sweep(ntrans, 1, tottime, "/")
    for (i in 1:nst) for (j in 1:nst) if ((paste(i) %in% rownames(tab)) && 
        (paste(j) %in% colnames(tab))) 
        estmat[paste(i), paste(j)] <- tab[paste(i), paste(j)]
    inits <- t(estmat)[t(qmatrix) == 1]
    names(inits) <- paste(t(row(qmatrix)), t(col(qmatrix)), sep = "-")[t(qmatrix) == 1]
    inits
}

# Function for graphical output for covariates and/or survival probabilities
plot.tdc.surv<-function(stanf,mydata,ncov,graphcov,model,newcov,covariate,formula,surv.plot,sobrevi,now,nstate,transition,plot.trans){

auxiliar<-0
if (graphcov==99) { 	graphcov<-1
			auxiliar<-1}
if (ncov!=0) {	rg1 <- range(mydata[,formula[graphcov]])
		covdiff <- (rg1[2] - rg1[1])/40
		cove <- seq(rg1[1], rg1[2], covdiff)
		legend.pos1 <- c(max(cove) - 15 * covdiff, 1)
		legend.pos5<-(rg1[2]+rg1[1])/2-abs(rg1[1]+(rg1[2]+rg1[1])/2)/4}

		rg2 <- range(mydata[,3])
    		timediff <- (rg2[2] - rg2[1])/40
    		times <- seq(rg2[1], rg2[2], timediff)
		legend.pos2 <- c(max(times) - 15 * timediff, 1)
		legend.pos3 <- max(times)/3
		legend.pos4 <- max(times)/2

if (auxiliar==1) graphcov<-0

if (model[1]==1) 	{ 	
	if (graphcov!=0) {		stanf1<-coxph(Surv(mydata$start,mydata$stop,mydata$event)~pspline(covariate[,graphcov]),data=mydata,na.action=na.exclude)
		  	 	tem<-predict(stanf1,type="terms",se.fit=T)
		   		tmat<-cbind(tem$fit,tem$fit-1.96*tem$se.fit,tem$fit+1.96*tem$se.fit)
 		   		jj<-match(sort(unique(covariate[,graphcov])),covariate[,graphcov])

cat(“ ”,"\n")	
cat(“Checking linearity assumption for ”, colnames(mydata)[formula[graphcov]] ,"\n")
print(stanf1)
cat(“ ”,"\n")

matplot(covariate[,graphcov][jj],tmat[jj,],type="l",lty=c(1,2,2),xlab="covariate",ylab="log hazard",main=list(paste("Smooth log-hazard (along with 95%CI) for",main=names(covariate)[graphcov])))
		   		legend(legend.pos1[1],1,legend="TDCM")}
	if (surv.plot==TRUE) {	#stanf<-coxph(Surv(mydata$start,mydata$stop,mydata$event)~covariate,data=mydata,na.action=na.exclude)
		   		plot(survfit(stanf),main="Cox regression estimates of survival",xlab="Survival time",ylab="Proportion Surviving",cex = 2, lty = 1:2,ylim = c(0, 1))
		  		legend(legend.pos4,legend.pos2[2],legend="TDCM")}}

if (model[2]==1) { 
	if (graphcov!=0) { 

cat(“ ”,"\n")	
cat(“Checking linearity assumption for ”, colnames(mydata)[formula[graphcov]] ,"\n")
print(stanf)
cat(“ ”,"\n")

		   	tem<-predict(stanf,type="terms",se.fit=T)
		   	tmat<-cbind(tem$fit,tem$fit-1.96*tem$se.fit,tem$fit+1.96*tem$se.fit)
 		   	jj<-match(sort(unique(newcov)),newcov)
		   	matplot(newcov[jj],tmat[jj,],type="l",lty=c(1,2,2),ylab="log hazard",main=list(paste("Smooth log-hazard (along with 95%CI) for",names(covariate)[graphcov],"(CMM)")))
		   	legend(legend.pos5,2,legend=list(paste("From stage",transition[1],"to stage",transition[2])))}
	if (surv.plot==TRUE & now==1) {	plot(sobrevi[,1],sobrevi[,2],type="l",main="Survival estimate (CMM)" ,xlab="Survival time",ylab="Survival probability")}
#	if (plot.trans==TRUE & now==1) { for (k in 3:(nstate+2)) {	plot(sobrevi[,1],sobrevi[,k],type="l",main="Estimates of transition probability (CMM)" ,xlab="Survival time")
#									legend(legend.pos3,0.22,legend=list(paste("From stage 1 to stage",k-2)))	 }}
		   }
if (model[3]==1) 	{  if (plot.trans==TRUE)	{ 	plot(stanf)
		   				title(main="Probability estimates for the mortality intensity (HMM)")
		   				#legend(0.7,0.3,legend="HMM")
					for (j in 1:(nstate-2)) { 
						timediff2<-(rg2[2]-rg2[1])/100
						time<-seq(rg2[1],rg2[2],timediff2)
						p1<-vector(length=101)
						for (k in 1:101) p1[k]<-pmatrix.msm(stanf,time[k])[1,j+1]
						plot(time,p1,type="l")
						title(main=list(paste("Transition Probability estimate from stage 1 to stage ",j+1,"(HMM)")))}
						#legend(0.7,0.22,legend=list(paste("From stage 1 to stage",j+1)))}
					}
		 if (surv.plot==TRUE)	{	timediff2<-(rg2[2]-rg2[1])/100
						time<-seq(rg2[1],rg2[2],timediff2)
						surv<-vector(length=101)
						for (k in 1:101) surv[k]<-1-pmatrix.msm(stanf,time[k])[1,nstate]
						plot(time,surv,type="l",main="Survival estimate (HMM)" ,xlab="Survival time",ylab="Survival probability")}}

if (model[4]==1 & plot.trans==TRUE) {  	plot(stanf)
		   			title(main="Probability estimates for the mortality intensity (NHM, first branch)") }
					#legend(0,0.5,legend="NHM")}
  		   			#legend(legend.pos2[1],legend.pos2[2],legend="NHM")}
if (model[5]==1) {  
	if (graphcov!=0) { 
cat(“ ”,"\n")	
cat(“Checking linearity assumption for ”, colnames(mydata)[formula[graphcov]] ,"\n")
print(stanf)
cat(“ ”,"\n")

			tem<-predict(stanf,type="terms",se.fit=T)
		   	tmat<-cbind(tem$fit,tem$fit-1.96*tem$se.fit,tem$fit+1.96*tem$se.fit)
 		   	jj<-match(sort(unique(newcov)),newcov)
		   	matplot(newcov[jj],tmat[jj,],type="l",lty=c(1,2,2),ylab="log hazard",main=list(paste("Smooth log-hazard (along with 95%CI) for",names(covariate)[graphcov],"(CSMM)")))
		   	legend(legend.pos5,2,legend=list(paste("From stage",transition[1],"to stage",transition[2])))}

	if (surv.plot==TRUE & now==TRUE) {	plot(sobrevi[,1],sobrevi[,2],type="l",main="Survival estimate (CSMM)" ,xlab="Survival time",ylab="Survival probability")}
#	if (plot.trans==TRUE & now==TRUE) { for (k in 3:(nstate+2)) {	plot(sobrevi[,1],sobrevi[,k],type="l",main="Estimates of transition probability (CSMM)" ,xlab="Survival time")
#									legend(legend.pos3,0.22,legend=list(paste("From stage 1 to stage",k-2)))	 }
#					        for (k in 2:(nstate-1)) {	plot(sobrevi[,1],sobrevi[,k+nstate+1],type="l",main="Estimates of probability transition (CSMM)",xlab="Survival time",ylab="Transition prob")
#									legend(legend.pos3,0.22,legend=list(paste("From stage", k, "to stage ", nstate)))	}}

		            }

if (model[6]==1)  {	if (surv.plot==TRUE) 	plot(sobrevi[,1],sobrevi[,2],type="l",main="Survival estimate (NPM)" ,xlab="Survival time",ylab="Survival probability")
					
			if (plot.trans==TRUE) {	for (k in 1:nstate) { 	plot(sobrevi[,1],sobrevi[,k+2],type="l",main=list(paste("Estimates of probability transition from state 1 to state ",k,"(NPM)")) ,xlab="Survival time",ylab="Transition probability")
									legend(legend.pos3,0.22,legend=list(paste("From stage 1 to stage ", legend=k)))}}
		}
}

# FUNCTION FOR CONTRASTES USING WALDS TEST
wald.tdc.surv<-function(stanf,nstate,ncov,fixed){

	#Wald test to check for differences between mortality intensities
	if (nstate==3){
		cat("***** Estimate of the ratio of the progression rate 1-3 *****","\n")  	
		cat("***** into   death   to   the  corresponding  rate  2-3 *****","\n")  	
		cat("Estimate:",qratio.msm(stanf,c(1,nstate),c(2,nstate))$estimate,"\n")		
		cat("SE:",qratio.msm(stanf,c(1,nstate),c(2,nstate))$se,"\n")			
		cat("\n")}

	V<-matrix(data=NA,nrow=(nstate-2),ncol=(nstate-2))
	for (j in 1:(nstate-2)) {for (k in 1:(nstate-2)) {V[j,k]<-stanf$covmat[nstate-2+j,nstate-2+k]-stanf$covmat[nstate-2+k,2*nstate-3]-stanf$covmat[nstate-2+j,2*nstate-3]+stanf$covmat[2*nstate-3,2*nstate-3]}}  
	R<-vector(length=(nstate-2))
	for (j in 1:(nstate-2)) {R[j]<-coef.msm(stanf)$logbaseline[j,nstate]-coef.msm(stanf)$logbaseline[nstate-1,nstate]}
	W<-t(R)%*%(solve(V))%*%R
	p<-1-pchisq(W,df=(nstate-2))	
	cat("\n")
	cat("**** WALDS TEST TO CHECK FOR DIFFERENCES BETWEEN MORTALITY INTENSITIES **** ","\n")							
	cat("LET H0 BE THE HYPOTHESIS OF NON-DIFFERENCES BETWEEN MORTALITY INTENSITIES","\n")	
	cat("H0 produces a chisquare statistic of:",W,"\n")					
	cat("with a p-value of:",p,"\n")	
	cat("\n")

	#Wald test to check for differences in each covariate on the transition intensities
	if (fixed==0){
	if (ncov<=4 & ncov!=0) {
	ntrans<-2*nstate-3
	Vcov<-matrix(data=NA,nrow=(ntrans-1),ncol=(ntrans-1))
	for (l in 1:ncov){
	for (j in 1:(ntrans-1)) {for (k in 1:(ntrans-1)) {Vcov[j,k]<-stanf$covmat[ntrans*l+1+j-1,ntrans*l+1+k-1]-stanf$covmat[ntrans*l+1+j-1,ntrans*(l+1)]-stanf$covmat[ntrans*l+1+k-1,ntrans*(l+1)]+stanf$covmat[ntrans*(l+1),ntrans*(l+1)]}}  #ERRO
	Rcov<-vector(length=(ntrans-1))

	x<-data.frame(coef(stanf)[l+1])
	for (j in 1:(nstate-2)) {Rcov[j]<-x[1,j+1]-x[nstate-1,nstate]}
	for (j in (nstate-2+1):(ntrans-1)) {Rcov[j]<-x[j-1,nstate]-x[nstate-1,nstate]}
	Wcov<-t(Rcov)%*%(solve(Vcov))%*%Rcov
	pcov<-1-pchisq(Wcov,df=(nstate-1))	
	cat("\n")
	cat("**** WALDS TEST TO CHECK FOR DIFFERENCES BETWEEN COVARIATE ",l," IN THE TRANSITION INTENSITIES **** ","\n")							
	cat("LET H0 BE THE HYPOTHESIS OF NON-DIFFERENCES BETWEEN COVARIATE ",l," IN THE TRANSITION INTENSITIES","\n")	
	cat("H0 produces a chi-square statistic of:",Wcov,"\n")					
	cat("with a p-value of:",pcov,"\n")	
	cat("\n")}}}

}

#FUNCTION TO FIT TIME-DEPENDENT COX REGRESSION MODELS
tdcrm.tdc.surv<-function(mydata,ncov,formula,covariate,newmodel,covmat,graphcov,surv.plot,plot.trans){
		   newmodel[1]<-1
		   now<-0
		   covariate<-matrix(data=NA,nrow=nrow(mydata),ncol=ncov)
		   covariate<-data.frame(covariate)
		   if (ncov==0) covariate<-c(rep(1,nrow(mydata)))
		   if (ncov!=0) { 	covariate[,1]<-mydata[,formula[1]]
				names(covariate)[1]<-c(names(mydata)[formula[1]])}
		   if (ncov>1) { for (k in 2:ncov) { 	covariate[,k]<-mydata[,formula[k]]
						names(covariate)[k]<-c(names(mydata)[formula[k]])}}
		   if (ncov!=0) {
			covu<-vector(length=ncov)
			covu<-names(covariate)
 			covu2<-c(covu,"treat")
			covar <- as.formula(paste(" Surv(mydata$start,mydata$stop,mydata$event)~ ", paste(covu, collapse= "+")))#}
			covar2 <- as.formula(paste(" Surv(mydata$start,mydata$stop,mydata$event)~ ", paste(covu2, collapse= "+")))}

 		   ifelse(ncov==0, stanf<-coxph(Surv(mydata$start,mydata$stop,mydata$event)~1,data=mydata),stanf<-coxph(covar,data=mydata))   
		   ifelse(ncov==0, stanf2<-coxph(Surv(mydata$start,mydata$stop,mydata$event)~1,data=mydata),stanf2<-coxph(covar2,data=mydata))   

		   cat("*****************************************************************************","\n")
		   cat("    ***** TIME-DEPENDENT COX REGRESSION MODEL *****    ","\n")
		   cat("*****************************************************************************","\n")
		   if (ncov!=0) {	print(summary(stanf2))}
		   cat("-2*Log-likelihood:",-2*stanf2$loglik[2],"\n")
	               if (ncov!=0 & covmat==TRUE) { cat("*** variance-covariance matrix ***","\n")
			       	         		print(stanf2$var)
				         		cat("\n") }
		   if (graphcov!=99 | surv.plot==TRUE) {	#stanf<-coxph(Surv(mydata$start,mydata$stop,mydata$event)~pspline(covariate[,graphcov]),data=mydata)
							#cat("*** P-spline fit for covariate 1 *** ","\n")
							#print(stanf)
							plot.tdc.surv(stanf,mydata,ncov,graphcov,newmodel,newcov,covariate,formula,surv.plot,sobrevi,now,nstate,transition,plot.trans)}
		   cat("\n")
		   newmodel[1]<-0}

#FUNCTION TO FIT COX MARKOV (SEMI-MARKOV) MODELS
cmm.tdc.surv<-function(mydata,msmdata,nstate,ncov,covariate,formula,graphcov,newmodel,covmat,model,writedata,coxdata,surv.plot,plot.trans){

lines<-nrow(mydata)-1+(nstate-2)*sum(mydata[,5]==0)
	l<-nstate
	p<-2
	coxdata<-matrix(data=NA,ncol=(ncol(mydata)+1),nrow=lines)
	for (k in 1:ncol(mydata)) {coxdata[1,k]<-mydata[1,k]}
	coxdata[1,ncol(mydata)+1]<-1

	for (j in 2:(nstate-1)) { 	coxdata[j,]<-coxdata[1,]
				coxdata[j,ncol(mydata)+1]<-j
				ifelse(mydata[2,5]==(j-1),coxdata[j,4]<-1,coxdata[j,4]<-0)}

	if (mydata[1,1]==mydata[2,1]) { 	for (k in 1:ncol(mydata)){coxdata[nstate,k]<-mydata[2,k]}
				            coxdata[nstate,ncol(mydata)+1]<-nstate+mydata[2,5]-1
					l<-l+1
					p<-p+1}

	for (k in p:(nrow(mydata)-1)) { 	if (mydata[k,1]!=mydata[k-1,1]) { 	for (r in 1:ncol(mydata)) {coxdata[l,r]<-mydata[k,r]}
									coxdata[l,ncol(mydata)+1]<-1
									l<-l+1
									for (j in 2:(nstate-1)) { 	coxdata[l+j-2,]<-coxdata[l-1,]
											    	coxdata[l+j-2,ncol(mydata)+1]<-j	
												ifelse(mydata[k+1,5]==(j-1) & k!=(nrow(mydata)-1),coxdata[l+j-2,4]<-1,coxdata[l+j-2,4]<-0)}	
									l<-l+nstate-2}					
													
				       	else { for (r in 1:ncol(mydata)){coxdata[l,r]<-mydata[k,r]}																			
				       	coxdata[l,ncol(mydata)+1]<-nstate+mydata[k,5]-1
				       	l<-l+1}}
						
	mydata2<-mydata[1:(nrow(mydata)-1),1:ncol(mydata)]

	coxdata<-data.frame(coxdata)
	names(coxdata)<-c(names(mydata),"trans")

	covariate<-matrix(data=NA,nrow=nrow(coxdata),ncol=ncov)
	covariate<-data.frame(covariate)
	if (ncov==0) covariate<-c(rep(1,nrow(coxdata)))
	
	if (ncov!=0) { 	covariate[,1]<-coxdata[,formula[1]]
			names(covariate)[1]<-c(names(mydata)[formula[1]])}
	if (ncov==1) {	covariate<-matrix(data=NA,nrow=nrow(coxdata),ncol=2)
			covariate<-data.frame(covariate)
			covariate[,1]<-coxdata[,formula[1]]
			covariate[,2]<-coxdata[,formula[1]]
			names(covariate)<-c(names(mydata)[formula[1]],names(mydata)[formula[1]])}

	if (ncov>1) { for (k in 2:ncov) { 	covariate[,k]<-coxdata[,formula[k]]
					names(covariate)[k]<-c(names(mydata)[formula[k]])}}
	if (ncov!=0) {	covu<-vector(length=ncov)
			covu<-names(covariate)
			covar <- as.formula(paste(" Surv(coxdata[,2],coxdata[,3],coxdata[,4])~ ", paste(covu, collapse= "+")))}

if (model[2]==1){
	newmodel[2]<-1
	t<-0
	q<-2
	now<-0			
	cat("*****************************************************************************","\n")
	cat("***** COX MARKOV MODEL. SEPARATED FITTED COX MODELS FOR EACH TRANSITION *****","\n")
	cat("*****************************************************************************","\n")	
	cat("\n")
	for (k in 2:(nstate-1)){cat("\n")
	t<-t+1
	cat("***** FROM STATE 1 TO STATE",k," *****","\n")
	ifelse(ncov==0, fit<-coxph(Surv(coxdata[,2],coxdata[,3],coxdata[,4])~1,data=coxdata,subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude),fit<-coxph(covar,data=coxdata,subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude))
	transition<-c(1,k)
	if (surv.plot==TRUE | plot.trans==TRUE) {
	if (k==2) {	y<-basehaz(fit)
			dim<-dim(y)[1]
			haz<-matrix(data=NA,ncol=2,nrow=dim)
			haz[1,2]<-y[1,1]
			haz[1,1]<-y[1,2]
			for (s in 2:dim) {	haz[s,2]<-y[s,1]-y[s-1,1]
		 			haz[s,1]<-y[s,2]}
			nhaz<-haz
			dim1<-dim}

	if (k>2) {	y<-basehaz(fit)
			dim<-dim(y)[1]
			haz<-matrix(data=NA,ncol=2,nrow=dim)
			haz[1,2]<-y[1,1]
			haz[1,1]<-y[1,2]
			for (s in 2:dim) {	haz[s,2]<-y[s,1]-y[s-1,1]
		 			haz[s,1]<-y[s,2]}
			q<-q+1
			nhaz<-newmatrix(q,dim,dim1,nhaz,haz)
			dim1<-dim(nhaz)[1]}}

	if (ncov!=0) { print(summary(fit))}
	cat("-2*Log-likelihood:",-2*fit$loglik[2],"\n")
	if (ncov!=0) { 
	if (covmat==TRUE) { cat("*** variance-covariance matrix *** ","\n")
			      print(fit$var)
			      cat("\n") }
	if (graphcov!=99) {	if (ncov==1) stanf<-coxph(Surv(coxdata[,2],coxdata[,3],coxdata[,4])~pspline(covariate[,1]),subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude)
				if (ncov!=1)  stanf<-coxph(Surv(coxdata[,2],coxdata[,3],coxdata[,4])~pspline(covariate[,graphcov]),subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude)
				#cat("*** P-spline fit for covariate 1 *** ","\n")
				#print(stanf)
				newcov<-vector(length=(stanf$n))
				p<-1
				for (j in 1:nrow(coxdata)) {if (coxdata[j,ncol(coxdata)]==t) { if (ncov==1) newcov[p]<-covariate[j,1]
											if (ncov!=1)  newcov[p]<-covariate[j,graphcov]
									   		p<-p+1}}
				plot.tdc.surv(stanf,mydata,ncov,graphcov,newmodel,newcov,covariate,formula,surv.plot,sobrevi,now,nstate,transition,plot.trans)}}
	cat("\n")}

	t<-t+1
	cat("***** FROM STATE 1 TO STATE",nstate," *****","\n")
	ifelse(ncov==0, fit<-coxph(Surv(coxdata[,2],coxdata[,3],coxdata[,4])~1,data=coxdata,subset=(coxdata[,ncol(coxdata)]==1),na.action=na.exclude),fit<-coxph(covar,data=coxdata,subset=(coxdata[,ncol(coxdata)]==1),na.action=na.exclude))
	transition<-c(1,nstate)
	if (surv.plot==TRUE | plot.trans==TRUE) {
	y<-basehaz(fit)
	dim<-dim(y)[1]
	haz<-matrix(data=NA,ncol=2,nrow=dim)
	haz[1,2]<-y[1,1]
	haz[1,1]<-y[1,2]
	for (s in 2:dim) {	haz[s,2]<-y[s,1]-y[s-1,1]
			haz[s,1]<-y[s,2]}

	q<-q+1
	nhaz<-newmatrix(q,dim,dim1,nhaz,haz)
	dim1<-dim(nhaz)[1]}

	if (ncov!=0) { print(summary(fit))}
	cat("-2*Log-likelihood:",-2*fit$loglik[2],"\n")
	if (ncov!=0) { 
	if (covmat==TRUE) { cat("*** variance-covariance matrix *** ","\n")
			      print(fit$var)
			      cat("\n") }
 	if (graphcov!=99) {	if (ncov==1) stanf<-coxph(Surv(coxdata[,2],coxdata[,3],coxdata[,4])~pspline(covariate[,1]),subset=(coxdata[,ncol(coxdata)]==1),na.action=na.exclude)
				if (ncov!=1)   stanf<-coxph(Surv(coxdata[,2],coxdata[,3],coxdata[,4])~pspline(covariate[,graphcov]),subset=(coxdata[,ncol(coxdata)]==1),na.action=na.exclude)
				#cat("*** P-spline fit for covariate 1 *** ","\n")
				#print(stanf)
				newcov<-vector(length=(stanf$n))
				p<-1
				for (j in 1:nrow(coxdata)) {if (coxdata[j,ncol(coxdata)]==t) {	if (ncov==1) newcov[p]<-covariate[j,1]
											if (ncov!=1)   newcov[p]<-covariate[j,graphcov]
									   		p<-p+1}}
				plot.tdc.surv(stanf,mydata,ncov,graphcov,newmodel,newcov,covariate,formula,surv.plot,sobrevi,now,nstate,transition,plot.trans)}}
	cat("\n")

	for (k in nstate:(2*nstate-3)){
	cat("\n")
	t<-t+1
	cat("***** FROM STATE",2+k-nstate," TO STATE",nstate," *****","\n")

	cat("Checking the Markov assumption:", "\n")
cat("testing if the time spent in state 1 (start) is important on transition from state", 2+k-nstate, " to state",nstate,"\n")
	TMA<-coxph(Surv(stop,event)~start,data=coxdata, subset=(coxdata[,ncol(coxdata)]==k))
print(TMA)
	cat("\n")

	ifelse(ncov==0, fit<-coxph(Surv(coxdata[,2],coxdata[,3],coxdata[,4])~1,data=coxdata,subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude),fit<-coxph(covar,data=coxdata,subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude))
	transition<-c(2+k-nstate,nstate)

	if (surv.plot==TRUE | plot.trans==TRUE) {
	y<-basehaz(fit)
	dim<-dim(y)[1]						
	haz<-matrix(data=NA,ncol=2,nrow=dim)
	haz[1,2]<-y[1,1]
	haz[1,1]<-y[1,2]
	for (s in 2:dim) {	haz[s,2]<-y[s,1]-y[s-1,1]
			haz[s,1]<-y[s,2]}
	q<-q+1
	nhaz<-newmatrix(q,dim,dim1,nhaz,haz)
	dim1<-dim(nhaz)[1]
	if (k==(2*nstate-3)) {	sobrevi<-survival(nhaz,nstate) 
				now<-1}
	}

	if (ncov!=0) { print(summary(fit))}
	cat("-2*Log-likelihood:",-2*fit$loglik[2],"\n")
	if (ncov!=0) { 
	if (covmat==TRUE) { cat("*** variance-covariance matrix *** ","\n")
			      print(fit$var)
			      cat("\n")}}

	if (graphcov!=99 | surv.plot==TRUE) { if (ncov!=0)	{	if (ncov==1) stanf<-coxph(Surv(coxdata[,2],coxdata[,3],coxdata[,4])~pspline(covariate[,1]),subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude)
						if (ncov!=1 & graphcov!=99)   stanf<-coxph(Surv(coxdata[,2],coxdata[,3],coxdata[,4])~pspline(covariate[,graphcov]),subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude)
						if (ncov!=1 & graphcov==99) stanf<-coxph(Surv(coxdata[,2],coxdata[,3],coxdata[,4])~1,subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude)
						#cat("*** P-spline fit for covariate 1 *** ","\n")
						#print(stanf)
						newcov<-vector(length=(stanf$n))
						p<-1
						for (j in 1:nrow(coxdata)) {if (coxdata[j,ncol(coxdata)]==t) {	if (ncov==1) newcov[p]<-covariate[j,1]
											if (ncov!=1 & graphcov!=99)   newcov[p]<-covariate[j,graphcov]
									   		p<-p+1}}}
						plot.tdc.surv(stanf,mydata,ncov,graphcov,newmodel,newcov,covariate,formula,surv.plot,sobrevi,now,nstate,transition,plot.trans)}
	now<-0
	cat("\n")}

	newmodel[2]<-0
	if (writedata==TRUE) { 	cat("*** FORMATED DATA *** ","\n")
			       	write.table(coxdata,file="C:\COXDATA.txt",sep=",",na = "NA",dec = ".",col.names = TRUE)}}

if (model[5]==1){
	newmodel[5]<-1
	t<-0
	q<-2
	now<-0			
	cat("******************************************************************************","\n")
	cat("*** COX SEMI-MARKOV MODEL. SEPARATED FITTED COX MODELS FOR EACH TRANSITION ***","\n")
	cat("******************************************************************************","\n")	
	cat("\n")

	if (ncov!=0) {	covu<-vector(length=ncov)
			covu<-names(covariate)
			covar <- as.formula(paste(" Surv(coxdata[,3]-coxdata[,2],coxdata[,4])~ ", paste(covu, collapse= "+")))}

	for (k in 2:(nstate-1)){cat("\n")
	t<-t+1
	cat("***** FROM STATE 1 TO STATE",k," *****","\n")
	ifelse(ncov==0, fit<-coxph(Surv(coxdata[,3]-coxdata[,2],coxdata[,4])~1,data=coxdata,subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude),fit<-coxph(covar,data=coxdata,subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude))
	transition<-c(1,k)
	if (surv.plot==TRUE | plot.trans==TRUE) {
	if (k==2) {	y<-basehaz(fit)
			dim<-dim(y)[1]
			haz<-matrix(data=NA,ncol=2,nrow=dim)
			haz[1,2]<-y[1,1]
			haz[1,1]<-y[1,2]
			for (s in 2:dim) {	haz[s,2]<-y[s,1]-y[s-1,1]
		 			haz[s,1]<-y[s,2]}
			nhaz<-haz
			dim1<-dim}

	if (k>2) {	y<-basehaz(fit)
			dim<-dim(y)[1]
			haz<-matrix(data=NA,ncol=2,nrow=dim)
			haz[1,2]<-y[1,1]
			haz[1,1]<-y[1,2]
			for (s in 2:dim) {	haz[s,2]<-y[s,1]-y[s-1,1]
		 			haz[s,1]<-y[s,2]}
			q<-q+1
			nhaz<-newmatrix(q,dim,dim1,nhaz,haz)
			dim1<-dim(nhaz)[1]}}

	if (ncov!=0) { print(summary(fit))
		      cat("-2*Log-likelihood:",-2*fit$loglik[2],"\n")}
	if (ncov==0) cat("-2*Log-likelihood:",-2*fit$loglik[1],"\n")
	if (ncov!=0) { 
	if (covmat==TRUE) { cat("*** variance-covariance matrix *** ","\n")
			      print(fit$var)
			      cat("\n") }
	if (graphcov!=99) {	if (ncov==1) stanf<-coxph(Surv(coxdata[,3]-coxdata[,2],coxdata[,4])~pspline(covariate[,1]),subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude)
				if (ncov!=1 & graphcov!=99)   stanf<-coxph(Surv(coxdata[,3]-coxdata[,2],coxdata[,4])~pspline(covariate[,graphcov]),subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude)
				#cat("*** P-spline fit for covariate 1 *** ","\n")
				#print(stanf)
				newcov<-vector(length=(stanf$n))
				p<-1
				for (j in 1:nrow(coxdata)) {if (coxdata[j,ncol(coxdata)]==t) {	if (ncov==1) newcov[p]<-covariate[j,1]
											if (ncov!=1)   newcov[p]<-covariate[j,graphcov]
									   		p<-p+1}}
				plot.tdc.surv(stanf,mydata,ncov,graphcov,newmodel,newcov,covariate,formula,surv.plot,sobrevi,now,nstate,transition,plot.trans)}}
	cat("\n")}

	t<-t+1
	cat("***** FROM STATE 1 TO STATE",nstate," *****","\n")
	ifelse(ncov==0, fit<-coxph(Surv(coxdata[,3]-coxdata[,2],coxdata[,4])~1,data=coxdata,subset=(coxdata[,ncol(coxdata)]==1),na.action=na.exclude),fit<-coxph(covar,data=coxdata,subset=(coxdata[,ncol(coxdata)]==1),na.action=na.exclude))
	transition<-c(1,nstate)
	if (surv.plot==TRUE | plot.trans==TRUE) {
	y<-basehaz(fit)
	dim<-dim(y)[1]
	haz<-matrix(data=NA,ncol=2,nrow=dim)
	haz[1,2]<-y[1,1]
	haz[1,1]<-y[1,2]
	for (s in 2:dim) {	haz[s,2]<-y[s,1]-y[s-1,1]
			haz[s,1]<-y[s,2]}

	q<-q+1
	nhaz<-newmatrix(q,dim,dim1,nhaz,haz)
	dim1<-dim(nhaz)[1]}

	if (ncov!=0) { print(summary(fit))
		      cat("-2*Log-likelihood:",-2*fit$loglik[2],"\n")}
	if (ncov==0) cat("-2*Log-likelihood:",-2*fit$loglik[1],"\n")

	if (ncov!=0) { 
	if (covmat==TRUE) { cat("*** variance-covariance matrix *** ","\n")
			      print(fit$var)
			      cat("\n") }
 	if (graphcov!=99) {	if (ncov==1) stanf<-coxph(Surv(coxdata[,3]-coxdata[,2],coxdata[,4])~pspline(covariate[,1]),subset=(coxdata[,ncol(coxdata)]==1),na.action=na.exclude)
				if (ncov!=1)   stanf<-coxph(Surv(coxdata[,3]-coxdata[,2],coxdata[,4])~pspline(covariate[,graphcov]),subset=(coxdata[,ncol(coxdata)]==1),na.action=na.exclude)
				#cat("*** P-spline fit for covariate 1 *** ","\n")
				#print(stanf)
				newcov<-vector(length=(stanf$n))
				p<-1
				for (j in 1:nrow(coxdata)) {if (coxdata[j,ncol(coxdata)]==t) {	if (ncov==1) newcov[p]<-covariate[j,1]
											if (ncov!=1 & graphcov!=99)   newcov[p]<-covariate[j,graphcov]
									   		p<-p+1}}
				plot.tdc.surv(stanf,mydata,ncov,graphcov,newmodel,newcov,covariate,formula,surv.plot,sobrevi,now,nstate,transition,plot.trans)}}
	cat("\n")

	for (k in nstate:(2*nstate-3)){cat("\n")
	t<-t+1
	cat("***** FROM STATE",2+k-nstate," TO STATE",nstate," *****","\n")

cat("Checking the Markov assumption:", "\n")
cat("testing if the time spent in state 1 (start) is important on transition from state", 2+k-nstate, " to state",nstate,"\n")
	TMA<-coxph(Surv(stop,event)~start,data=coxdata, subset=(coxdata[,ncol(coxdata)]==k))
print(TMA)
	cat("\n")

	ifelse(ncov==0, fit<-coxph(Surv(coxdata[,3]-coxdata[,2],coxdata[,4])~1,data=coxdata,subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude),fit<-coxph(covar,data=coxdata,subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude))
	transition<-c(2+k-nstate,nstate)
	if (surv.plot==TRUE | plot.trans==TRUE) {
	y<-basehaz(fit)
	dim<-dim(y)[1]
	haz<-matrix(data=NA,ncol=2,nrow=dim)
	haz[1,2]<-y[1,1]
	haz[1,1]<-y[1,2]
	for (s in 2:dim) {	haz[s,2]<-y[s,1]-y[s-1,1]
			haz[s,1]<-y[s,2]}

	q<-q+1
	nhaz<-newmatrix(q,dim,dim1,nhaz,haz)
	dim1<-dim(nhaz)[1]
	if (k==(2*nstate-3)) {	sobrevi<-survival(nhaz,nstate) 
				now<-1}
	}
#	now<-1

	if (ncov!=0) { print(summary(fit))
		      cat("-2*Log-likelihood:",-2*fit$loglik[2],"\n")}
	if (ncov==0) cat("-2*Log-likelihood:",-2*fit$loglik[1],"\n")

	if (ncov!=0) { 
	if (covmat==TRUE) { cat("*** variance-covariance matrix *** ","\n")
			      print(fit$var)
			      cat("\n") }}
	if (graphcov!=99 | surv.plot==TRUE) { if (ncov!=0) {		if (ncov==1) stanf<-coxph(Surv(coxdata[,3]-coxdata[,2],coxdata[,4])~pspline(covariate[,1]),subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude)
				if (ncov!=1 & graphcov!=99)   stanf<-coxph(Surv(coxdata[,3]-coxdata[,2],coxdata[,4])~pspline(covariate[,graphcov]),subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude)
				if (ncov!=1 & graphcov==99) stanf<-coxph(Surv(coxdata[,2],coxdata[,3],coxdata[,4])~1,subset=(coxdata[,ncol(coxdata)]==k),na.action=na.exclude)
				#cat("*** P-spline fit for covariate 1 *** ","\n")
				#print(stanf)
				newcov<-vector(length=(stanf$n))
				p<-1
				for (j in 1:nrow(coxdata)) {if (coxdata[j,ncol(coxdata)]==t) { if (ncov==1) newcov[p]<-covariate[j,1]
											if (ncov!=1 & graphcov!=99)   newcov[p]<-covariate[j,graphcov]
									   		p<-p+1}}}
				plot.tdc.surv(stanf,mydata,ncov,graphcov,newmodel,newcov,covariate,formula,surv.plot,sobrevi,now,nstate,transition,plot.trans)}
	cat("\n")}
	newmodel[5]<-0
}
}

#FUNCTION TO FIT HOMOGENEOUS MARKOV MODEL OR NON-HOMOGENEOUS PIECEWISE MODEL
multistate<-function(mydata,msmdata,ncov,ncovdata,nstate,formula,covmat,prevalences,writedata,graphcov,model,newmodel,fixedpars,fixed,cut,exacttimes,surv.plot,plot.trans,timeunit){
					     
#	msmdata<-msmdata(mydata,ncov,ncovdata,nstate)
	msmdata<-data.frame(msmdata)
	names(msmdata)<-c("ptnum","dtime","state",names(mydata)[5:ncol(mydata)])
	covariate<-matrix(data=NA,nrow=nrow(msmdata),ncol=ncov)
	covariate<-data.frame(covariate)

		if (ncov!=0) { 	covariate[,1]<-msmdata[,formula[1]-1]
				names(covariate)[1]<-c(names(msmdata)[formula[1]-1])}
		if (ncov>1) { for (k in 2:ncov) { 	covariate[,k]<-msmdata[,formula[k]-1]
						names(covariate)[k]<-c(names(msmdata)[formula[k]-1])}}
		   if (ncov!=0) {
			covu<-vector(length=ncov)
			covu<-names(covariate)
			covar <- as.formula(paste("~", paste(covu, collapse= "+")))}

	qmat<-matrix(data=NA,nrow=nstate,ncol=nstate)
	for (k in 2:(nstate-1)) { 	qmat[1,]<-c(0,rep(1,nstate-1))
				qmat[k,]<-c(rep(0,nstate-1),1)
				qmat[nstate,]<-c(rep(0,nstate))}

	init<-inits.tdc.surv(msmdata[,3], msmdata[,2], msmdata[,1],data=msmdata,qmatrix=qmat) 
	inits<-c(init,rep(0,(2*nstate-3)*ncov))
	
# MULTI-STATE HOMOGENEOUS MARKOV MODEL	
if (model[3]==1){
	newmodel[3]<-1
	cat("*****************************************************************************","\n")
	cat("   ***** MULTI-STATE HOMOGENEOUS MARKOV MODEL *****    ","\n")
	cat("*****************************************************************************","\n")
	cat("\n")
	cat("***** Undergoing Transitions *****","\n")
	print(statetable.msm(msmdata[,3],msmdata[,1],data=msmdata))
	cat("\n")
	cat("*** convergence ***","\n")

	if (ncov==0) { 	if (exacttimes==TRUE) stanf<-msm(msmdata[,3]~msmdata[,2],subject=msmdata[,1],qmatrix=qmat,inits=init,death=TRUE,method="BFGS",control=list(trace=2),exacttimes=T)
			else stanf<-msm(msmdata[,3]~msmdata[,2],subject=msmdata[,1],qmatrix=qmat,inits=init,death=TRUE,method="BFGS",control=list(trace=2),exacttimes=F,tunit=timeunit)}
	if (ncov!=0) { 	if (exacttimes==TRUE) stanf<-msm(msmdata[,3]~msmdata[,2],subject=msmdata[,1],qmatrix=qmat,data=msmdata,covariates=covar,inits=inits,death=TRUE,method="BFGS",control=list(trace=2),exacttimes=T,fixedpars=fixedpars)
			else stanf<-msm(msmdata[,3]~msmdata[,2],subject=msmdata[,1],qmatrix=qmat,data=msmdata,covariates=covar,inits=inits,death=TRUE,method="BFGS",control=list(trace=2),exacttimes=F,fixedpars=fixedpars,tunit=timeunit)}
	
	cat("\n")
	cat("-2*Log-likelihood:",stanf$minus2loglik,"\n")
	cat("\n")
	cat("***** Estimated coefficients *****","\n")
	print(coef.msm(stanf))
	print(hazard.msm(stanf))
	
	wald.tdc.surv(stanf,nstate,ncov,fixed)			
	cat("\n")
	cat("***** Sojourn times *****","\n")					
	sojourn.msm(stanf)
	print(sojourn.msm(stanf))
	cat("\n")	
		
	if (covmat==TRUE) { cat("*** variance-covariance matrix *** ","\n")
			      print(stanf$covmat)
			      cat("\n") }
	if (prevalences==TRUE) { cat("*** PREVALENCES *** ","\n")
			         	print(prevalence.msm(stanf,times=seq(0,max(msmdata[,2]),length=10))[3:4])}
	if (writedata==TRUE) { 	cat(" *** FORMATED DATA *** ","\n")
			       	write.table(msmdata,file="c:\HMM.txt",sep=",",na = "NA",dec = ".",col.names = colnames(msmdata))}
	if (surv.plot==TRUE | plot.trans==TRUE) plot.tdc.surv(stanf,mydata,ncov,graphcov,newmodel,newcov,covar,formula,surv.plot,sobrevi,now,nstate,transition,plot.trans)
	newmodel[3]<-0}
		
# MULTI-STATE NON-HOMOGENEOUS PIECEWISE MODEL
	if (model[4]==1) { if (missing(cut)) stop("Missing cut")
			  newmodel[4]<-1
			  msmdata2<-rbind(msmdata,c(rep(0,ncol(msmdata))))
			  q<-0
			  r<-0
			  msmcutone<-msmdata2[1,]
			  for (k in 2:(nrow(msmdata2)-1)) { x<-msmdata2[k,]

						       if (q==1) { 	y<-msmdata2[k-1,]
								y[1]<-max(msmdata2[,1])+r+1
								x[1]<-y[1]
								if (msmdata2[k,2]-msmdata2[k-1,2]>=cut) {x[2]<-msmdata2[k-1,2]+cut
													  x[3]<-msmdata2[k-1,3]}
								msmcutone<-rbind(msmcutone,y,x)
								r<-r+1
								q<-0}
								 
						       else {
						             	if (msmdata2[k,1]==msmdata2[k-1,1]) {
			       		               		if (msmdata2[k,2]-msmdata2[k-1,2]>cut & msmdata2[k-1,3]==1) { 	x[3:ncol(msmdata2)]<-msmdata2[k-1,3:ncol(msmdata2)] 
															  	x[2]<-cut 
															   	if (msmdata2[k,1]==msmdata2[k+1,1]) q<-1 }
															    
						       		if (msmdata2[k,2]-msmdata2[k-1,2]>cut & msmdata2[k-1,3]!=1) { 	x[3:ncol(msmdata2)]<-msmdata2[k-1,3:ncol(msmdata2)] 
															   	x[2]<-msmdata2[k-1,2]+cut }}
							     	msmcutone<-rbind(msmcutone,x) }}	
			  
			  cat("*****************************************************************************","\n")
			  cat("           ***** MULTI-STATE NON-HOMOGENEOUS MODEL *****              ","\n")
	                          cat("*****************************************************************************","\n")
			  cat("\n")
			  cat("PRIMEIRO RAMO","\n")
			  msmcutone<-data.frame(msmcutone)
			  covariate<-matrix(data=NA,nrow=nrow(msmcutone),ncol=ncov)

			covariate<-data.frame(covariate)

			if (ncov!=0) { 	covariate[,1]<-msmcutone[,formula[1]-1]
					names(covariate)[1]<-c(names(msmdata)[formula[1]-1])}
			if (ncov>1) { for (k in 2:ncov) { 	covariate[,k]<-msmcutone[,formula[k]-1]
							names(covariate)[k]<-c(names(msmdata)[formula[k]-1])}}
		   	if (ncov!=0) {
					covu<-vector(length=ncov)
					covu<-names(covariate)
					covar <- as.formula(paste("~", paste(covu, collapse= "+")))}	  
			  cat("\n")
			  cat("*** convergence ***","\n")	                       
			  if (ncov==0) {	if (exacttimes==TRUE) stanf<-msm(msmcutone[,3]~msmcutone[,2],subject=msmcutone[,1],qmatrix=qmat,inits=init,death=TRUE,method="BFGS",control=list(trace=2),exacttimes=T)
					else stanf<-msm(msmcutone[,3]~msmcutone[,2],subject=msmcutone[,1],qmatrix=qmat,inits=init,death=TRUE,method="BFGS",control=list(trace=2),exacttimes=F,tunit=timeunit)}
			  if (ncov!=0) { if (exacttimes==TRUE) stanf<-msm(msmcutone[,3]~msmcutone[,2],subject=msmcutone[,1],qmatrix=qmat,data=msmcutone,covariates=covar,inits=inits, death=TRUE, method="BFGS", control=list(trace=2),exacttimes=T,fixedpars=fixedpars)
					else stanf<-msm(msmcutone[,3]~msmcutone[,2],subject=msmcutone[,1],qmatrix=qmat,data=msmcutone, covariates=covar, inits=inits, death=TRUE,method="BFGS", control=list(trace=2),exacttimes=F,fixedpars=fixedpars,tunit=timeunit)}
			  
			  cat("\n")
			  cat("-2*Log-likelihood:",stanf$minus2loglik,"\n")
			  cat("\n")
			  cat("***** Estimated coefficients *****","\n")
			  print(coef.msm(stanf))
		              print(hazard.msm(stanf))
			  cat("\n")
			  
			  wald.tdc.surv(stanf,nstate,ncov,fixed)			
			  cat("\n")

			  cat("***** Sojourn times *****","\n")					
			  sojourn.msm(stanf)
			  print(sojourn.msm(stanf))
			  cat("\n")					

			  if (covmat==TRUE) { cat("*** variance-covariance matrix *** ","\n")
			       		        print(stanf$covmat)
					        cat("\n") }
			  if (prevalences==TRUE) { cat("*** PREVALENCES *** ","\n")
			        		               print(prevalence.msm(stanf,times=seq(min(msmcutone[,2]),max(msmcutone[,2]),length=10))[3:4])}
			  cat("\n")
			  if (surv.plot==TRUE | plot.trans==TRUE ) plot.tdc.surv(stanf,mydata,ncov,graphcov,newmodel,newcov,covariate,formula,surv.plot,sobrevi,now,nstate,transition,plot.trans)
			  cat("\n")

			  # SEGUNDO RAMO

			   msmdata3<-rbind(msmdata2,c(rep(0,ncol(msmdata2))))
			   y<-msmdata3[1,]
			   for (k in 1:(nrow(msmdata3)-2)) { 
				if (msmdata3[k,2]==0){
				    if (msmdata3[k+1,2]-msmdata3[k,2]>cut & msmdata3[k,1]!=msmdata3[k+2,1]) { 
								x<-msmdata3[k,]
								x[2]<-cut
								y<-rbind(y,x)
								x<-msmdata3[k+1,]
								y<-rbind(y,x)}
				    if (msmdata3[k,1]==msmdata3[k+2,1]) {				
					if (msmdata3[k+1,2]-msmdata3[k,2]>cut & msmdata3[k+2,2]-msmdata3[k+1,2]>cut) {
								x<-msmdata3[k,]
								x[2]<-cut
								y<-rbind(y,x)
								x<-msmdata3[k+1,]
								y<-rbind(y,x)
								x<-msmdata3[k+2,]
								x[2]<-msmdata3[k+2,2]-cut
								y<-rbind(y,x)}					
					if (msmdata3[k+1,2]-msmdata3[k,2]>cut & msmdata3[k+2,2]-msmdata3[k+1,2]<=cut) {	
								x<-msmdata3[k,]
								x[2]<-cut
								y<-rbind(y,x)
								x<-msmdata3[k+1,]
								y<-rbind(y,x)}									
					if (msmdata3[k+1,2]-msmdata3[k,2]<=cut & msmdata3[k+2,2]-msmdata3[k+1,2]>cut) {
								x<-msmdata3[k+1,]
								x[2]<-cut
								y<-rbind(y,x)
								x<-msmdata3[k+2,]
								y<-rbind(y,x)}}}}	

			   msmcutwo<-y[2:nrow(y),]
 			   cat("\n")
			   cat("ULTIMO RAMO","\n")
			   msmcutwo<-data.frame(msmcutwo)
			   covariate<-matrix(data=NA,nrow=nrow(msmcutwo),ncol=ncov)
			   covariate<-data.frame(covariate)

			if (ncov!=0) { 	covariate[,1]<-msmcutwo[,formula[1]-1]
					names(covariate)[1]<-c(names(msmdata)[formula[1]-1])}
			if (ncov>1) { for (k in 2:ncov) { 	covariate[,k]<-msmcutwo[,formula[k]-1]
							names(covariate)[k]<-c(names(msmdata)[formula[k]-1])}}
		   	if (ncov!=0) {
					covu<-vector(length=ncov)
					covu<-names(covariate)
					covar <- as.formula(paste("~", paste(covu, collapse= "+")))}
			   cat("\n")
			   cat("*** convergence ***","\n")
			   if (ncov==0) { 	if (exacttimes==TRUE) stanf<-msm(msmcutwo[,3]~msmcutwo[,2],subject=msmcutwo[,1],qmatrix=qmat,inits=init,death=TRUE,method="BFGS",control=list(trace=2),exacttimes=T)
					else stanf<-msm(msmcutwo[,3]~msmcutwo[,2],subject=msmcutwo[,1],qmatrix=qmat,inits=init,death=TRUE,method="BFGS",control=list(trace=2),exacttimes=F,tunit=timeunit)}

			  if (ncov!=0) { if (exacttimes==TRUE) stanf<-msm(msmcutwo[,3]~msmcutwo[,2],subject=msmcutwo[,1],qmatrix=qmat,data=msmcutwo,covariates=covar, inits=inits, death=TRUE, method="BFGS", control=list(trace=2),exacttimes=T,fixedpars=fixedpars)
					else stanf<-msm(msmcutwo[,3]~msmcutwo[,2],subject=msmcutwo[,1],qmatrix=qmat,data=msmcutwo,covariates=covar, inits=inits,death=TRUE,method="BFGS", control=list(trace=2),exacttimes=F,fixedpars=fixedpars,tunit=timeunit)}

			   cat("\n")
			   cat("-2*Log-likelihood:",stanf$minus2loglik,"\n")
			   cat("\n")
			   cat("***** Estimated coefficients *****","\n")
			   print(coef.msm(stanf))
		               print(hazard.msm(stanf))
			   cat("\n")
			  
			   wald.tdc.surv(stanf,nstate,ncov,fixed)			
			   cat("\n")	
			   
			   cat("***** Sojourn times *****","\n")					
			   sojourn.msm(stanf)
			   print(sojourn.msm(stanf))
			   cat("\n")				
					
			   if (covmat==TRUE) { cat("*** variance-covariance matrix *** ","\n")
			       		         print(stanf$covmat)
					         cat("\n") }
			   newmodel[4]<-0
			   cat("\n")}
return(msmdata)
}

tdc.msm<-function(mydata = data.frame(),formula=parent.frame(),fixedpars=parent.frame(), models= “TDCM”, cut=0,prevalences=FALSE, writedata=FALSE,covmat=FALSE,graphcov=0, exacttimes=T, surv.plot=FALSE,plot.trans=FALSE,timeunit=1)
{
# Function for the analysis of survival data with time-dependent covariates.
# given any database in the form of counting Processes, this function adapt the data
# to use any of the following models: 
# (1) Time-dependent Cox regression model; 
# (2) Cox Markov model; 
# (3) multi-state homogeneous Markov model; 
# (4) Non-homogeneous model; 
# (5) Cox semi-Markov model. 
# This software can accommodate any multi-state model with one initial state, a finite number of 
# transient states corresponding to different values of the time-dependent covariate and one 
# absorbent state. In the particular situation of only one transient state, we obtain the 
# so-called progressive illness-death model.

require(survival)
require(splines)
require(msm)

aux1<-as.matrix(formula)
ncov<-nrow(aux1)

if (missing(ncov)) stop("Missing number of covariates")
if (missing(mydata)) stop("Missing data")
fixed<-1
if (missing(fixedpars)) {	fixedpars<-NULL
			fixed<-0}
model<-vector(length=6)
if (models == “TDCM”) model<-c(1,0,0,0,0,0)
if (models == “CMM”) model<-c(0,1,0,0,0,0)
if (models == “HMM”) model<-c(0,0,1,0,0,0)
if (models == “NHM”) model<-c(0,0,0,1,0,0)
if (models == “CSMM”) model<-c(0,0,0,0,1,0)
if (models == “NPM”) model<-c(0,0,0,0,0,1)
if (models == “ALL”) model<-c(1,1,1,1,1,1)

if (missing(model)) model<-c(1,0,0,0,0,0)
if (missing(prevalences)) prevalences<-FALSE
if (missing(writedata)) writedata<-FALSE
if (missing(covmat)) covmat<-FALSE
if (ncov>7) stop("Too many covariates")
if (missing(exacttimes)) exacttimes<-TRUE
if (missing(graphcov)) graphcov<-99
if (graphcov==0) graphcov<-0
if (missing(surv.plot)) surv.plot<-FALSE
if (missing(cut) & model[4]==1) stop("Missing cut")

nstate<-max(mydata[,5])+2
newmodel<-vector(length=6)
nmsmdata<-nrow(mydata)+max(mydata[,1])   
ngr<-vector(length=6)
ngr<-c(0,0,0,0,0,0)

if ( graphcov!=99 & surv.plot==TRUE) {
	ifelse (model[1]==1,ngr[1]<-2,ngr[1]<-0)
	ifelse (model[2]==1,ngr[2]<-2*nstate-2,ngr[2]<-0)
#	if (model[2]==1 & plot.trans==TRUE) ngr[2]<-ngr[2]+nstate
	ifelse (model[3]==1, ngr[3]<-1,ngr[3]<-0)
	if (model[3]==1 & plot.trans==TRUE) ngr[3]<-ngr[3]+nstate-1
	if (model[4]==1 & plot.trans==TRUE) ngr[4]<-ngr[4]+1
	ifelse (model[5]==1,ngr[5]<-2*nstate-2,ngr[5]<-0)
#	if (model[5]==1 & plot.trans==TRUE) ngr[5]<-ngr[5]+2*nstate-2
	ifelse (model[6]==1, ngr[6]<-1,ngr[6]<-0)
	if (model[6]==1 & plot.trans==TRUE) ngr[6]<-ngr[6]+nstate
	ngrsum<-ngr[1]+ngr[2]+ngr[3]+ngr[4]+ngr[5]+ngr[6]
	ngraph<-ceiling((ngrsum)/2)
	if (ngrsum==9 | ngrsum==8) op<-par(mfrow=c(3,3),mar=.1+c(2,2,3,1))
	if (ngrsum==1 | ngrsum==0) op<-par(mfrow=c(1,1),mar=.1+c(2,2,3,1))
	if (ngrsum==2) op<-par(mfrow=c(2,1),mar=.1+c(2,2,3,1))
	if (ngrsum!=1 & ngrsum!=2 & ngrsum!=8 & ngrsum!=9) op<-par(mfrow=c(ngraph,2),mar=.1+c(2,2,3,1))}

if ( graphcov!=99 & surv.plot!=TRUE) {
	ifelse (model[1]==1,ngr[1]<-1,ngr[1]<-0)
	ifelse (model[2]==1,ngr[2]<-2*nstate-3,ngr[2]<-0)
	if (model[3]==1 & plot.trans==TRUE) ngr[3]<-nstate-1
	if (model[4]==1 & plot.trans==TRUE) ngr[4]<-1
	ifelse (model[5]==1,ngr[5]<-2*nstate-3,ngr[5]<-0)
	ngrsum<-ngr[1]+ngr[2]+ngr[3]+ngr[4]+ngr[5]+ngr[6]
	ngraph<-ceiling((ngrsum)/2)
	if (ngrsum==9 | ngrsum==8) op<-par(mfrow=c(3,3),mar=.1+c(2,2,3,1))
	if (ngrsum==1 | ngrsum==0) op<-par(mfrow=c(1,1),mar=.1+c(2,2,3,1))
	if (ngrsum==2) op<-par(mfrow=c(2,1),mar=.1+c(2,2,3,1))
	if (ngrsum!=1 & ngrsum!=2 & ngrsum!=8 & ngrsum!=9 & ngrsum!=0) op<-par(mfrow=c(ngraph,2),mar=.1+c(2,2,3,1))}

if ( graphcov==99 & surv.plot==TRUE) {
	ifelse (model[1]==1,ngr[1]<-1,ngr[1]<-0)
	ifelse (model[2]==1,ngr[2]<-1,ngr[2]<-0)
#	if (model[2]==1 & plot.trans==TRUE) ngr[2]<-ngr[2]+nstate
	ifelse (surv.plot==TRUE & model[3]==1, ngr[3]<-1,ngr[3]<-0) 
	if (model[3]==1 & plot.trans==TRUE) ngr[3]<-ngr[3]+nstate-1
	ifelse (plot.trans==TRUE & model[4]==1, ngr[4]<-1,ngr[4]<-0)
	ifelse (model[5]==1,ngr[5]<-1,ngr[5]<-0)
#	if (model[5]==1 & plot.trans==TRUE) ngr[5]<-ngr[5]+2*nstate-2
	ifelse (model[6]==1,ngr[6]<-1,ngr[6]<-0)
	if (model[6]==1 & plot.trans==TRUE) ngr[6]<-ngr[6]+nstate
	ngrsum<-ngr[1]+ngr[2]+ngr[3]+ngr[4]+ngr[5]+ngr[6]
	ngraph<-ceiling((ngrsum)/2)
	if (ngrsum==9 | ngrsum==8) op<-par(mfrow=c(3,3),mar=.1+c(2,2,3,1))
	if (ngrsum==1 | ngrsum==0) op<-par(mfrow=c(1,1),mar=.1+c(2,2,3,1))
	if (ngrsum==2) op<-par(mfrow=c(2,1),mar=.1+c(2,2,3,1))
	if (ngrsum!=1 & ngrsum!=2 & ngrsum!=8 & ngrsum!=9 & ngrsum!=0) op<-par(mfrow=c(ngraph,2),mar=.1+c(2,2,3,1))}

if ( graphcov==99 & surv.plot==FALSE & plot.trans==TRUE) {
#	if (model[2]==1) ngr[2]<-nstate
#	if (model[5]==1) ngr[5]<-2*nstate-2
	if (model[3]==TRUE) ngr[3]<-nstate-1
	if (model[4]==TRUE) ngr[4]<-1
	if (model[6]==1) ngr[6]<-nstate
	ngrsum<-ngr[1]+ngr[2]+ngr[3]+ngr[4]+ngr[5]+ngr[6]
	ngraph<-ceiling((ngrsum)/2)
	if (ngrsum==9 | ngrsum==8) op<-par(mfrow=c(3,3),mar=.1+c(2,2,3,1))
	if (ngrsum==1 | ngrsum==0) op<-par(mfrow=c(1,1),mar=.1+c(2,2,3,1))
	if (ngrsum==2) op<-par(mfrow=c(2,1),mar=.1+c(2,2,3,1))
	if (ngrsum!=1 & ngrsum!=2 & ngrsum!=8 & ngrsum!=9 & ngrsum!=0) op<-par(mfrow=c(ngraph,2),mar=.1+c(2,2,3,1))}

if ( graphcov!=99 & surv.plot==FALSE & plot.trans==TRUE) {
	ifelse (model[1]==1,ngr[1]<-1,ngr[1]<-0)
#	ifelse (model[2]==1,ngr[2]<-3*nstate-2,ngr[2]<-0)
	ifelse (model[2]==1,ngr[2]<-2*nstate-2,ngr[2]<-0)
	ifelse (model[3]==1,ngr[3]<-nstate-1,ngr[3]<-0)
	ifelse (model[4]==1,ngr[4]<-1,ngr[4]<-0)
#	ifelse (model[5]==1,ngr[5]<-4*nstate-4,ngr[5]<-0)
	ifelse (model[5]==1,ngr[5]<-2*nstate-2,ngr[5]<-0)
	if (model[6]==1) ngr[6]<-nstate
	ngrsum<-ngr[1]+ngr[2]+ngr[3]+ngr[4]+ngr[5]+ngr[6]
	ngraph<-ceiling((ngrsum)/2)
	if (ngrsum==9 | ngrsum==8) op<-par(mfrow=c(3,3),mar=.1+c(2,2,3,1))
	if (ngrsum==1 | ngrsum==0) op<-par(mfrow=c(1,1),mar=.1+c(2,2,3,1))
	if (ngrsum==2) op<-par(mfrow=c(2,1),mar=.1+c(2,2,3,1))
	if (ngrsum!=1 & ngrsum!=2 & ngrsum!=8 & ngrsum!=9 & ngrsum!=0) op<-par(mfrow=c(ngraph,2),mar=.1+c(2,2,3,1))}

#print(date())

# TIME-DEPENDENT COX REGRESSION MODEL
if (model[1]==1) {tdcrm.tdc.surv(mydata,ncov,formula,covariate,newmodel,covmat,graphcov,surv.plot,plot.trans)}

msmdata<-matrix(data=NA,nrow=nrow(mydata)+max(mydata[,1]),ncol=15)
covariate<-matrix(data=NA,nrow=nrow(mydata)+max(mydata[,1]),ncol=ncov)
mydata[nrow(mydata)+1,]<-c(rep(1,ncol(mydata)))
vec<-vector(length=ncov)
ncovdata<-ncol(mydata)-5

if (model[2]==1 | model[3]==1 | model[4]==1 | model[5]==1 | model[6]==1) msmdata<-msmdata(mydata,ncov,ncovdata,nstate)
if (model[2]==1 | model[5]==1 | model[6]==1) sobrevi<-nelsonomp(msmdata,nstate)

# COX MARKOV (SEMI-MARKOV) MODEL 
if (model[2]==1|model[5]==1){	cmm.tdc.surv(mydata,msmdata,nstate,ncov,covar,formula,graphcov,newmodel,covmat,model,writedata,coxdata,surv.plot,plot.trans)}

if (model[3]==1|model[4]==1){	multistate(mydata,msmdata,ncov,ncovdata,nstate,formula,covmat,prevalences,writedata,graphcov,model,newmodel,fixedpars,fixed,cut,exacttimes,surv.plot,plot.trans,timeunit)}

if (model[6]==1) {	#msmdata<-msmdata(mydata,ncov,ncovdata,nstate)
			newmodel[6]<-1
			cat("*****************************************************************************","\n")
			cat("           ***** NON-PARAMETRIC MULTI-STATE MODEL *****              ","\n")
	                        cat("*****************************************************************************","\n")
			cat("\n")
			#sobrevi<-nelsonomp(msmdata,nstate)
			if (plot.trans==TRUE | surv.plot==TRUE) plot.tdc.surv(stanf,mydata,ncov,graphcov,newmodel,newcov,covariate,formula,surv.plot,sobrevi,now,nstate,transition,plot.trans) 
			nonparam<-matrix(ncol=ncol(sobrevi),nrow=((nrow(sobrevi)+1)/2))
			for (k in 1:nrow(nonparam)) nonparam[k,]<-sobrevi[2*k-1,] 
			nonparam<-data.frame(nonparam)
			x<-matrix(ncol=9,nrow=2)
			names(x)<-c("P12","P13","P14","P15","P16","P17","P18","P19","P1,10")
			names(nonparam)<-c("Times","Survival","P11",names(x)[1:(nstate-1)])
			print(nonparam)
			newmodel<-0
			if (writedata==TRUE) { 	cat("*** FORMATED DATA *** ","\n")
			       			write.table(nonparam,file="C:\nonparam.txt",sep=",",na = "NA",dec = ".",col.names = TRUE)}
		}
if (ngr[1]+ngr[2]+ngr[3]+ngr[4]+ngr[5]+ngr[6]!=0) { if ( graphcov!=99 | surv.plot==TRUE | plot.trans==TRUE) par(op) }
}

#Example: 

data(stanford)

tdc.msm(stanford,formula=c(6,7,8),models=“ALL”,graphcov=1,cut=60,prevalences=T)

tdc.msm(stanford,formula=c(6,7,8),models=“CMM”,graphcov=1,cut=60,prevalences=F, surv.plot=T)

tdc.msm(stanford,formula=c(6,7,8),models=“ALL”,graphcov=1,cut=90,prevalences=T, fixedpars=c(5,7,8,9))

tdc.msm(stanford,formula=c(6,7,8),models= “CSMM”,surv.plot=TRUE,graphcov=1,cut=90, plot.trans=TRUE)

tdc.msm(stanford ,formula=c(6,7,8),models= “CMM”,surv.plot=TRUE,graphcov=1)