#Fitting survival model with covaraite subject to limit of detection
library(foreign)

# Uses Simpson's rule; function values are elements of vec
# length of vec should be odd
simpson <- function(vec, a, b)
{

        n <- length(vec)
        m <- n - 1      # m is the number of intervals

        h <- (b - a)/m  # h is the length of each interval
        d <- c(1, rep(c(4, 2), m/2 - 1), 4, 1)

        integral <- (h/3) * sum(d * vec)
        integral
}

surv.like.corrected <- function(dat=data0, initial=c(0, 0, 0, 0, 0, 0, 0), c0=2)
{
	y <- dat$Time
	x <- dat$x
	x1<- dat$age
	x2<- dat$sex
	x3<- dat$black
	x4<- dat$charlson
	x5<- dat$psi

	n <- length(y)
	delta <- dat$Event
	r <- dat$r
	
		obj.fn <- function(theta)	
	{
		beta0 <- theta[1]
		beta1 <- theta[2]
		
		beta2 <- theta[3]
		beta3 <- theta[4]
		beta4 <- theta[5]
		beta5 <- theta[6]
		beta6 <- theta[7]

		pi.r <- 1 - pnorm((c0-theta[8])/theta[9])
		log.like  <- 0
		for (i in 1:n)
		{
			if (r[i]==1)
			{
			
				# hazard:        lambda * exp(x[i] * beta0)
				# survival:      exp(- exp(x[i] * beta0) * lambda * y[i])
				# likelihood:    {h.y ^ delta[i]} * S.y * (1 - pi.r) * f.x
					
				f.x <- dnorm(x[i], mean=theta[7], sd=theta[8])
				eta.i <- beta0 + beta1 * x[i] + beta2 * x1[i] + beta3 * x2[i] + beta4 * x3[i] + beta5 * x4[i] + beta6 * x5[i] 
 
				log.like.i <- - (delta[i] * eta.i - exp(eta.i) * y[i] + log((pi.r)) + log(f.x))

			} else {

				low.x <- -20
				x0 <- seq(low.x, c0, length=51)
				f.x0 <- dnorm(x0, mean=theta[8], sd=theta[9]) 
				eta0 <- beta0 + beta1 * x0 + beta2 * x1[i] + beta3 * x2[i] + beta4 * x3[i] + beta5 * x4[i] + beta6 * x5[i] 
				#browser()
				log.like.i <- -log(simpson( exp(delta[i] * eta0 - exp(eta0) * y[i] + log((1 - pi.r)) + log(f.x0)),  low.x, c0))
			}
			log.like <- log.like + log.like.i
			
 		}
		log.like
		print(log.like)

	}
	initial= c(initial, mean(x[r==1]), sqrt(var(x[r==1])))
	surv.fit <- optim(par=initial, fn=obj.fn, hessian=TRUE, method = c("L-BFGS-B"))
	cat(surv.fit$message, "\n")
	estimate <- surv.fit$par
	std.err <- sqrt(diag(solve(surv.fit$hessian)))	
	objective <- surv.fit$value
	list(estimate=estimate, std.err=std.err, ci=cbind(estimate-qnorm(.975)*std.err,estimate+qnorm(.975)*std.err),
 		wald.pval=pchisq((estimate/std.err)^2,df=1, lower.tail =FALSE),objective=objective)
}

dta=read.dta("I:\\_Stat.Research\\CurrentResearch\\SurvivalModel_LOD\\RealData\\Data\\GenIMS_IL6_day1.dta")
formatdta=data.frame(dta[,c("day", "age", "sex", "black", "charlson", "psi")], x=log(dta$il6),r=ifelse(dta$il6<0,0,1), Time=dta$Days90, Event=dta$day90_status)
formatdta=formatdta[!(is.na(dta$il6) | is.na(dta$Days90)),] # removing rows with missing data
corected=surv.like.corrected(dat=formatdta, c0=log(2))

corected