emplikH.disc.R
########################################################
############ emplikH.disc() #########################
########################################################
emplikH.disc <- function(x, d, y= -Inf, K, fun,
tola=.Machine$double.eps^.25, theta)
{
n <- length(x)
if(n <= 2) stop("Need more observations")
if(length(d) != n ) stop("length of x and d must agree")
if(any((d!=0)&(d!=1))) stop("d must be 0/1's for censor/not-censor")
if(!is.numeric(x)) stop("x must be numeric values --- observed times")
#temp<-summary(survfit(Surv(x,d),se.fit=F,type="fleming",conf.type="none"))
#
newdata <- Wdataclean2(x,d)
temp <- DnR(newdata$value, newdata$dd, newdata$weight,y=y)
otime <- temp$times # only uncensored time? Yes.
orisk <- temp$n.risk
odti <- temp$n.event
###if the last jump is of size 1, we need to drop last jump from computation
last <- length(orisk)
if (orisk[last] == odti[last]) {
otime <- otime[-last]
orisk <- orisk[-last]
odti <- odti[-last]
}
######## compute the function g(ti, theta)
gti <- fun(otime,theta)
Kcenter <- sum(gti * log(1- odti/orisk))
### the constrain function. To be solved in equation later.
constr <- function(x, Konst, gti, rti, dti, n) {
rtiLgti <- rti + x*n*gti
OneminusdH <- (rtiLgti - dti)/rtiLgti
if( any(OneminusdH <= 0) ) stop("estimator not well defined")
sum(gti*log(OneminusdH)) - Konst }
##############################################################
differ <- constr(0, Konst=K, gti=gti, rti=orisk, dti=odti, n=n)
if( abs(differ) < tola ) { lam <- 0 } else {
step <- 0.2/sqrt(n)
if(abs(differ) > 200*log(n)*step ) #Why 200 ?
{ print( Kcenter )
stop("the given K value is too far away from K_0. \n Move K closer to the above value")
}
mini<-0
maxi<-0
######### assume the constrain function is increasing in lam (=x)
if(differ > 0) {
mini <- -step
while(constr(mini, Konst=K, gti=gti, rti=orisk, dti=odti, n=n) > 0
&& mini > -200*log(n)*step )
mini <- mini - step
}
else {
maxi <- step
while(constr(maxi, Konst=K, gti=gti, rti=orisk, dti=odti, n=n) < 0
&& maxi < 200*log(n)*step )
maxi <- maxi+step
}
if(constr(mini, Konst=K, gti=gti, rti=orisk, dti=odti, n=n)*constr(maxi,
Konst=K, gti=gti, rti=orisk, dti=odti, n=n) > 0 )
stop("given theta/K is/are too far away from theta0/K0")
# Now we solve the equation to get lambda, to satisfy the constraint of Ho
temp2 <- uniroot(constr,c(mini,maxi), tol = tola,
Konst=K, gti=gti, rti=orisk, dti=odti, n=n)
lam <- temp2$root
}
####################################################################
rPlgti <- orisk + n*lam*gti
loglik <- 2*sum(odti*log(rPlgti/orisk) +
(orisk-odti)*log(((orisk-odti)*rPlgti)/(orisk*(rPlgti-odti)) ) )
#?is that right? YES the -2log lik ratio.
# Notice the output time and jumps has less the last point.
list("discrete.-2LLR"=loglik, lambda=lam, times=otime,
jumps=odti/rPlgti)
}