.packageName <- "intcox"
"intcox" <-
function (formula = formula(data), data = parent.frame(), subset, 
    na.action, x = FALSE, y = TRUE, epsilon = 1e-04, itermax = 10000, 
    no.warnings = FALSE) 
{
    require(survival)
#
#   preparation
#
    copy.data <- data
    call <- match.call()
    m <- match.call(expand = FALSE)
    temp <- c("", "formula", "data", "copy.data", "subset", "na.action")
    m <- m[match(temp, names(m), nomatch = 0)]
    m[[1]] <- as.name("model.frame")
    Terms <- if (missing(copy.data)) 
        terms(formula)
    else terms(formula, data = copy.data)
    m$formula <- Terms
    m <- eval(m, parent.frame())
    Y <- model.extract(m, "response")
    if (!inherits(Y, "Surv")) 
        stop("Response must be a survival object")
    attr(Terms, "intercept") <- 1
    type <- attr(Y, "type")
    if (type != "interval") 
        stop("Invalid survival type (only interval censored data)")
    else {
        if (any(Y[, 3] != 3 && Y[, 3] != 0)) 
            stop("Invalid cens status")
        else {
            if (any(Y[, 3] == 3)) 
                Y <- cbind(Y[, 1:2], Y[, 3])
            else Y <- cbind(Y[, 1], Y[, 3])
        }
    }
    copy.data$left <- Y[, 1]
    copy.data$right <- Y[, 2]
    copy.data$cens <- Y[, 3]
    copy.data$mix <- ifelse(copy.data$cens == 3, copy.data$right, 
        copy.data$left)                             # start with the right ends if observed 
    ord <- order(copy.data$mix, 3 - copy.data$cens) #sort, first by right interval limit, then by censoring, that the Breslow-estimator can work
    copy.data <- copy.data[ord, ]
    X <- model.matrix(Terms, m)
    X <- X[, -1, drop = FALSE]
    covar <- X[ord, , drop = FALSE]
#
#   calculation
#
    if (length(attr(Terms, "term.label")) == 0) 
        stop("Invalid type of model")
    else intcox.new <- intcox.fit(formula, copy.data, covar, epsilon, 
        itermax)
#
#   output
#
    if (is.character(intcox.new)) {
        intcox.new <- list(fail = intcox.new)
        class(intcox.new) <- "coxph"
    }
    else {
        intcox.new$n <- nrow(Y)
        class(intcox.new) <- "coxph"
        intcox.new$terms <- Terms
        intcox.new$assign <- assign
        intcox.new$var <- matrix(NA, ncol = 1)
        intcox.new$wald.test <- NA
        intcox.new$score <- NA
        intcox.new$loglik <- intcox.new$likeli.vec[length(intcox.new$likeli.vec)]
        intcox.new$linear.predictors <- as.vector(intcox.new$coefficients * 
            t(covar))
        intcox.new$residuals <- NA
        intcox.new$means <- apply(covar, 2, mean)
        na.action <- attr(m, "na.action")
        if (length(na.action)) 
            intcox.new$na.action <- na.action
        if (x) {
            intcox.new$x <- covar
        }
        if (y) 
            intcox.new$y <- Y
    }
    intcox.new$formula <- formula(Terms)
    intcox.new$call <- call
    intcox.new$method <- NA
    if (no.warnings==FALSE) {
        if (intcox.new$termination == 4) 
            cat("inside precondition(s) at iteration = ", intcox.new$iter, 
                "not fulfilled\n")
        if (intcox.new$termination == 2) 
            cat("no improvement of likelihood possible, iteration = ", 
                intcox.new$iter, "\n")
        if (intcox.new$termination == 3) 
            cat("algorithm did not converge - maximum number of iteration reached, itermax = ", 
                intcox.new$iter, "\n")
    }
    return(intcox.new)
}
"intcox.breslow" <-
function (formula, data, covar)             # Breslow-Estimator
{
    lokal.cens <- data$cens/3
    formula.covar <- formula[[3]]
    formula <- as.formula(Surv(data$mix, lokal.cens) ~ .)
    formula[[3]] <- formula.covar
    fit <- coxph(formula, data)$coef        # Cox-regression coefficient
    e <- exp(t(as.matrix(fit)) %*% t(covar))
    event.num <- cumsum(data$cens == 3)
    hazard.rate0 <- NULL
    for (i in 1:max(event.num)) hazard.rate0 <- c(hazard.rate0, 
        1/sum(e[i <= event.num]))
    cumhaz <- cumsum(hazard.rate0)
    breslow.ret <- list(cumhaz = cumhaz, fit = fit)
    return(breslow.ret)
}
"intcox.derivs" <-
function (data, covar, lambda0u, lambda0v, beta)        # derivatives and likelihood
{
    eps <- 10^(-10)                                     # for numerical stability
    cn <- data$cens
    cens.nr <- (1:length(cn))[cn == 0]
    ncens.nr <- (1:length(cn))[cn == 3]
    e <- exp(beta %*% t(covar))                         # often used terms
    qu <- lambda0u * e                                  #
    qu.n <- qu[cn == 3]                                 #
    qu.c <- qu[cn == 0]                                 #
    qv <- lambda0v * e[cn == 3]                         #
    equ <- exp(-qu)                                     #
    equ.n <- equ[cn == 3]                               #
    equ.c <- equ[cn == 0]                               #
    eqv <- exp(-qv)                                     #
    likeli <- sum(log(equ.n - eqv)) - sum(qu.c)         # Likelihood
    l1u.n <- -qu.n * equ.n/(equ.n - eqv)                # d L/d F0
    l1u.c <- -qu.c                                      #
    l1u <- NULL                                         #
    l1u[ncens.nr] <- l1u.n                              #
    l1u[cens.nr] <- l1u.c                               #
    l1v <- qv * eqv/(equ.n - eqv)                       #
    l1 <- c(l1u, l1v)                                   #
    temp.l2 <- (-qu.n * equ.n + qv * eqv)/(equ.n - eqv) # d L/d beta
    if (sum(3 - cn) == 0) {
        l2 <- as.vector(apply(as.matrix(temp.l2 * covar[cn == 
            3, ]), 2, sum))
    }
    else {
        l2 <- as.vector(apply(as.matrix(temp.l2 * covar[cn == 
            3, ]), 2, sum) - apply(as.matrix(qu.c * covar[cn == 
            0, ]), 2, sum))
    }
    g1u.n <- qu.n * equ.n/(equ.n - eqv) * (1 - qu.n + qu.n * 
        equ.n/(equ.n - eqv))                            # diagonal elements of -d^2 L/d F0^2
    g1u.c <- qu.c
    g1u <- NULL
    g1u[ncens.nr] <- g1u.n
    g1u[cens.nr] <- g1u.c
    g1v <- -qv * eqv/(equ.n - eqv) * (1 - qv - qv * eqv/(equ.n - 
        eqv))
    g1 <- c(g1u, g1v)
    g1 <- g1 + eps                                      # for numerical stability
    temp.g2 <- (-qu.n * equ.n * (1 - qu.n) + qv * eqv * (1 - 
        qv) - (qu.n * equ.n - qv * eqv)^2/(equ.n - eqv))/(equ.n - 
        eqv)                                            # diagonal elements of -d^2 L/d beta^2
    if (sum(3 - cn) == 0) {
        g2 <- -as.vector(apply(as.matrix(temp.g2 * covar[cn == 
            3, ] * covar[cn == 3, ]), 2, sum))
    }
    else {
        g2 <- -as.vector(apply(as.matrix(temp.g2 * covar[cn == 
            3, ] * covar[cn == 3, ]), 2, sum) - apply(as.matrix(qu.c * 
            covar[cn == 0, ] * covar[cn == 0, ]), 2, sum))
    }
    g2 <- g2 + eps                                      # for numerical stability
    derivs.ret <- list(l1 = l1, l2 = l2, g1 = g1, g2 = g2, likeli = likeli)
    return(derivs.ret)
}
"intcox.fit" <-
function (formula, data, covar, eps, itermax)                       # Iterated Convex Minorant Algorithm 
{
    folge <- sort.list(c(data$left, data$right[data$cens == 3]))    # sorting the pooled ends
    rang <- rank(c(data$left, data$right[data$cens == 3]))          # ranks of the pooled ends
    lamb.beta <- intcox.breslow(formula, data, covar)
    lambda0v <- intcox.hazard0(data$right[data$cens == 3], data,
        lamb.beta)                                                  # first estimation of the baseline-hazard for the interval ends
    beta <- lamb.beta$fit                                           # first estimation for beta
    lambda0u <- intcox.hazard0.beg(data$left, data, lamb.beta)      # derived baseline-hazard for the interval beginnings
    lambda0g <- c(lambda0u, lambda0v)
    likeli.vec <- NULL
    null.vector <- rep(0, length(lambda0g))
    itmax <- itermax
    it <- 0                                                         # iteration counter
    abbruch <- 0                                                    # break condition not fulfilled
    while (abbruch == 0) {                                          # until the break condition is fulfilled
        it <- it + 1
        if (it > itmax) 
            abbruch <- 3
        derivs.wert <- intcox.derivs(data, covar, lambda0u, lambda0v, 
            beta)                                                   # first and second derivatives and likelihood
        if (any(derivs.wert$g1 <= 0)) {
            abbruch <- 4
        }
        if (any(!is.finite(derivs.wert$l1))) {
            abbruch <- 4
        }
        if (any(derivs.wert$g2 <= 0)) {
            abbruch <- 4
        }
        if (any(!is.finite(derivs.wert$l2))) {
            abbruch <- 4
        }
        likeli.vec <- c(likeli.vec, derivs.wert$likeli)
        g1.inv <- 1/(derivs.wert$g1)
        g2.inv <- 1/(derivs.wert$g2)
        lambda0g <- c(lambda0u, lambda0v)
        alpha <- 1                                                  # step size
        if (any(!is.finite(lambda0g + alpha * g1.inv * derivs.wert$l1))) {
            abbruch <- 4
        }
        ind <- rep(1, length(lambda0g))
        if (abbruch == 0) {
            repeat {                                                # step size trials
                lambda0g.new <- intcox.pavaC(derivs.wert$g1[folge], 
                  pmax((lambda0g + alpha * g1.inv * derivs.wert$l1)[folge], 
                    null.vector))[rang]                             # PAVA with C
                beta.neu <- beta + alpha * g2.inv * derivs.wert$l2
                lambda0u.neu <- lambda0g.new[1:length(data$left)]
                lambda0v.neu <- lambda0g.new[-(1:length(data$left))]
                likeli.neu <- intcox.derivs(data, covar, lambda0u.neu, 
                  lambda0v.neu, beta.neu)$likeli
                likeli.diff <- likeli.neu - derivs.wert$likeli
                if (is.na(likeli.diff)) {
                  abbruch <- 4
                  break
                }
                if (likeli.diff > 0) {
                  break
                }
                else {
                  alpha <- alpha * 1/2                              # try the half step size
                }
                if (alpha < 0.5^40) {
                  abbruch <- 2
                  break
                }
            }
        }
        if (abbruch == 0) {
            alpha.u <- alpha * 1/2
            lambda0g.u <- intcox.pavaC(derivs.wert$g1[folge], pmax((lambda0g + 
                alpha.u * g1.inv * derivs.wert$l1)[folge], null.vector))[rang]  # PAVA with C
            beta.u <- beta + alpha.u * g2.inv * derivs.wert$l2
            lambda0u.u <- lambda0g.u[1:length(data$left)]
            lambda0v.u <- lambda0g.u[-(1:length(data$left))]
            likeli.u <- intcox.derivs(data, covar, lambda0u.u, lambda0v.u, 
                beta.u)$likeli
            likeli.diff.u <- likeli.u - derivs.wert$likeli
        }
        if (abbruch == 0) {
            alpha.o <- alpha * 3/2
            lambda0g.o <- intcox.pavaC(derivs.wert$g1[folge], pmax((lambda0g + 
                alpha.o * g1.inv * derivs.wert$l1)[folge], null.vector))[rang]  # PAVA with C
            beta.o <- beta + alpha.o * g2.inv * derivs.wert$l2
            lambda0u.o <- lambda0g.o[1:length(data$left)]
            lambda0v.o <- lambda0g.o[-(1:length(data$left))]
            likeli.o <- intcox.derivs(data, covar, lambda0u.o, lambda0v.o, 
                beta.o)$likeli
            likeli.diff.o <- likeli.o - derivs.wert$likeli
        }
        zweig <- 0
        if (abbruch == 0) {
            if (likeli.neu >= likeli.o && likeli.neu >= likeli.u) {
                lambda0u <- lambda0u.neu
                lambda0v <- lambda0v.neu
                beta <- beta.neu
                zweig <- 2
            }
            if (likeli.u >= likeli.o && likeli.u >= likeli.neu) {
                lambda0u <- lambda0u.u
                lambda0v <- lambda0v.u
                likeli.diff <- likeli.diff.u
                beta <- beta.u
                zweig <- 1
            }
            if (likeli.o >= likeli.u && likeli.o >= likeli.neu) {
                lambda0u <- lambda0u.o
                lambda0v <- lambda0v.o
                likeli.diff <- likeli.diff.o
                beta <- beta.o
                zweig <- 3
            }
            likeli.rel <- abs(likeli.diff/derivs.wert$likeli)
            if (is.na(likeli.rel)) {
                abbruch <- 4
            }
            else {
                if (likeli.rel < eps) 
                  abbruch <- 1                                  # if no break
            }
        }
    }
    if (abbruch == 2) 
        likeli.vec <- likeli.vec
    else likeli.vec <- c(likeli.vec, likeli.neu)
    time.point <- c(data$left, data$right[data$cens == 3])
    if (abbruch == 2) 
        lambda0 <- lambda0g
    else lambda0 <- lambda0g.new
    ordnung <- sort.list(time.point)
    time.point <- time.point[ordnung]
    lambda0 <- lambda0[ordnung]
    time.point <- time.point[!duplicated(lambda0)]
    lambda0 <- lambda0[!duplicated(lambda0)]
    intcox.ret <- list(coefficients = beta, lambda0 = lambda0, time.point = time.point, 
        likeli.vec = likeli.vec, iter = it, termination = abbruch)
    return(intcox.ret)
}
"intcox.hazard0" <-
function (t, data, rate)        # baseline hazard
{
    mini <- min(data$left)
    time.intv <- c(mini, data$right[data$cens == 3])
    cumhaz <- c(0, rate$cumhaz)
    hazard0 <- NULL
    for (i in 1:length(t)) {
        hazard0 <- c(hazard0, cumhaz[time.intv <= t[i]][length(cumhaz[time.intv <= 
            t[i]])])
    }
    return(hazard0)
}
"intcox.hazard0.beg" <-
function (t, data, rate)        # baseline hazard for the lower interval ends
{
    mini <- min(data$left)
    time.intv <- c(mini, data$right[data$cens == 3])
    cumhaz <- c(0, rate$cumhaz)
    hazard0 <- NULL
    for (i in 1:length(t)) {
        if (t[i] == mini) {
            hazard0 <- c(hazard0, 0)
        }
        else {
            hazard0 <- c(hazard0, cumhaz[time.intv < t[i]][length(cumhaz[time.intv < 
                t[i]])])
        }
    }
    return(hazard0)
}
"intcox.pavaC" <-
function (w, y) 
{
    n <- length(y)
    index <- rep(0, n)
    weight <- rep(0, n)
    ghat <- rep(0, n)
    .C("pavaC", as.double(w), as.double(y), as.integer(n), as.integer(index), 
        as.double(weight), as.double(ghat), PACKAGE = "intcox" )[[6]]
}
".First.lib" <-
function (lib, pkg) 
{
    library.dynam("intcox", pkg, lib)
}
