.packageName <- "nltm"
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpjnicQE/R.INSTALL10d63af1/nltm/R/anova.nltm.R"
update.nltm.formula <- function(fit,form,pred)
{
  if(pred=="long"){
    if(fit$npred==1){
      fit$call$formula1 <- update.formula(fit$formula,form)
    }else{
      fit$call$formula1 <- update.formula(fit$formula$pred.long,form)
      fit$call$formula2 <- update.formula(fit$formula$pred.short,~1)
    }
  }else{
    if(!is.na(match("formula2",names(fit$call))))
      fit$call$formula2 <- update.formula(fit$formula$pred.short,form)
    else
      fit$call$formula2 <- update.formula(fit$formula$pred.long,form)
  }
  fit
}

anova.nltm <- function(object, ..., test=FALSE)
{
  cat(gettextf("Authors: G. Garibotti, A. Tsodikov\n"))
  dotargs <- list(...)
  named <- if (is.null(names(dotargs))) 
    rep(FALSE, length(dotargs))
  else (names(dotargs) != "")
  if (any(named)) 
    warning(paste("The following arguments to anova.nltm(..)", 
                  "are invalid and dropped:",
                  paste(deparse(dotargs[named]), collapse=", ")))
  dotargs <- dotargs[!named]
  is.nltm <- unlist(lapply(dotargs, function(x) inherits(x,"nltm")))
  dotargs <- dotargs[is.nltm]
  if(length(dotargs)>0)
    return(anova.nltmlist(c(list(object), dotargs), test=test))

  if(object$npred==1){
    if(object$nvar==0)
      stop(paste("\nModel with no covariantes,",
                 "it is not possible to build an anova table.", sep=" "))
  }else{
    if(object$nvar$pred.long==0 & object$nvar$pred.short==0)
      stop(paste("\nModel with no covariantes,",
                 "it is not possible to build an anova table.", sep=" "))
  }
  
  cure <- !is.na(match("cure",names(object$coef)))
  row.name <- character()
  resdev <- resdf <- NULL

  if(object$npred>1){
    if(object$nvar$pred.short>0){
      form <- "~."
      termlist <- attr(object$terms$pred.short,"term.label")
      if(object$nvar$pred.long==0 & !cure)
        termlist <- termlist[-1]
      nvars <- length(termlist)
      if(nvars>0){
        for(i in rev(termlist)){
          form <- paste(form,i,sep="-")
          fit <- update.nltm.formula(object,form,"short")
          fit <- update(fit)
          resdev <- c(resdev, -2*fit$loglik[2])
          resdf <- c(resdf,object$n-sum(!is.na(coef(fit))))
        }
      }
      row.name <- c(paste(attr(object$terms$pred.short,"term.labels"),
                           "(short term)",sep=" "))
    }
    termlist <- attr(object$terms$pred.long,"term.labels")
  }else{
    termlist <- attr(object$terms,"term.labels")
  }
  nvars <- length(termlist)
  form <- "~."
  if(nvars>1){
    for(i in rev(termlist[-1])){
      form <- paste(form,i,sep="-")
      fit <- update.nltm.formula(object,form,"long")
      fit <- update(fit)
      resdev <- c(resdev, -2*fit$loglik[2])
      resdf <- c(resdf,object$n-sum(!is.na(coef(fit))))
    }
    if(cure){
      form <- paste(form,termlist[1],sep="-")
      fit <- update.nltm.formula(object,form,"long")
      fit <- update(fit)
      resdev <- c(resdev, -2*fit$loglik[2])
      resdf <- c(resdf,object$n-sum(!is.na(coef(fit))))
    }
    row.name <- c(paste(termlist, "(long term)", sep=" "),row.name)
  }else{
    if(nvars==1){
      if(cure){
        form <- paste(form,termlist[1],sep="-")
        fit <- update.nltm.formula(object,form,"long")
        fit <- update(fit)
        resdev <- c(resdev, -2*fit$loglik[2])
        resdf <- c(resdf,object$n-sum(!is.na(coef(fit))))
      }
      row.name <- c(paste(termlist, "(long term)", sep=" "),row.name)
    }
  }
  if(cure)
    row.name <- c("cure (long term)",row.name)
  row.name <- c("NULL",row.name)

  resdf <- c(object$n,rev(resdf),object$n-sum(!is.na(coef(object))))
  resdev <- c(-2*object$loglik[1],rev(resdev),-2*object$loglik[2])
  table <- data.frame(c(NA,-diff(resdf)),c(NA,pmax(0,-diff(resdev))),resdf,
                      resdev)
  
  dimnames(table) <- list(row.name,c("Df","Deviance","Resid. Df",
                                      "Resid. Dev"))
  title <- paste("Analysis of Deviance Table\nNLT model: ",
                 object$call$nlt.model, "\nResponse is ",
                 ifelse(object$npred==1,deparse(object$terms[[2]]),
                        deparse(object$terms$pred.long[[2]])),
                 "\nTerms added sequentially (first to last)\n", sep = "")
  dispersion <- 1
  if (test)
    table <- stat.anova(table=table, test="Chisq", scale=dispersion)
  
  structure(table, heading=title, class=c("anova","data.frame"))
}


anova.nltmlist <- function(object, ..., test=FALSE)
{
  responses <- as.character(lapply(object,
                                   function(x) {x$call$formula1[[2]]}))
  same.resp <- responses==responses[1]
  if(!all(same.resp)){
    object <- object[same.resp]
    warning(paste("Models with response", deparse(responses[!same.resp]), 
                  "removed because response differs from model 1."))
  }
  same.n <- sapply(object, function(x) x$n)
  if(any(same.n!=same.n[1]))
    stop("Models were not all fitted to the same size of dataset.")
  nlt.models <- sapply(object, function(x) x$call$nlt.model)
  same.nltm <- nlt.models==nlt.models[1]
  if(!all(same.nltm)){
    object <- object[same.nltm]
    warning(paste("Models with nltm", deparse(responses[!same.nltm]), 
                  "removed because nltm differs from", "model 1"))
  }
  nmodels <- length(object)
  
  if(nmodels==1)
    stop(paste("Only one model left after removing models with different",
               "sample sizes and nltm", sep=" "))
  
  resdf <- as.numeric(lapply(object, function(x) x$n-sum(!is.na(coef(x)))))
  resdev <- as.numeric(lapply(object, function(x) -2*x$loglik[2]))
  table <- data.frame(resdf, resdev, c(NA,-diff(resdf)), c(NA,-diff(resdev)))
  variables1 <- lapply(object, function(x) x$call$formula1)
  dimnames(table) <- list(1:nmodels,
                          c("Resid. Df","Resid. Dev","Df","Deviance"))
  title <- "Analysis of Deviance Table\n"

  if(object[[1]]$npred==1)
    topnote <- paste("Model ", format(1:nmodels), ": ", variables1, 
                     sep="", collapse="\n")
  else{
    num <- character()
    num[seq(1,nmodels*2,2)] <- format(1:nmodels)
    num[seq(2,nmodels*2,2)] <- " "
    variables1[seq(1,nmodels*2,2)] <- variables1
    variables1[seq(2,nmodels*2,2)] <- lapply(object,
                                             function(x) x$call$formula2)
    topnote <- paste(c("Model ","      "), num,
                     c(": Long term:  ","  Short term: "), variables1, sep="",
                     collapse="\n")
  }
  if (test) {
    dispersion <- 1
    table <- stat.anova(table=table, test="Chisq", scale=dispersion)
  }
  structure(table, heading=c(title,topnote), class=c("anova","data.frame"))
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpjnicQE/R.INSTALL10d63af1/nltm/R/nltm.R"
# Predictor 1: long term effect
# Predictor 2: short term effect

nltm <- function(formula1=formula(data), formula2=formula(data),
                 data=parent.frame(), subset, na.action, init=NULL,
                 control, nlt.model=c("PH","PHC","PO","PHPHC","PHPOC",
                            "GFM","PHPO"),
                 model=FALSE, x=FALSE, y=FALSE, verbose=FALSE, ...)
{
  if(sys.parent()==0)
    cat(gettextf("Authors: G. Garibotti, A. Tsodikov\n"))
  if(!nlt.model %in% eval(formals()[["nlt.model"]]))
    stop(gettextf("nlt.model should be one of %s",
                  paste(dQuote(eval(formals()[["nlt.model"]])),collapse=", ")),
         domain=NA)

  call <- match.call()
  m <- match.call(expand=FALSE)
  # this is necessary because otherwise eval(m, parent.frame()) doesn't work
  names(m)[names(m)=="formula1"] <- "formula"  
  temp <- c("","formula","data","subset","na.action")
  m <- m[match(temp, names(m), nomatch=0)]
  Terms1 <- if(missing(data)) terms(formula1)
            else              terms(formula1, data=data)
  m$formula <- Terms1
  m[[1]] <- as.name("model.frame")
  m <- eval(m, parent.frame())
  if(NROW(m)==0)
    stop("No (non-missing) observations")
  
  Y <- model.extract(m, "response")
  if (!inherits(Y, "Surv")) stop("Response must be a survival object")

  attr(Terms1,"intercept") <- 1  # model always has \Lambda_0
  newTerms <- Terms1 
  
  X1 <- model.matrix(newTerms,m)
  assign <- lapply(attrassign(X1,newTerms)[-1],function(x) x-1)
  X1 <- X1[,-1,drop=FALSE]

  npred <- nPredictor(nlt.model)
  if(model & npred>1)
    m1 <- m
  
  if(!missing(formula2)){
    if(npred==1){
      cat(gettextf("\nWarning message:\n"))
      cat(gettextf(paste("Model", nlt.model, sep=" ")))
      cat(gettextf(" has only one predictor however there are two formulas,\nformula2 will not be used.\n\n"))
    }else{
      m <- match.call(expand=FALSE)
      names(m)[names(m)=="formula2"] <- "formula"  
      temp <- c("","formula","data","subset","na.action")
      m <- m[match(temp, names(m), nomatch=0)]
      Terms2 <- if(missing(data)) terms(formula2)
                else              terms(formula2, data=data)
      m$formula <- Terms2
      m[[1]] <- as.name("model.frame")
      m <- eval(m, parent.frame())
      if(NROW(m)==0)
        stop("No (non-missing) observations")

      attr(Terms2,"intercept")<- 1  # model always has \Lambda_0
  
      X2 <- model.matrix(Terms2,m)
      assign <- lapply(attrassign(X2,Terms2)[-1],function(x) x-1)
      X2 <- X2[,-1,drop=FALSE]
    }
  }else{
    if(npred>1){
      X2 <- NULL
      Terms2 <- Terms1
    }
  }
  
  extraArgs <- list(...)
  if (length(extraArgs)) {
    controlargs <- names(formals(nltm.control))  #legal arg names
    indx <- match(names(extraArgs), controlargs, nomatch=0)
    if (any(indx==0))
      stop("Argument ", names(extraArgs)[indx==0], "not matched")
  }
  controls <- nltm.control(...)
  if(!missing(control)) controls[names(control)] <- control
  
  if(verbose!=FALSE){
    res <- .C("openDebug", verbose)
    verbose <- TRUE
  }
  fit <- nltm.fit(X1, X2, Y, nlt.model, init, controls, verbose)
  if(verbose==TRUE) res <- .C("closeDebug")

  if(npred==1){
    fit$formula <- formula(Terms1)
    fit$terms <- Terms1
    fit$nvar <- dim(X1)[2]
  }else{
    fit$formula$pred.long <- formula(Terms1)
    fit$formula$pred.short <- formula(Terms2)
    fit$terms$pred.long <- Terms1
    fit$terms$pred.short <- Terms2
    fit$nvar$pred.long <- dim(X1)[2]
    if(!is.null(X2))
      fit$nvar$pred.short <- dim(X2)[2]
    else
      fit$nvar$pred.short <- dim(X1)[2]
  }
  fit$call <- call
  if(is.null(fit$call$nlt.model))
    fit$call$nlt.model <- "PH"
     
  na.action <- attr(m, "na.action")
  if (length(na.action))
    fit$na.action <- na.action
  fit$npred <- npred
  if(x)
    if(npred==1){
      fit$x <- X1
    }else{
      fit$x$pred.long <- X1
      if(!is.null(X2))
        fit$x$pred.short <- X2
      else
        fit$x$pred.short <- X1
    }
  if(model)
    if(npred==1){
      fit$model <- m
    }else{
      fit$model$pred.long <- m1
      fit$model$pred.short <- m
    }
  if(y)
    fit$y <- Y
  
  class(fit) <- "nltm"
  
  fit
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpjnicQE/R.INSTALL10d63af1/nltm/R/nltm.fit.R"
# y[,1]: times
# y[,2]: status
# x1: regression variables for first predictor (long term)
# x2: regression variables for second predictor (short term)

nltm.fit <- function(x1, x2, y, model, init, control, verbose)
{
  n <-  nrow(y)
  if(is.matrix(x1))
    nvar1 <- ncol(x1)   
  else
    if (length(x1)==0)
      nvar1 <- 0
    else
      nvar1 <- 1

  npred <- nPredictor(model)
  cure <- cureModel(model)

  if(npred>1){
    if(is.null(x2)){
      nbeta <- npred*nvar1+cure
      x2 <- x1
      nvar2 <- ncol(x2)
    }else{
      nvar2 <- ncol(x2)
      nbeta <- nvar1+nvar2+cure
    }
  }else{
    nbeta <- nvar1+cure
    x2 <- 0
    nvar2 <- 0
  }

  # Eliminate censored observations prior to first event time.
  # Individuals with censored time before the first event time are not
  # used in the estimation
  sorted <- order(y[,1])
  nc.small <- length(which(y[,1]<min(y[y[,2]==1,1])))+1
  # All observations are censored
  if(nc.small>n)
    stop("All observations are censored. The model cannot be fitted.\n")
  sorted <- sorted[nc.small:n]
  
  y <- y[sorted,]

  if(nvar1>0){
    names <- dimnames(x1)
    names[[1]] <- names[[1]][sorted]
    x1 <- matrix(x1[sorted,], nrow=length(sorted), dimnames=names)
  }
  if(npred>1 & nvar2>0){
    names <- dimnames(x2)
    names[[1]] <- names[[1]][sorted]
    x2 <- matrix(x2[sorted,], nrow=length(sorted), dimnames=names)
  }
  time <- eventTimes(y)  # time and status correspond to stime and sstat in 
  status <- y[,2]        # coxph.fit
    
  count <- counts(time, status)
  isurv <- initSurvival(count, cure)
  s0 <- isurv$s0
  
  if(!is.null(init)){
    if(length(init)!=nbeta)
      stop("Wrong length for initial values")
  }else{
    init <- rep(0,nbeta)
    if(cure)
      init[nbeta] <- log(-log(isurv$tailDefect))
  }

  if(nbeta>0){
    bound <- boundary(x1, x2, npred, cure, control$bscale)
    
    fit <- optim(par=init, fn=profileLikR, gr=NULL, method="L-BFGS-B",
                 lower=-bound, upper=bound,
                 control=control[-match(c("bscale","s0.tol"),names(control))],
                 hessian=FALSE, x1, x2, status, count, s0, model, cure,
                 control$s0.tol, nvar1, nvar2, as.integer(n-nc.small+1), npred,
                 verbose)
    
    coef <- fit$par
    
    if(npred==1)
      names(coef) <- dimnames(x1)[[2]]
    else
      names(coef) <- c(dimnames(x1)[[2]],dimnames(x2)[[2]])
    if(cure)
      names(coef)[nbeta] <- "cure"

    # Find survival function at betaMLE
    survMLE <- .C("profileLik", coef=as.double(coef), t(x1), t(x2),
                  as.integer(status), count$dd, count$rr, surv=as.double(s0),
                  model, as.integer(cure), control$s0.tol, nvar1, nvar2,
                  nrow(count), as.integer(n-nc.small+1), as.integer(npred),
                  as.integer(verbose), plik=double(1), PACKAGE="nltm")$surv
  
    # Find covariance matrix
    imat1 <- .C("informationMatrix", coef, t(x1), t(x2), as.integer(status),
                count$dd, count$rr, survMLE, model, as.integer(cure), nvar1,
                nvar2, nrow(count), as.integer(n-nc.small+1),
                as.integer(npred), as.integer(verbose),
                imat=double(nbeta*nbeta), PACKAGE="nltm")$imat

    # Invert information matrix
    cov <- solve(matrix(imat1, nbeta, nbeta))

    counts <- fit$count[1]
    names(counts) <- NULL
    res <- list(coefficients=coef, surv=survMLE, var=cov, n=n,
                maxit=control$maxit, counts=counts,
                convergence=fit$convergence, message=fit$message)
  }
  # Find profile likelihood of null model
  init <- rep(0,nbeta)
  lik0 <- .C("profileLik", coef=as.double(init), t(x1), t(x2),
             as.integer(status), count$dd, count$rr, surv=as.double(s0),
             model, as.integer(cure), control$s0.tol, nvar1, nvar2,
             nrow(count), as.integer(n-nc.small+1), as.integer(npred),
             as.integer(verbose), plik=double(1), PACKAGE="nltm")$plik

  if(nbeta>0)
    res$loglik <- c(lik0, fit$value)
  else
    res <- list(loglik=lik0, n=n)
  res
}


# The actual value of the censored times doesn't matter, people with
# censored time in [t_i, t_{i+1}) get t_i
eventTimes <- function(y)
{
  td <- unique(y[y[,2]==1,1])  # td is sorted because y[,1] is
  time <- array(NA, nrow(y)) 
  
  for(i in 1 : length(td))
    time[y[,1]>=td[i]] <- td[i]
  time
}

# count$dd: number of deaths at t_i
# count$rr: number of dead or censored people at t_i
counts <- function(time, status)
{
  count <- data.frame(table(time[status==1]),table(time))
  count <- count[,c(2,4)]
  colnames(count) <- c("dd","rr")
  count
}

reverseCumsum <- function(a)
{
  res <- cumsum(a[length(a):1])[length(a):1]
  res
}


# use the Nelson-Aalen etimator (Klein, Moeschberger, page 86) to
# initialize the survival jumps
# Note: cumprod(res) is the KM in [SP]
initSurvival <- function(count, cure)
{
  if(cure){
    nt <- nrow(count)
    aux <- cumsum(count$dd/reverseCumsum(count$rr))
    res <- aux[nt]-aux
    res <- list(s0=res/c(aux[nt],res[1:(nt-1)]), tailDefect=exp(-aux[nt]))
  }else
    res <- list(s0=exp(-count$dd/reverseCumsum(count$rr)))
  res
}

# compute boundaries for maximization
# Note: this function assumes that either length(bscale)=1 or
# bscale=(bscale_theta, bscale_eta, bscale_cure); that is, long-term
# predictor, short-term predictor, cure term
boundary <- function(x1, x2, npred, cure, bscale)
{
  r <- numeric()
  if(ncol(x1)>0)
    r <- apply(abs(x1), MARGIN=2, max)
  if(npred>1)
    r <- c(r, apply(abs(x2), MARGIN=2, max))
  r <- ifelse(r<1e-10, 1e-10, r)
  
  bound <- bscale/r
  
  if(cure)
    bound <- c(bound, ifelse(length(bscale)>1, bscale[length(bscale)], bscale))
  bound
}


# control parameters for nltm
nltm.control <- function(fnscale=-1, maxit=1000, reltol, factr=1e7, pgtol=0,
                         s0.tol=1e-5, bscale=5)
{
  # control parameters for optimization (see optim help)
  if(fnscale>0) warning("Minimization will take place")
  if(maxit<=0) stop("Invalid value for iterations")
  if(!missing(reltol)) if(reltol<0) stop("Invalid convergence tolerance")
  if(factr<0) stop("Invalid factor convergence tolerance")
  if(pgtol<0) stop("Invalid tolerance on the projected gradient")

  if(s0.tol<0) stop("Invalid convergence tolerance of hazard self-consistency")
  if(sum(bscale<=0)>0) stop("Invalid bound scale")
  
  ifelse(missing(reltol), 
         control <- list(fnscale=fnscale, maxit=maxit, factr=factr,
                         pgtol=pgtol, s0.tol=s0.tol, bscale=bscale),
         control <- list(fnscale=fnscale, maxit=maxit, reltol=reltol,
                         factr=factr, pgtol=pgtol, s0.tol=s0.tol,
                         bscale=bscale))
  control
}


profileLikR <- function(beta, x1, x2, status, count, s0, model, cure, tol,
                        nvar1, nvar2, nobs, npred, verbose)
{
  # need to transpose x otherwise when passed to C it is stored in a vector
  # by columns and I need it by rows because of dmat
  res <- .C("profileLik", coef=as.double(beta), t(x1), t(x2),
            as.integer(status), count$dd, count$rr, surv=as.double(s0), model,
            as.integer(cure), tol, nvar1, nvar2, nrow(count), nobs,
            as.integer(npred), as.integer(verbose), plik=double(1),
            PACKAGE="nltm")$plik 

  res
}

nPredictor <- function(model)
{
  switch(model, PH=1, PHC=1, PO=1, PHPHC=2, PHPOC=2, GFM=2, PHPO=2)
}

cureModel <- function(model)
{
  switch(model, PH=FALSE, PHC=TRUE, PO=FALSE, PHPHC=TRUE, PHPOC=TRUE,
         GFM=FALSE, PHPO=FALSE)
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpjnicQE/R.INSTALL10d63af1/nltm/R/print.nltm.R"
print.nltm <- function(x, digits=max(options()$digits - 4, 3), ...)
{
  cat(gettextf("Authors: G. Garibotti, A. Tsodikov\n"))

  if(!is.null(cl<- x$call)) {
    cat("Call:\n")
    dput(cl)
    cat("\n")
  }

  savedig <- options(digits = digits)
  on.exit(options(savedig))

  if(!is.null(x$coef)){
    coef <- x$coef
    negVar <- which(diag(x$var)<=0)
    if(length(negVar)==0)
      se <- sqrt(diag(x$var))
    else{
      se <- array(NA, dim=length(coef))
      se[-negVar] <- sqrt(diag(x$var)[-negVar])
    }
    
    if(is.null(coef) | is.null(se))
      stop("Input is not valid")

    if(x$convergence==1){
      cat("Iteration limit maxit: ", x$maxit, " has been reached.\n")
      cat("Consider a new initial value or optimization parameters such as",
          "bscale.\n", sep=" ")
      cat("See nltm.control help.\n")
    }else{
      if(x$convergence==51){
        cat("Warning from optimization method: ", x$message, "\n")
        cat("Consider a new initial value or optimization parameters such as",
            "bscale.\n", sep=" ")
        cat("See nltm.control help.\n")
      }else{
        if(x$convergence==52){
          cat("Error in optimization method: ", x$message, "\n")
          cat("Consider a new initial value or optimization parameters such",
              "as bscale.\n", sep=" ")
          cat("See nltm.control help.\n")
        }
      }
    }

    if(length(negVar)>0){
      cat(gettextf("\nWarning message:\n"))
      cat(gettextf("Problem with covariance matrix of coefficients.\n"))
      cat(gettextf("Some diagonal terms are negative.\n"))
      cat(gettextf("p-values and confidence intervals are unreliable.\n\n\n"))
    }
    
    cat("Non Linear Transformation Model: ", x$call$nlt.model,
        ", fit by maximum likelihood\n\n",sep="")
    tmp <- cbind(coef, exp(coef), se, coef/se,
                 signif(1 - pchisq((coef/ se)^2, 1), digits-1))
    dimnames(tmp) <- list(names(coef),
                          c("coef","exp(coef)","se(coef)","z","p"))
    
    if(x$npred==1){
      cat("\n")
      prmatrix(tmp)
    }else{
      ind1 <- numeric()
      if(x$nvar$pred.long>0)
        ind1 <- 1:x$nvar$pred.long
      if(!is.na(match("cure",names(coef))))
        ind1 <- c(ind1,length(coef))
      if(length(ind1)>0){
        cat("Long term predictor\n")
        prmatrix(matrix(tmp[ind1,],nrow=length(ind1)),
                 rowlab=rownames(tmp)[ind1],collab=colnames(tmp))
      }
      ind1 <- setdiff(1:length(coef),ind1)
      if(length(ind1>0)){
        cat("\nShort term predictor\n")
        prmatrix(matrix(tmp[ind1,],nrow=length(ind1)),
                 rowlab=rownames(tmp)[ind1],collab=colnames(tmp))
      }
    }
    logtest <- -2 * (x$loglik[1] - x$loglik[2])
    if (is.null(x$df)) df <- sum(!is.na(coef))
    else  df <- round(sum(x$df),2)
    cat("\n")
    cat("Likelihood ratio test=", format(round(logtest, digits)), " on ",
        df, " df,", " p=", format(1 - pchisq(logtest, df)),  sep="")
    cat("\n")
    omit <- x$na.action
    if(length(omit))
      cat("n=", x$n, " (", naprint(omit), ")\n", sep="")
    else cat("n=", x$n, "\n", sep="")
  }else{
    cat("Null model: ",x$call$nlt.model,"\n",sep="")
    cat("Likelihood=", format(round(x$loglik[1], digits)), " n=", x$n, "\n",
        sep="")
  }
  invisible()
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpjnicQE/R.INSTALL10d63af1/nltm/R/summary.nltm.R"
summary.nltm <- function(object,  coef = TRUE, conf.int = 0.95, 
                         digits = max(options()$digits - 4, 3),...)
{
  if(!is.null(object$coef)){
    beta <- object$coef
    nabeta <- !(is.na(beta))          # non-missing coefs
    beta2 <- beta[nabeta]
    if(is.null(beta) | is.null(object$var))
      stop("Input is not valid")
    
    negVar <- which(diag(object$var)<=0)
    if(length(negVar)==0)
      se <- sqrt(diag(object$var))
    else{
      se <- array(NA, dim=length(beta))
      se[-negVar] <- sqrt(diag(object$var)[-negVar])
    }
    rval <- list(call=object$call, convergence=object$convergence,
                 na.action=object$na.action, n=object$n)
    
    if(coef){
      rval$coef <- cbind(beta, exp(beta), se, beta/se,
                         signif(1 - pchisq((beta/se)^2, 1), digits -1))
      dimnames(rval$coef) <- list(names(beta),
                                  c("coef", "exp(coef)", "se(coef)", "z", "p"))
    }
    if(conf.int){
      z <- qnorm((1 + conf.int)/2, 0, 1)
      tmp <- cbind(exp(beta), exp(-beta), exp(beta-z*se), exp(beta+z*se))
      dimnames(tmp) <- list(names(beta),
                            c("exp(coef)", "exp(-coef)",
                              paste("lower .",round(100*conf.int,2),sep = ""),
                              paste("upper .",round(100*conf.int,2),sep = "")))
      rval$conf.int <- tmp
    }
    df <- length(beta2)
    logtest <- -2 * (object$loglik[1] - object$loglik[2])
    rval$logtest <- c(test=logtest, df=df, pvalue=1 - pchisq(logtest, df))
  }else{
    rval <- list(call=object$call, n=object$n, loglik=object$loglik[1])
  }
  rval$npred <- object$npred
  rval$nvar <- object$nvar
  class(rval) <- "summary.nltm"
  rval
}


print.summary.nltm <- function(x, digits = max(options()$digits - 4, 3), ...)
{
  cat(gettextf("Authors: G. Garibotti, A. Tsodikov\n"))
  if(!is.null(x$call)) {
    cat("Call:\n")
    dput(x$call)
    cat("\n")
  }

  if(!is.null(x$convergence)){
    if(x$convergence==1){
      cat("Iteration limit maxit: ", x$maxit, " has been reached.\n")
      cat("Consider a new initial value or optimization parameters such",
          "as bscale.\n",sep=" ")
      cat("See nltm.control help\n")
    }else{
      if(x$convergence==51){
        cat("Warning from optimization method: ", x$message, "\n")
        cat("Consider a new initial value or optimization parameters such",
            "as bscale.\n", sep=" ")
        cat("See nltm.control help\n")
      }else{
        if(x$convergence==52){
          cat("Error in optimization method: ", x$message, "\n")
          cat("Consider a new initial value or optimization parameters such",
              "as bscale.\n", sep=" ")
          cat("See nltm.control help\n")
        }
      }
    }

    negVar <- which(diag(x$var)<=0)
    if(length(negVar)>0){
      cat(gettextf("\nWarning message:\n"))
      cat(gettextf("Problem with covariance matrix of coefficients.\n"))
      cat(gettextf("Some diagonal terms are negative.\n"))
      cat(gettextf("p-values and confidence intervals are unreliable.\n\n\n"))
    }
  
    savedig <- options(digits = digits)
    on.exit(options(savedig))
  
    cat("Non Linear Transformation Model: ", x$call$nlt.model,
        ", fit by maximum likelihood\n\n", sep="")

    if(x$npred==1){
      if(!is.null(x$coef)) {
        prmatrix(x$coef)
      }
      if(!is.null(x$conf.int)) {
        cat("\n")
        prmatrix(x$conf.int)
      }
    }else{
      ind1 <- numeric()
      if(x$nvar$pred.long>0)
        ind1 <- 1:x$nvar$pred.long
      if(!is.na(match("cure",dimnames(x$coef)[[1]])))
        ind1 <- c(ind1,dim(x$coef)[1])    
      if(length(ind1)>0){
        cat("Long term predictor\n")
        if(!is.null(x$coef[ind1,])) 
          prmatrix(matrix(x$coef[ind1,],nrow=length(ind1)),
                   rowlab=rownames(x$coef)[ind1],collab=colnames(x$coef))
        if(!is.null(x$conf.int[ind1,])) {
          cat("\n")
          prmatrix(matrix(x$conf.int[ind1,],nrow=length(ind1)),
                   rowlab=rownames(x$conf.int)[ind1],
                   collab=colnames(x$conf.int))
        }
      }
      ind1 <- setdiff(1:dim(x$coef)[1],ind1)    
      if(length(ind1>0)){
        cat("\nShort term predictor\n")
        if(!is.null(x$coef[ind1,]))
          prmatrix(matrix(x$coef[ind1,],nrow=length(ind1)),
                   rowlab=rownames(x$coef)[ind1],collab=colnames(x$coef))
        if(!is.null(x$conf.int[ind1,])) {
          cat("\n")
          prmatrix(matrix(x$conf.int[ind1,],nrow=length(ind1)),
                   rowlab=rownames(x$conf.int)[ind1],
                   collab=colnames(x$conf.int))
        }
      }
    }
    cat("\n")
    cat("Likelihood ratio test=", format(round(x$logtest["test"], 2)),
        " on ", x$logtest["df"], " df, p=", format(x$logtest["pvalue"]),
        sep="")
    cat("\n")
    
    omit <- x$na.action
    if(length(omit))
      cat("n=", x$n, " (", naprint(omit), ")\n", sep="")
    else cat("n=", x$n, "\n",sep="")
  }else{
    cat("Null model: ",x$call$nlt.model,"\n",sep="")
    cat("Likelihood=", format(round(x$loglik[1], digits)), " n=", x$n, "\n",
        sep="")
  }
  invisible()
}
