https://github.com/cran/emplik
Tip revision: 1a00a3dce7a5ac3854ad4978598dbc8ab597c49f authored by Mai Zhou on 23 October 2010, 00:00:00 UTC
version 0.9-6
version 0.9-6
Tip revision: 1a00a3d
ROCnp.R
#### Modified by Mai Zhou, 9/2007, 10/2007.
#### For two samples of right censored data.
#### Both samples use nonparametric likelihood = EL.
#### We test the Ho: ROC curve R(t0) = b0
#### Or testing the (1-b0)th quantile of sample one equals to
#### the (1-t0)th quantile of sample two.
#### There could be several versions: either use optimize() or optim() to
#### minimize over the value of commom quantile (the cstar value in output)
#### and either use el.cen.EM2() or use emplikH.disc() to compute EL.
#### It seems the el.cen.EM2() version is more stable.
####
#### input: t1; right censored times, sample 1.
#### d1; censoring status, d1=1 means uncensored.
ROCnp <- function(t1, d1, t2, d2, b0, t0) {
if ( length(b0) != 1) stop ("check length of b0")
if ( length(t0) != 1) stop ("check length of t0")
if ( b0 >=1 | b0 <=0) stop ("check the value of b0")
if ( t0 >=1 | t0 <=0) stop ("check the value of t0")
tempnp2 <- WKM(x=t2, d=d2)
place2 <- sum(tempnp2$surv >= t0)
c2 <- tempnp2$times[place2]
tempnp1 <- WKM(x=t1, d=d1)
place1 <- sum(tempnp1$surv > b0)
c1 <- tempnp1$times[place1]
if(c2 <= c1) c1 <- tempnp1$times[place1 +1] #### need to do this??
#llr <- function(const, t1, d1, t2, d2, b0, t0) {
# npllik2 <- emplikH.disc(x=t2, d=d2, K=log(t0),
# fun=function(x,theta){as.numeric(x <= theta)},
# theta=const)$"discrete.-2LLR"
# npllik1 <- emplikH.disc(x=t1, d=d1, K=log(b0),
# fun=function(x,theta){as.numeric(x <= theta)},
# theta=const)$"discrete.-2LLR"
# return( npllik2 + npllik1 )
# }
#
### distribution version
llr <- function(const, t1, d1, t2, d2, b0, t0) {
npllik1 <- el.cen.EM2(x=t1, d=d1,
fun= function(x,theta){as.numeric(x <= theta)},
mu=1-b0, theta=const)$"-2LLR"
npllik2 <- el.cen.EM2(x=t2, d=d2,
fun= function(x,theta){as.numeric(x <= theta)},
mu=1-t0, theta=const)$"-2LLR"
return(npllik1+npllik2)
}
##temp <- optim(par=(c2+c1)/2, fn=llr, method="L-BFGS-B", lower=min(c2,c1),
## upper=max(c2,c1), t1=t1,d1=d1,t2=t2,d2=d2,b0=b0,t0=t0)
temp <- optimize(f=llr, lower=min(c2,c1), upper=max(c2,c1),
t1=t1, d1=d1, t2=t2, d2=d2, b0=b0, t0=t0)
##cstar <- temp$par
cstar <- temp$minimum
##val <- temp$value
val <- temp$objective
list("-2LLR"=val, cstar=cstar)
}