# 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)
{
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
}
```

## Try the thregI package in your browser

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

thregI documentation built on May 2, 2019, 6:27 a.m.