R/crrstep.r

In crrstep: Stepwise Covariate Selection for the Fine & Gray Competing Risks Regression Model

### crrstep function with custom class output
### Updated 12/31/25

crrstep <- function(formula, scope.min = ~1, etype, ..., subset, data, 
                    direction = c("backward", "forward"), 
                    criterion = c("AIC", "BICcr", "BIC"), 
                    crr.object = FALSE, trace = TRUE, steps = 100) {
  
  # Improved version of crrstep:
  # 1. Properly handles factors, interactions, and polynomial terms (e.g., I(x^2):factor)
  # 2. Fixes bug where null model returned wrong log-likelihood
  # 3. Avoids formula deprecation warnings
  # 4. Returns custom crrstep class with print and summary methods
  
  # Helper function to safely convert to formula
  safe_as_formula <- function(x) {
    if (inherits(x, "formula")) {
      return(x)
    }
    if (inherits(x, "terms")) {
      # Convert terms to formula safely
      f <- reformulate(attr(x, "term.labels"))
      return(f)
    }
    if (is.character(x)) {
      if (length(x) > 1) {
        x <- paste(x, collapse = "")
      }
      return(as.formula(x))
    }
    return(as.formula(x))
  }
  
  # Helper function to drop a term from formula safely
  safe_drop_term <- function(current_formula, term_to_drop) {
    current_terms <- attr(terms(current_formula), "term.labels")
    term_idx <- which(current_terms == term_to_drop)
    
    if (length(term_idx) == 0) {
      return(current_formula)
    }
    
    remaining_terms <- current_terms[-term_idx]
    
    if (length(remaining_terms) == 0) {
      return(as.formula("~ 1"))
    } else {
      return(as.formula(paste("~", paste(remaining_terms, collapse = " + "))))
    }
  }
  
  # Helper function to add a term to formula safely
  safe_add_term <- function(current_formula, term_to_add) {
    current_terms <- attr(terms(current_formula), "term.labels")
    
    if (length(current_terms) == 0) {
      return(as.formula(paste("~", term_to_add)))
    } else {
      new_terms <- c(current_terms, term_to_add)
      return(as.formula(paste("~", paste(new_terms, collapse = " + "))))
    }
  }
  
  # Helper function to extract AIC
  crr.extractAIC <- function(object, cov1, k = 2) {
    round(-2 * object$loglik + k * ncol(cov1), 2)
  }
  
  # Improved function to make covariate matrix with proper tracking
  crr.makecov1 <- function(formula, data) {
    # Extract just the RHS if it's a two-sided formula
    if (length(formula) == 3) {
      formula <- formula[-2]  # Remove LHS, keep only RHS
    }
    
    # Ensure formula is a proper formula object
    formula <- safe_as_formula(formula)
    
    mf <- model.frame(formula, data = as.data.frame(data), na.action = na.pass)
    mt <- attr(mf, "terms")
    cov1 <- model.matrix(mt, data = mf)
    
    # Store the assign attribute before removing intercept
    assign_attr <- attr(cov1, "assign")
    term_labels <- attr(mt, "term.labels")
    
    # Remove intercept column
    if ("(Intercept)" %in% colnames(cov1)) {
      cov1 <- cov1[, -1, drop = FALSE]
      assign_attr <- assign_attr[-1]
    }
    
    # Store metadata as attributes
    attr(cov1, "assign") <- assign_attr
    attr(cov1, "term.labels") <- term_labels
    
    return(cov1)
  }
  
  # Function to get column indices for a term
  get_term_cols <- function(term_name, cov_matrix) {
    assign_vec <- attr(cov_matrix, "assign")
    term_labels <- attr(cov_matrix, "term.labels")
    
    # Find which term number this is
    term_idx <- which(term_labels == term_name)
    
    if (length(term_idx) == 0) {
      return(integer(0))
    }
    
    # Return all columns assigned to this term
    which(assign_vec == term_idx)
  }
  
  # Function to remove columns for a specific term
  remove_term_cols <- function(cov_matrix, term_name) {
    cols_to_remove <- get_term_cols(term_name, cov_matrix)
    
    if (length(cols_to_remove) == 0) {
      return(cov_matrix)
    }
    
    if (length(cols_to_remove) == ncol(cov_matrix)) {
      # Return empty matrix but maintain structure
      result <- cov_matrix[, integer(0), drop = FALSE]
      attr(result, "assign") <- integer(0)
      attr(result, "term.labels") <- character(0)
      return(result)
    }
    
    # Remove columns
    result <- cov_matrix[, -cols_to_remove, drop = FALSE]
    
    # Update assign attribute
    old_assign <- attr(cov_matrix, "assign")
    new_assign <- old_assign[-cols_to_remove]
    
    # Renumber assignments
    term_labels <- attr(cov_matrix, "term.labels")
    term_idx <- which(term_labels == term_name)
    new_assign[new_assign > term_idx] <- new_assign[new_assign > term_idx] - 1
    
    # Update term labels
    new_term_labels <- term_labels[-term_idx]
    
    attr(result, "assign") <- new_assign
    attr(result, "term.labels") <- new_term_labels
    
    return(result)
  }
  
  # Function to add columns for a specific term
  add_term_cols <- function(base_matrix, full_matrix, term_name) {
    cols_to_add <- get_term_cols(term_name, full_matrix)
    
    if (length(cols_to_add) == 0) {
      return(base_matrix)
    }
    
    # Extract columns to add
    new_cols <- full_matrix[, cols_to_add, drop = FALSE]
    
    # Combine matrices
    if (ncol(base_matrix) == 0) {
      result <- new_cols
      attr(result, "assign") <- rep(1, ncol(new_cols))
      attr(result, "term.labels") <- term_name
    } else {
      result <- cbind(base_matrix, new_cols)
      
      # Update assign attribute
      base_assign <- attr(base_matrix, "assign")
      n_terms_base <- length(attr(base_matrix, "term.labels"))
      new_assign <- c(base_assign, rep(n_terms_base + 1, length(cols_to_add)))
      
      # Update term labels
      base_labels <- attr(base_matrix, "term.labels")
      new_term_labels <- c(base_labels, term_name)
      
      attr(result, "assign") <- new_assign
      attr(result, "term.labels") <- new_term_labels
    }
    
    return(result)
  }
  
  # Helper to format formula for printing
  format_formula <- function(f) {
    paste(deparse(f, width.cutoff = 500), collapse = "")
  }
  
  # Drop term function
  crr.dropterm <- function(scope.max, scope.min, ftime, fstatus, data, ...) {
    Terms <- terms(scope.max)
    scope <- attr(Terms, "term.labels")
    
    if (missing(scope.min)) 
      scope.min <- as.formula(~1)
    else if (is.null(scope.min)) 
      scope.min <- as.formula(~1)
    
    Terms1 <- terms(scope.min)
    scope1 <- attr(Terms1, "term.labels")
    
    # Get covariate matrices
    cov1 <- crr.makecov1(formula = scope.max, data = data)
    cov11 <- crr.makecov1(formula = scope.min, data = data)
    
    # Determine which terms can be dropped
    if (length(scope1) == 0) {
      terms_to_test <- scope
      current_cov <- cov1
    } else {
      terms_to_test <- scope[!(scope %in% scope1)]
      current_cov <- cov1
    }
    
    ns <- length(terms_to_test)
    
    # Initialize results matrix
    ans <- matrix(nrow = ns + 1, ncol = 1, 
                  dimnames = list(c("<none>", terms_to_test), "AIC"))
    
    # Fit current model
    if (ncol(current_cov) == 0) {
      init <- crr(ftime, fstatus, cov1 = matrix(0, nrow = length(ftime), ncol = 1), 
                  variance = FALSE, ...)
      ans[1, ] <- round(-2 * init$loglik.null, 2)
    } else {
      init <- crr(ftime, fstatus, current_cov, variance = FALSE, ...)
      ans[1, ] <- crr.extractAIC(init, current_cov, k)
    }
    
    # Test dropping each term
    for (i in seq_along(terms_to_test)) {
      term_to_drop <- terms_to_test[i]
      
      # Remove this term
      test_cov <- remove_term_cols(current_cov, term_to_drop)
      
      if (ncol(test_cov) == 0) {
        # Null model
        ans[i + 1, ] <- round(-2 * init$loglik.null, 2)
      } else {
        object2 <- crr(ftime, fstatus, test_cov, variance = FALSE, ...)
        ans[i + 1, ] <- crr.extractAIC(object2, test_cov, k)
      }
    }
    
    aod <- data.frame(ans[, 1])
    colnames(aod) <- criterion
    
    # Create heading with safely formatted formula
    head <- c("Single term deletions", "\nModel:", format_formula(scope.max))
    class(aod) <- "data.frame"
    attr(aod, "heading") <- head
    aod
  }
  
  # Add term function
  crr.addterm <- function(scope.max, scope.min, ftime, fstatus, data, ...) {
    Terms <- terms(scope.max)
    scope <- attr(Terms, "term.labels")
    
    if (missing(scope.min)) 
      scope.min <- as.formula(~1)
    else if (is.null(scope.min)) 
      scope.min <- as.formula(~1)
    
    Terms1 <- terms(scope.min)
    scope1 <- attr(Terms1, "term.labels")
    
    # Get covariate matrices
    cov_full <- crr.makecov1(scope.max, data)
    cov_current <- crr.makecov1(scope.min, data)
    
    # Determine which terms can be added
    if (length(scope1) == 0) {
      terms_to_test <- scope
    } else {
      terms_to_test <- scope[!(scope %in% scope1)]
    }
    
    ns <- length(terms_to_test)
    
    # Initialize results matrix
    ans <- matrix(nrow = ns + 1, ncol = 1, 
                  dimnames = list(c("<none>", terms_to_test), "AIC"))
    
    # Fit current model
    if (ncol(cov_current) == 0) {
      init <- crr(ftime, fstatus, cov1 = matrix(0, nrow = length(ftime), ncol = 1),
                  variance = FALSE, ...)
      ans[1, ] <- round(-2 * init$loglik.null, 2)
    } else {
      init <- crr(ftime, fstatus, cov_current, variance = FALSE, ...)
      ans[1, ] <- crr.extractAIC(init, cov_current, k)
    }
    
    # Test adding each term
    for (i in seq_along(terms_to_test)) {
      term_to_add <- terms_to_test[i]
      
      # Add this term
      test_cov <- add_term_cols(cov_current, cov_full, term_to_add)
      
      object2 <- crr(ftime, fstatus, test_cov, variance = FALSE, ...)
      ans[i + 1, ] <- crr.extractAIC(object2, test_cov, k)
    }
    
    aod <- data.frame(ans[, 1])
    colnames(aod) <- criterion
    
    # Create heading with safely formatted formula
    head <- c("Single term additions", "\nModel:", format_formula(scope.min))
    class(aod) <- "data.frame"
    attr(aod, "heading") <- head
    aod
  }
  
  cut.string <- function(string) {
    if (length(string) > 1) 
      string[-1] <- paste("\n", string[-1], sep = "")
    string
  }
  
  # Main function body
  
  stopifnot(identical(formula[[1]], as.name("~")))
  
  # Extract LHS and RHS - handle long formulas properly
  if (length(formula) == 3) {
    lhs <- formula[[2]]
    rhs_deparsed <- paste(deparse(formula[[3]], width.cutoff = 500), collapse = "")
    scope.max <- as.formula(paste("~", rhs_deparsed))
  } else {
    lhs <- NULL
    scope.max <- formula
  }
  
  stopifnot(identical(scope.min[[1]], as.name("~")))
  scope.min <- if (length(scope.min) == 3) {
    rhs_min <- paste(deparse(scope.min[[3]], width.cutoff = 500), collapse = "")
    as.formula(paste("~", rhs_min))
  } else {
    as.formula(scope.min)
  }
  
  fvars <- sort(attr(terms(scope.max), "term.labels"))
  svars <- sort(attr(terms(scope.min), "term.labels"))
  
  if (length(svars) > length(fvars)) 
    stop("`scope.min' contains more variables than `formula'")
  if (!all(svars %in% fvars)) 
    stop("`scope.min' contains variables not in `formula'")
  
  # Handle subset - matching original behavior
  if (!missing(subset)) {
    data <- data[subset, ]
  }
  
  # Extract time and status
  # Time from LHS of formula
  if (!is.null(lhs)) {
    ftime <- eval(lhs, envir = data)
  } else {
    stop("Formula must have a left-hand side specifying the time variable")
  }
  
  # Status from etype argument - matching original behavior
  fstatus <- eval(substitute(etype), envir = data)
  
  # Get model frame for covariates to check for missing
  Xmat <- model.frame(scope.max, data, na.action = na.pass)
  
  # Combine time, status, and covariates to find complete cases
  all_data <- cbind(time = ftime, status = fstatus, Xmat)
  nomiss <- complete.cases(all_data)
  
  if (sum(!nomiss) > 0 && trace) {
    cat(sum(!nomiss), "cases omitted due to missing values\n")
  }
  
  # Subset to complete cases
  data <- data[nomiss, ]
  ftime <- ftime[nomiss]
  fstatus <- fstatus[nomiss]
  
  # Store sample size information
  n_obs <- length(ftime)
  n_events <- sum(fstatus == 1)
  
  # Set criterion
  criterion <- match.arg(criterion, choices = c("AIC", "BICcr", "BIC"))
  if (criterion == "AIC") 
    k <- 2
  if (criterion == "BICcr") 
    k <- log(sum(fstatus == 1))
  if (criterion == "BIC") 
    k <- log(nrow(data))
  
  rm(Xmat)
  invisible(gc())
  
  call <- match.call()
  if (trace) 
    print(call)
  
  direction <- match.arg(direction)
  backward <- direction == "backward"
  forward <- direction == "forward"
  
  # Initialize for backward selection
  if (backward) {
    current_formula <- scope.max
    current_cov <- crr.makecov1(formula = current_formula, data = data)
    
    # Check if already at minimum
    if (identical(fvars, svars)) {
      object <- crr(ftime, fstatus, cov1 = current_cov, variance = TRUE, ...)
      fit <- object
      
      std.error <- sqrt(diag(fit$var))
      outmat <- cbind(signif(fit$coef, 3), signif(std.error, 3), 
                      signif(abs(fit$coef)/std.error, 3))
      colnames(outmat) <- c("estimate", "std.error", "t-stat")
      rownames(outmat) <- colnames(current_cov)
      
      if (crr.object) 
        return(fit)
      else return(list(coefficients = outmat, 
                       log.likelihood = round(fit$loglik, 2)))
    }
    
    object <- crr(ftime, fstatus, cov1 = current_cov, variance = TRUE, ...)
    fit <- object
  }
  
  # Initialize for forward selection
  if (forward) {
    current_formula <- scope.min
    vars.in <- attr(terms(scope.min), "term.labels")
    
    if (length(vars.in) == 0) {
      current_cov <- crr.makecov1(scope.min, data = data)
      if (ncol(current_cov) == 0) {
        # Start from null model
        object <- crr(ftime, fstatus, 
                      cov1 = matrix(0, nrow = length(ftime), ncol = 1),
                      variance = FALSE, ...)
      } else {
        object <- crr(ftime, fstatus, cov1 = current_cov, variance = TRUE, ...)
      }
    } else {
      current_cov <- crr.makecov1(scope.min, data = data)
      object <- crr(ftime, fstatus, cov1 = current_cov, variance = TRUE, ...)
    }
    fit <- object
  }
  
  # Stepwise loop
  while (steps > 0) {
    steps <- steps - 1
    aod <- NULL
    change <- NULL
    
    if (backward) {
      # Check if at minimum scope
      current_terms <- attr(terms(current_formula), "term.labels")
      min_terms <- attr(terms(scope.min), "term.labels")
      
      if (setequal(current_terms, min_terms)) {
        break
      }
      
      if (ncol(current_cov) == 0) 
        break
      
      # Get drop term table
      aod <- crr.dropterm(scope.max = current_formula, scope.min = scope.min, 
                          ftime = ftime, fstatus = fstatus, data = data, ...)
      
      rn <- rownames(aod)
      row.names(aod) <- c(rn[1], paste("-", rn[-1], sep = ""))
      
      if (is.null(aod) || ncol(aod) == 0) 
        break
      
      o <- order(aod[, 1])
      aod.o <- aod[o, , drop = FALSE]
      
      if (trace) {
        print(aod.o)
        utils::flush.console()
      }
      
      # Check if best is to keep current model
      if (o[1] == 1) {
        break
      }
      
      # Drop the best term
      term_to_drop <- rownames(aod)[o[1]]
      term_to_drop <- sub("^-", "", term_to_drop)
      
      if (trace) {
        cat("Dropping:", term_to_drop, "\n")
      }
      
      # Update formula by dropping term using safe function
      current_formula <- safe_drop_term(current_formula, term_to_drop)
      current_cov <- crr.makecov1(formula = current_formula, data = data)
      
      if (ncol(current_cov) == 0) {
        fit <- crr(ftime, fstatus, 
                   cov1 = matrix(0, nrow = length(ftime), ncol = 1),
                   variance = FALSE, ...)
        
        if (trace) {
          cat("\nStep: ", criterion, "= ", 
              format(round(-2 * fit$loglik.null, 2)), "\n")
          cat("Null model\n\n")
          utils::flush.console()
        }
        break
      } else {
        fit <- crr(ftime, fstatus, cov1 = current_cov, variance = TRUE, ...)
        
        if (trace) {
          cat("\nStep: ", criterion, "= ", 
              format(round(crr.extractAIC(fit, current_cov, k), 2)), "\n")
          cat(format_formula(current_formula), "\n\n")
          utils::flush.console()
        }
      }
    }
    
    if (forward) {
      # Get add term table
      aod <- crr.addterm(scope.max, current_formula, ftime, fstatus, data, ...)
      
      # Filter out already included terms
      rn <- rownames(aod)
      invar <- if (is.null(vars.in)) {
        rep(FALSE, length(rn))
      } else {
        (rn %in% vars.in) | (rn == "<none>")
      }
      
      # Keep <none> row
      invar[1] <- FALSE
      aod_filtered <- aod[!invar, , drop = FALSE]
      
      if (nrow(aod_filtered) == 1) {
        # Only <none> left
        break
      }
      
      rn <- rownames(aod_filtered)
      row.names(aod_filtered) <- c(rn[1], paste("+", rn[-1], sep = ""))
      
      if (is.null(aod_filtered) || nrow(aod_filtered) <= 1) 
        break
      
      o <- order(aod_filtered[, 1])
      aod.o <- aod_filtered[o, , drop = FALSE]
      
      if (trace) {
        print(aod.o)
        utils::flush.console()
      }
      
      # Check if best is to keep current model
      if (o[1] == 1) {
        break
      }
      
      # Add the best term
      term_to_add <- rownames(aod_filtered)[o[1]]
      term_to_add <- sub("^\\+", "", term_to_add)
      
      if (trace) {
        cat("Adding:", term_to_add, "\n")
      }
      
      # Update formula by adding term using safe function
      current_formula <- safe_add_term(current_formula, term_to_add)
      current_cov <- crr.makecov1(formula = current_formula, data = data)
      fit <- crr(ftime, fstatus, cov1 = current_cov, variance = TRUE, ...)
      
      vars.in <- c(vars.in, term_to_add)
      
      if (trace) {
        cat("\nStep: ", criterion, "= ", 
            format(round(crr.extractAIC(fit, current_cov, k), 2)), "\n")
        cat(format_formula(current_formula), "\n\n")
        utils::flush.console()
      }
      
      # Check if we've reached the maximum scope
      max_terms <- attr(terms(scope.max), "term.labels")
      current_terms <- attr(terms(current_formula), "term.labels")
      if (setequal(current_terms, max_terms)) {
        break
      }
    }
  }
  
  # Prepare final output
  is_null_model <- ncol(current_cov) == 0
  
  # Set default confidence level
  conf.level <- 0.95
  
  if (is_null_model) {
    # Null model output
    loglik <- fit$loglik.null
    loglik_null <- fit$loglik.null
    criterion_value <- round(-2 * loglik, 2)
    coefficients <- NULL
    std.error <- NULL
    zvalue <- NULL
    pvalue <- NULL
    conf.int <- NULL
    converged <- TRUE
    
    if (trace) {
      cat("Final model: Null model\n")
      cat("Log-likelihood:", round(loglik, 2), "\n")
    }
    
  } else {
    # Non-null model output
    loglik <- fit$loglik
    loglik_null <- fit$loglik.null
    criterion_value <- crr.extractAIC(fit, current_cov, k)
    
    coefficients <- fit$coef
    std.error <- sqrt(diag(fit$var))
    zvalue <- coefficients / std.error
    pvalue <- 2 * (1 - pnorm(abs(zvalue)))
    
    # Confidence intervals (95%)
    z_crit <- qnorm((1 + conf.level) / 2)
    conf.int <- cbind(
      lower = coefficients - z_crit * std.error,
      upper = coefficients + z_crit * std.error
    )
    rownames(conf.int) <- names(coefficients)
    
    converged <- fit$converged
    
    if (trace) {
      cat("Final model:\n")
      outmat <- cbind(
        estimate = signif(coefficients, 3),
        std.error = signif(std.error, 3),
        z = signif(zvalue, 3),
        `p-value` = signif(pvalue, 3)
      )
      rownames(outmat) <- colnames(current_cov)
      print(outmat)
      cat("\nLog-likelihood:", round(loglik, 2), "\n")
    }
  }
  
  # Return crr object if requested
  if (crr.object) {
    return(fit)
  }
  
  # Build result object
  result <- list(
    call = call,
    formula = current_formula,
    direction = direction,
    criterion = criterion,
    criterion_value = criterion_value,
    coefficients = coefficients,
    std.error = std.error,
    zvalue = zvalue,
    pvalue = pvalue,
    conf.int = conf.int,
    conf.level = conf.level,
    loglik = loglik,
    loglik_null = loglik_null,
    n = n_obs,
    nevent = n_events,
    converged = converged,
    is_null_model = is_null_model,
    crr_fit = fit
  )
  
  class(result) <- "crrstep"
  return(result)
}


#' @export
print.crrstep <- function(x, digits = max(3L, getOption("digits") - 3L), ...) {
  
  cat("\nCall:\n")
  print(x$call)
  
  cat("\n")
  cat("Direction:", x$direction, "\n")
  cat("Criterion:", x$criterion, "\n")
  
  cat("\nSelected model:\n")
  if (x$is_null_model) {
    cat("  ~ 1 (Null model)\n")
  } else {
    cat(" ", deparse(x$formula, width.cutoff = 500), "\n")
  }
  
  cat("\n", x$criterion, ": ", format(round(x$criterion_value, 2), nsmall = 2), "\n", sep = "")
  cat("Log-likelihood:", format(round(x$loglik, 2), nsmall = 2), "\n")
  
  cat("\nN =", x$n, " Events =", x$nevent, "\n")
  
  if (!x$is_null_model) {
    cat("\nCoefficients:\n")
    coef_table <- cbind(
      coef = x$coefficients,
      `exp(coef)` = exp(x$coefficients),
      `se(coef)` = x$std.error,
      z = x$zvalue,
      `p` = x$pvalue
    )
    printCoefmat(coef_table, digits = digits, 
                 P.values = TRUE, has.Pvalue = TRUE,
                 signif.stars = getOption("show.signif.stars"), ...)
  }
  
  if (!x$converged) {
    cat("\nWarning: Algorithm did not converge\n")
  }
  
  invisible(x)
}


#' @export
print.summary.crrstep <- function(x, digits = max(3L, getOption("digits") - 3L),
                                  signif.stars = getOption("show.signif.stars"), ...) {
  
  cat("\nStepwise Fine-Gray Competing Risks Regression\n")
  cat(paste(rep("=", 50), collapse = ""), "\n")
  
  cat("\nCall:\n")
  print(x$call)
  
  cat("\nSelection method:", x$direction, "\n")
  cat("Selection criterion:", x$criterion, "\n")
  
  cat("\nSelected model:\n")
  if (x$is_null_model) {
    cat("  ~ 1 (Null model)\n")
  } else {
    cat(" ", deparse(x$formula, width.cutoff = 500), "\n")
  }
  
  cat("\n  n =", x$n, ", number of events =", x$nevent, "\n")
  
  if (!x$is_null_model) {
    cat("\nCoefficients:\n")
    printCoefmat(x$coefficients, digits = digits,
                 P.values = TRUE, has.Pvalue = TRUE,
                 signif.stars = signif.stars, ...)
    
    cat("\n", format(x$conf.level * 100), "% Confidence intervals for hazard ratios:\n", sep = "")
    print(x$conf.int, digits = digits)
  }
  
  cat("\n")
  cat(paste(rep("-", 50), collapse = ""), "\n")
  cat(x$criterion, ":", format(round(x$criterion_value, 2), nsmall = 2), "\n")
  cat("Log-likelihood:", format(round(x$loglik, 2), nsmall = 2), "\n")
  cat("Null log-likelihood:", format(round(x$loglik_null, 2), nsmall = 2), "\n")
  
  if (!x$is_null_model) {
    cat("\nLikelihood ratio test: ", 
        format(round(x$lr_test$statistic, 2), nsmall = 2), " on ",
        x$lr_test$df, " df, p = ", 
        format.pval(x$lr_test$pvalue, digits = digits), "\n", sep = "")
  }
  
  if (!x$converged) {
    cat("\nWarning: Algorithm did not converge\n")
  }
  
  invisible(x)
}

#' @export
summary.crrstep <- function(object, conf.level = 0.95, ...) {
  
  # Recalculate confidence intervals if different level requested
  if (!object$is_null_model && conf.level != object$conf.level) {
    z_crit <- qnorm((1 + conf.level) / 2)
    conf.int <- cbind(
      lower = object$coefficients - z_crit * object$std.error,
      upper = object$coefficients + z_crit * object$std.error
    )
    rownames(conf.int) <- names(object$coefficients)
  } else {
    conf.int <- object$conf.int
  }
  
  # Build coefficient table
  if (!object$is_null_model) {
    coef_table <- cbind(
      coef = object$coefficients,
      `exp(coef)` = exp(object$coefficients),
      `se(coef)` = object$std.error,
      z = object$zvalue,
      `Pr(>|z|)` = object$pvalue
    )
    
    # Confidence intervals for hazard ratios
    conf.int.hr <- cbind(
      `exp(coef)` = exp(object$coefficients),
      `exp(-coef)` = exp(-object$coefficients),
      `lower` = exp(conf.int[, "lower"]),
      `upper` = exp(conf.int[, "upper"])
    )
    colnames(conf.int.hr)[3:4] <- paste0(c("lower ", "upper "), 
                                         format(conf.level * 100), "%")
  } else {
    coef_table <- NULL
    conf.int.hr <- NULL
  }
  
  # Likelihood ratio test vs null model
  if (!object$is_null_model) {
    lr_stat <- 2 * (object$loglik - object$loglik_null)
    lr_df <- length(object$coefficients)
    lr_pvalue <- 1 - pchisq(lr_stat, df = lr_df)
  } else {
    lr_stat <- 0
    lr_df <- 0
    lr_pvalue <- 1
  }
  
  result <- list(
    call = object$call,
    formula = object$formula,
    direction = object$direction,
    criterion = object$criterion,
    criterion_value = object$criterion_value,
    coefficients = coef_table,
    conf.int = conf.int.hr,
    conf.level = conf.level,
    loglik = object$loglik,
    loglik_null = object$loglik_null,
    lr_test = list(statistic = lr_stat, df = lr_df, pvalue = lr_pvalue),
    n = object$n,
    nevent = object$nevent,
    converged = object$converged,
    is_null_model = object$is_null_model
  )
  
  class(result) <- "summary.crrstep"
  return(result)
}



#' @export
coef.crrstep <- function(object, ...) {
  object$coefficients
}


#' @export
vcov.crrstep <- function(object, ...) {
  if (object$is_null_model) {
    return(NULL)
  }
  object$crr_fit$var
}


#' @export
confint.crrstep <- function(object, parm, level = 0.95, ...) {
  if (object$is_null_model) {
    return(NULL)
  }
  
  z_crit <- qnorm((1 + level) / 2)
  cf <- object$coefficients
  se <- object$std.error
  
  ci <- cbind(
    cf - z_crit * se,
    cf + z_crit * se
  )
  
  pct <- format(100 * c((1 - level) / 2, (1 + level) / 2), trim = TRUE, digits = 3)
  colnames(ci) <- paste(pct, "%")
  
  if (!missing(parm)) {
    ci <- ci[parm, , drop = FALSE]
  }
  
  ci
}


#' @export
logLik.crrstep <- function(object, ...) {
  val <- object$loglik
  attr(val, "df") <- if (object$is_null_model) 0 else length(object$coefficients)
  attr(val, "nobs") <- object$n
  class(val) <- "logLik"
  val
}


#' @export
AIC.crrstep <- function(object, ..., k = 2) {
  ll <- object$loglik
  
  df <- if (object$is_null_model) 0 else length(object$coefficients)
  -2 * ll + k * df
}


#' @export
BIC.crrstep <- function(object, ...) {
  AIC(object, k = log(object$n))
}