# Program calculates a Cox model parameter estimate for the effect of single    #
# covariate on survival time, when the survival time is interval censored.             #
# To use this program the data for each observation should consist of,               #
# in order: obs. no. , left endpoint for survival time, right endpoint for                   #
# survival time, covariate value.  The program calls a C subroutine from              #                               # the object file : met99.o					             # 
# The methodology upon which this program is based treats right censoring      #
# as a special case of interval censoring. Right censored observations should   #
# be given a left endpoint endpoint equal to the censoring time and a right        #
# endpoint equal to a number greater than the largest right endpoint or             #    
# censoring time for any observation. The number on line 62 ( if rend[i] == )      #
# should be changed to this number.                                                                   #  
# To run program:	                                                                                           #
# (1) Copy onecov.s and met99.o into your directory                                         #
# (2) Have data in a file called int.in                                                                   #
# (3) Go into Splus	                                                                                          #
# (4) At the prompt type: >!dyn.load("met99.o")		                          #
# (5)  >motif()						          #
# (6) > myfunc_source("onecov.s")				          #
# (7) > myfunc(s,m)					                          #
# s & m are the number of shuffles and number of imputations required by the   #
# method. In general s should be higher with data sets for data with larger       # 
# sample size and wider intervals. m should be higher for data with wider intervals. #
# Appropriate adjustments in these parameters can be made by examining plots  #
# of the convergence of betahat and the serial correlation of the score statistics.  #
# If there is too much between iteration fluctuation in the parameter estimates m  #
# should be increased, if too much serial correlation in the score statistics,            #
# s should be increased.						#    	
	
function(s, m)
{
# m is the no. of iterations, s is the no. of shuffles,#
# bo is the initial value for beta#

# Read in data #	
		
zzlist <- list(o=0,l=0,r=0,c=0)
	
zz <- scan("int.in", what = zzlist)	
	
	n <- length(zz$l)
	diff <- NULL
	diff[1] <- 0.002
	par(mfrow = c(2, 2), oma = c(0, 0, 4, 0))
	betahat <- NULL
	midpt <- NULL
	expo <- NULL
	inf <- NULL
	cloglik <- rep(0, m)
	cdloglik <- rep(0, m)
	totll <- NULL
	totdll <- NULL
	totd2ll <- NULL
	totdllsq <- NULL
	var <- NULL
	sterr <- NULL
	inft1 <- NULL
	inft2 <- NULL
	inft3 <- NULL
                zstat <- NULL
                norpr <- NULL
	beta <- NULL
                retmat <- NULL
	tloglik <- 1.0
	tdloglik <- 1.0
	td2loglik <- 1.0
	tdlogliksq <- 1.0

	q <- 1
	rend <- zz$r
	lend <- zz$l
	x <- zz$c
	cens <- NULL

######### Get initial estimate using midpoint imputation #######	
	
	
	for(i in 1:n) {
		if(rend[i] == 99) {
			midpt[i] <- lend[i]
	                                 cens[i] <- 0
		}
		else {
			midpt[i] <- (rend[i] + lend[i])/2
	                                cens[i] <- 1
		}
	}
	midorder <- order(midpt)
	reno <- rend[midorder]
	leno <- lend[midorder]
	xo <- x[midorder]
	midpto <- midpt[midorder]

     middat <- matrix(0,3,n)
     middat <- cbind(midpto,1,xo)
     midout <- coxreg(middat[,1],status = cens,middat[,3])

     cat("Cox analysis using midpoint imputation method")
     cat("\n")
     print(midout)
     cat("\n\n")

     midbeta <- midout$coef


######################################################

	
	dloglik <- rep(0,m)
		
	beta[1] <- midbeta


######### BEGIN MONTE CARLO EM PROCEDURE ####################


	
	while(q < 9) {
		expo <- exp(beta[q] * xo)
		rando <- runif(1) * (100000)
	        betahat <- beta[q]
		zzz <- .C("execute",
			betahat = as.double(betahat),
			as.integer(n),
			as.integer(m),
			as.integer(s),
			as.integer(q),
			leno = as.integer(leno),
			reno = as.integer(reno),
			xo = as.double(xo),
			expo = as.double(expo),
			as.double(rando),
			dloglik = as.double(dloglik),
	                tloglik = as.double(tloglik),
	                tdloglik = as.double(tdloglik),
	                td2loglik = as.double(td2loglik))
	               
		leno <- zzz$leno
		reno <- zzz$reno
		xo <- zzz$xo
	                dloglik <- zzz$dloglik
		expo <- zzz$expo
		totll[q] <- zzz$tloglik
		totdll[q] <- zzz$tdloglik
		totd2ll[q] <- zzz$td2loglik
		totdllsq[q] <- sum((dloglik)^2)
	        beta[q+1] <- zzz$betahat

		inft1[q] <- -(1/m) * totd2ll[q]
		inft2[q] <- (1/m) * totdllsq[q]
		inft3[q] <- (((1/m)*(totdll[q]))^2)
	#	cat("information term1        info term2            info term3")
	#	cat("\n")
	#	cat(inft1[q], inft2[q], inft3[q], sep = "        ")
		cat("\n")
                                inf[q] <- inft1[q] - inft2[q] + inft3[q]
		var[q] <- 1/(inf[q])
		sterr[q] <- (var[q])^0.5
		diff[q + 1] <- beta[q + 1] - beta[q]
		cat("\n")
		cat("iteration     betahat                 change                   std. err.")
		cat("\n")
		cat(q, beta[q + 1], diff[q + 1], sterr[q], sep = "       ")
		cat("\n")
		
	 
		acf.lik <- acf(dloglik, 10, "correlation")
		cat("\n")
                zstat[q] <- beta[q+1]/sterr[q]
                norpr[q] <- pnorm(zstat[q])
        if(q > 7){  
	
               cat("z = ", zstat[q], sep = "     ")
                cat("\n")
                cat("p-value = ", norpr[q], sep = "     ")
                cat("\n")
		cat("correlation of log-likelihoods, lags 1,2,3 ...:")
		cat("\n")
		cat(acf.lik$acf)
		cat("\n")}
		plot(beta, type = "o", main = "beta vs. iterations")
		plot(sterr, type = "o", main = 
			"std. error of beta vs. iterations")
	

		plot(diff, type = "o", main = "difference from preceding beta vs. iterations")
		abline(0, 0)

		q <- q + 1
	}


}