R/lrt.confints.R In km.ci: Confidence intervals for the Kaplan-Meier estimator

```"lrt.confints" <-
function(time,status,t0,alpha=0.05) {
ftimes<-sort(unique(time[status==1]))
k<-length(ftimes)
dvec<-rep(1,k)
nvec<-dvec
jmax<-0
for(i in 1:k) {
dvec[i]<-sum(time==ftimes[i] & status==1)
nvec[i]<-sum(time>=ftimes[i])
if(t0 >= ftimes[i])
jmax<-i
}
nsub<-nvec[1:jmax]
dsub<-dvec[1:jmax]
theta<-qchisq(1-alpha,1)

#handle no deaths specially
if( jmax == 0 || sum(dsub) == 0 )
return(list(lower=exp(-theta/(2*max(nsub))),upper=1))
#else
lrt<-function(lambda,opar) {
nsub<-opar\$nsub
dsub<-opar\$dsub
opar\$theta -2*sum(nsub*log(1+lambda/nsub)-(nsub-dsub)*log(1
+ lambda/(nsub-dsub)))
}
plambda<-function(lambda,opar) {
nsub<-opar\$nsub
dsub<-opar\$dsub
prod(1-dsub/(nsub+lambda))
}
t1<-sum(1/(nsub-dsub)-1/nsub)
l1<- -sqrt(theta/t1)
#the likelihood is undefined if l1 < lbd
lbd<- dvec[jmax]-nvec[jmax]
tol2 <- 1/100
#if( nvec[jmax]-dvec[jmax]+l1 < 0 ) {
# l1<-(dvec[jmax]-nvec[jmax]-0.001)/10
#}
parlist<-list(nsub=nsub,dsub=dsub,theta=theta)
for(i in 1:10) {
if( l1 < lbd ) {
l1<-lbd + tol2
tol2<- tol2/2
}
v1 <- lrt(l1,parlist)
if(v1 < 0 ) break
else {
slope<-(theta - v1)/l1
l1<-(theta + 0.5 )/slope
}
}
lower<-bisect(lrt,parlist,l1,0)
l1<- -l1
for(i in 1:10) {
v1<-lrt(l1,parlist)
if(v1 < 0 ) break
else {
slope<-(theta - v1)/l1
l1<-(theta + 0.5 )/slope
}
}
upper<-bisect(lrt,parlist,0,l1)
lp<-plambda(lower,parlist)
up<-plambda(upper,parlist)
return(list(lower=lp,upper=up))
}
```

Try the km.ci package in your browser

Any scripts or data that you put into this service are public.

km.ci documentation built on May 2, 2019, 2:46 a.m.