setClass("tvcCoxph", contains="mle2")
## setClass("tvcCoxphList", contains="list")

cox.tvc <- function(obj,var=NULL,method="logt") {
    stopifnot(attr(obj$y,"type") == "right")
    if (is.null(var))
        return(as(lapply(attr(obj$terms,"term.labels"),
                      function(name) cox.tvc(obj, var=name,method=method)),
                  "tvcCoxList"))
    method <- match.arg(method)
    y <- as.matrix(obj$y)
    time <- y[,1]
    status <- y[,2]
    X <- as.matrix(model.matrix(obj))
    index <- order(time,-status)
    X <- X[index, , drop = FALSE]
    time <- time[index]
    status <- status[index]
    k <- match(var,attr(obj$terms,"term.labels"))
    beta <- c(coef(obj),0)
    names(beta) <- c(names(coef(obj)),sprintf("%s:log(t)",var))
    minuslogl <- function(beta) -.Call("test_cox_tvc2",
                                   list(time=time,event=status,X=X,beta=beta,k=k-1),PACKAGE="rstpm2")
    gr <- function(beta) -.Call("test_cox_tvc2_grad",
                                list(time=time,event=status,X=X,beta=beta,k=k-1),PACKAGE="rstpm2")
    parnames(minuslogl) <- parnames(gr) <- names(beta)
    fit <- mle2(start=beta,
                 minuslogl = minuslogl,
                 gr = gr,
                 method="BFGS")
    fit@data <- list(object=obj,k=k,var=var)
    ## return the mle2 object
    as(fit,"tvcCoxph")
}

setMethod("show", "tvcCoxph", function(object) print(summary(object)))

setMethod("plot", signature(x="tvcCoxph", y="missing"),
          function(x,y,
                   add=FALSE,rug=!add,
                   type="l",lty=c(1,2,2),col=1,ylab="Effect",log="y",...) {
              obj <- x@data$object
              y <- as.matrix(obj$y)
              time <- y[,1]
              status <- y[,2]
              timek <- seq(min(time[status==1]), max(time[status==1]), length=301)
              Xk <- cbind(1,log(timek))
              k <- x@data$k
              index <- c(k,length(coef(x)))
              betak <- coef(x)[index]
              fitted <- as.vector(Xk %*% betak)
              gd <- t(Xk)
              se.fit <- sqrt(colSums(gd * (vcov(x)[index,index] %*% gd)))
              if (add) {
                  matlines(timek,exp(fitted+cbind(0,-1.96*se.fit,1.96*se.fit)),lty=lty,col=col,...)
              } else {
                  matplot(timek,exp(fitted+cbind(0,-1.96*se.fit,1.96*se.fit)),type=type,lty=lty,col=col,ylab=ylab,log=log,...)
              }
              abline(h=exp(coef(obj)[k]),lty=3)
              if (rug)
                  rug(time[status==1])
          })