R/thregI.R

In thregI: Threshold Regression for Interval-Censored Data with a Cure Rate Option

## The following is the definition of the thregI function for interval-censored data
"thregI" <-
function (formula,data)
{
 ### Read in all arguments
 cl <- match.call()
 indx <- match(c("formula", "data"), names(cl), nomatch=0)
 if (indx[1] ==0) stop("A formula argument is required")
 mf<- cl[c(1, indx)]
 f <- Formula(formula)
 f1<-formula(f,lhs=1)
 f1<-Formula(f1)
 mf[[1]] <- as.name("model.frame")
 mf$formula <- if(missing(data)) terms(f1) else terms(f1, data=data)
 mf$formula <- f1
 mf <- eval(mf, parent.frame())
 if (nrow(mf) ==0) stop("No (non-missing) observations")
 Terms <- terms(mf)
 Y <- model.extract(mf, "response")
 if (!inherits(Y, "Surv")) stop("Response must be a survival object")
 type <- attr(Y, "type")
 if (type!='interval') stop(paste("thregI package only support \"", type, "\" survival data", sep=''))
 f2<-formula(f1, lhs = 0)
 if (length(f2[[2]])!=3) stop(paste("Predictors for both lny0 and mu should be specified"))
 x_lny0<-model.matrix(f1, data, rhs=1)
 x_mu<-model.matrix(f1, data, rhs=2)
 left  <- Y[,1]
 right <- Y[,2]
 delta <- Y[,3] #delta=1:event, delta=3:left/interval delta=0:right
                #no delta=2 for left=NA since our dateset no NA
 ####################################
 #(0,R]---->delta=3, left=0, right=R
 #(a,b]---->delta=3, left=a, right=b
 #(L,Inf]-->delta=0, left=L, right=1
 #(a,a]---->delta=1, left=a, right=1
 ###################################
 delta3=matrix(0,length(left),1)
 # fix the exact times .............................................................
 for (i in 1 :length(left))
 {
   if (delta[i]==1) {right[i]=left[i]  #right=a instead of 1
   delta3[i]=1}      #delta3=1 for exact time
 }
 #right_max=10*max(right[right!=Inf])
 right_max=10*max(right)
 delta1=matrix(0,length(left),1)
 delta2=matrix(0,length(left),1)
 for (i in 1 :length(left)){
   #if (right[i]==1) { right[i]=right_max }
   if (delta[i]==0) {right[i]=right_max}
   else if (left[i]==0) { delta1[i]=1 }   #delta1=1 for left censoring
   else if (delta[i]==3){ delta2[i]=1 }   #delta2=1 for interval censoring
 }                                        #delta1=delta2=delta3=0 for right censoring
 ####################################
 #(0,R]---->delta1=1, left=0, right=R
 #(a,b]---->delta2=1, left=a, right=b
 #(L,Inf]-->delta1=delta2=delta3=0, left=L, right=Inf
 #(a,a]---->delta3=1, left=a, right=a
 ###################################

 lny0<-function(para_lny0){x_lny0%*%para_lny0}
 mu<-function(para_mu){x_mu%*%para_mu}
 d<-function(para){
  para_lny0=para[1:length(dimnames(x_lny0)[[2]])]
  para_mu=para[(length(dimnames(x_lny0)[[2]])+1):(length(dimnames(x_lny0)[[2]])+length(dimnames(x_mu)[[2]]))]
  -mu(para_mu)/exp(lny0(para_lny0))
 }
 v<-function(para){
  para_lny0=para[1:length(dimnames(x_lny0)[[2]])]
  exp(-2*lny0(para_lny0))
 }
 su<-function(para){
  pnorm((1-d(para)*left)/sqrt(v(para)*left))-exp(2*d(para)/v(para))*pnorm(-(1+d(para)*left)/sqrt(v(para)*left))
 }
 sv<-function(para){
  pnorm((1-d(para)*right)/sqrt(v(para)*right))-exp(2*d(para)/v(para))*pnorm(-(1+d(para)*right)/sqrt(v(para)*right))
 }
 logdf<-function(para){
  -.5*(log(2*pi*v(para)*(right^3))+(d(para)*right-1)^2/(v(para)*right))
 }
 logf<-function(para) {
-sum(delta1*log(1-sv(para)))-sum(delta2*log(su(para)-sv(para)), na.rm = TRUE)-sum((1-delta1-delta2-delta3)*log(su(para)))-sum(delta3*logdf(para))
 }

 p<-rep(0,(length(dimnames(x_lny0)[[2]])+length(dimnames(x_mu)[[2]])))
 est<-nlm(logf, p, hessian = TRUE)

 names(est$estimate) <-c(paste("lny0:",dimnames(x_lny0)[[2]]),paste("  mu:",dimnames(x_mu)[[2]]))
 loglik = (-1)*est$minimum

 fit<-list(coefficients  = est$estimate,
          var    = solve(est$hessian),
          loglik = loglik,
          AIC    = (-2)*loglik+2*(length(dimnames(x_lny0)[[2]])+length(dimnames(x_mu)[[2]])),
          iter   = est$iterations,
          call   = cl,
          mf     = mf,
          lny0   = dimnames(x_lny0)[[2]],
          mu     = dimnames(x_mu)[[2]])
 class(fit) <- 'thregI'
 fit
}