https://github.com/cran/emplik
Revision a09ab78cfdae5c86da702f74c604cda1abde8ddf authored by Mai Zhou on 08 August 1977, 00:00:00 UTC, committed by Gabor Csardi on 08 August 1977, 00:00:00 UTC
1 parent 1d6b37a
Tip revision: a09ab78cfdae5c86da702f74c604cda1abde8ddf authored by Mai Zhou on 08 August 1977, 00:00:00 UTC
version 0.5-1
version 0.5-1
Tip revision: a09ab78
emplikdisc.test.R
########################################################
############ emplikdisc.test() #########################
########################################################
emplikdisc.test <- function(x, d, 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)
otime <- temp$time # 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)
### 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(" too far away ")
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 60 ?
stop("given theta value is too far away from theta0")
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.-2logemlikRatio"=loglik, lambda=lam, times=otime,
jumps=odti/rPlgti)
}
Computing file changes ...