R/plsi.log.auto.constraint.R

In EPLSIM: Partial Linear Single Index Models for Environmental Mixture Analysis

#' Partial linear single index log-linear (count) regression with automatic smoothness selection and built-in identifiability constraint
#'
#' Count-outcome analog of \code{plsi.lr.auto.constraint()}: the same model as
#' \code{plsi.log.auto()} (penalized regression spline link via
#' \code{mgcv::gam} with a log link, negative binomial or Poisson family, REML
#' smoothing-parameter selection), but with the single-index direction
#' searched over a reparameterization that satisfies the identifiability
#' constraint by construction instead of normalizing after the fact.
#'
#' The single-index coefficient \eqn{\beta} is constrained to \eqn{\|\beta\| = 1}
#' and \eqn{\beta_1 > 0} for identifiability. The search is reparameterized in
#' terms of \eqn{\sigma = (\sigma_1, \dots, \sigma_{p-1})} via
#' \eqn{\beta(\sigma) = \{(1 - \|\sigma\|^2)^{1/2}, \sigma_1, \dots,
#' \sigma_{p-1}\}^T}. Both constraints hold automatically for any
#' \eqn{\sigma} with \eqn{\|\sigma\| < 1}, so \code{optim()} runs unconstrained
#' over \eqn{\sigma}. Standard errors for \eqn{\beta} are obtained from the
#' Hessian in \eqn{\sigma}-space via the delta method.
#'
#' The anchor exposure, the one assigned \eqn{\beta_1}, the coordinate
#' constrained \eqn{> 0}, is chosen as the candidate exposure with the
#' largest \emph{absolute} correlation with the count outcome, so an exposure
#' strongly but negatively associated with the outcome can still be selected
#' as the anchor. Which exposure was selected is reported in
#' \code{attr(si.coefficient, "note")} and in \code{anchor.note} of the
#' returned list.
#'
#' @param data A data set including all needed variables
#' @param Y.name Variable name for the count outcome. Must be non-negative
#'   integers (0, 1, 2, ...).
#' @param X.name Variable name vector for exposures (at least 2).
#' @param Z.name Variable name vector for confounders
#' @param family Count distribution to use. \code{"nb"} (default) fits a
#'   negative binomial GAM (\code{mgcv::nb()}), which estimates its own
#'   dispersion parameter and is the safer default for real, usually
#'   overdispersed, count data. \code{"poisson"} fits a Poisson GAM, valid
#'   only if the mean-equals-variance assumption holds.
#' @param k Upper bound on basis dimension for the smooth link function.
#'   Default 10.
#' @param bs Smooth basis type passed to \code{mgcv::s()}; default \code{"cr"}
#'   (cubic regression spline).
#' @param initial.random.num Number of random initials for the single-index
#'   direction search
#' @param seed A single integer used to seed the RNG so results (random initials
#'   and downstream optim() convergence) are reproducible. Set to NULL to skip
#'   seeding and use whatever RNG state is already active in the calling
#'   environment. Default is 2026.
#' @param maxit Maximum number of Nelder-Mead iterations per starting value
#'   for the \eqn{\sigma}-space search. Default 100, matching
#'   \code{plsi.lr.auto.constraint()}. With many exposures this is often too
#'   few for the simplex to converge, which can leave the numerical
#'   information matrix indefinite and the single-index SEs inflated; check
#'   \code{optim.convergence} (and the warning issued when the selected start
#'   did not converge) and increase \code{maxit} if needed. Ignored when there
#'   are exactly two exposures (Brent search).
#' @return A list of model estimation and prediction results, structured
#'   analogously to `plsi.log.auto()`'s output, plus:
#'   - `anchor.note`: Which exposure was used as the anchor (\eqn{\beta_1 > 0})
#'     and why; also attached as `attr(si.coefficient, "note")`.
#'   - `initial.sigma.table`: The starting values actually used by the
#'     \eqn{\sigma}-space search.
#'   - `optim.convergence`: For every starting value, the `optim()`
#'     convergence code (0 = converged, 1 = `maxit` reached), the number of
#'     log-likelihood evaluations, and the final log-likelihood.
#'
#'   The remaining elements are as in `plsi.log.auto()`:
#'   - `si.coefficient`: Single-index direction estimates and inference.
#'   - `confounder.coefficient`: Confounder log-rate coefficients, plus rate
#'     ratios and their 95% CIs.
#'   - `si.fun`: The estimated single-index link function, on both the log
#'     scale (`fit`/`se`/`lwr`/`upr`) and the count scale
#'     (`count.fit`/`count.lwr`/`count.upr`).
#'   - `si.fun.model`: A confounder-free count `gam` on the single index
#'     alone (confounder contribution held as a fixed offset); predict() from
#'     this needs only `single_index_estimated` (and `confounder_offset = 0`).
#'   - `full.model`: The full fitted `gam` (smooth + confounders).
#'
#' @details
#' The confounder-free single-index model is built exactly as in
#' \code{plsi.log.auto()}: the confounder contribution from the full model is
#' fixed as a known \code{offset()} and a second count \code{gam} (same
#' \code{family}) is refit directly against the observed \code{y}.
#'
#' Standard errors are computed from a numerical Hessian in a local chart of
#' the unit sphere anchored on the exposure with the largest \eqn{|\beta_j|}
#' (not necessarily the \eqn{\beta_1 > 0} anchor), then mapped to
#' \eqn{\beta} by the delta method. This keeps the finite-difference steps
#' away from the \eqn{\|\sigma\| = 1} boundary when the estimated
#' \eqn{\beta_1} is small; the point estimates are unaffected.
#'
#' With exactly two exposures, \eqn{\sigma} is one-dimensional; the search then
#' uses \code{optim(method = "Brent")} on \eqn{(-1, 1)} instead of Nelder-Mead,
#' which \code{optim()} flags as unreliable in one dimension.
#'
#' The returned object works with all plotting functions using
#' \code{type = "log"}.
#'
#' @importFrom stats qnorm
#' @importFrom mgcv gam nb
#'
#' @examples
#' \dontrun{
#' data(nhanes.new)
#' data <- nhanes.new
#'
#' X.name <- c("X1_trans.b.carotene", "X2_retinol", "X3_g.tocopherol", "X4_a.tocopherol",
#'             "X5_PCB99", "X6_PCB156", "X7_PCB206",
#'             "X8_3.3.4.4.5.pncb", "X9_1.2.3.4.7.8.hxcdf", "X10_2.3.4.6.7.8.hxcdf")
#' Z.name <- c("AGE.c", "SEX.Female", "RACE.NH.Black",
#'            "RACE.MexicanAmerican", "RACE.OtherRace", "RACE.Hispanic" )
#'
#' # demo count outcome (illustrative only as nhanes.new has no native count
#' # variable), simulated from a true single-index combination of the exposures
#' set.seed(2026)
#' beta_true <- c(0.30, -0.20, 0.10, 0.40, -0.30, 0.20, -0.10, 0.25, -0.15, 0.35)
#' beta_true <- beta_true / sqrt(sum(beta_true^2))
#' x_std <- scale(data[, X.name])
#' single_index_true <- as.vector(x_std %*% beta_true)
#' log_rate <- 0.3 + 0.4 * sin(single_index_true) + 0.05 * data$AGE.c
#' data$n.events <- stats::rpois(nrow(data), lambda = exp(log_rate))
#'
#' model_log_con <- plsi.log.auto.constraint(data = data, Y.name = "n.events",
#'                    X.name = X.name, Z.name = Z.name, family = "nb",
#'                    k = 10, bs = "cr", initial.random.num = 1, seed = 2026)
#' model_log_con$si.coefficient
#' si.fun.plot(model_log_con$si.fun, type = "log")
#' }
#'
#' @keywords partial linear single index
#' @keywords log-linear regression
#' @keywords automatic smoothness
#' @keywords count outcome
#' @author Yuyan Wang
#' @export
plsi.log.auto.constraint <- function(data, Y.name, X.name, Z.name,
                                     family = c("nb", "poisson"),
                                     k = 10, bs = "cr", initial.random.num = 5,
                                     seed = 2026, maxit = 100) {

  family <- match.arg(family)

  stopifnot(initial.random.num >= 1, maxit >= 1)
  if (length(X.name) < 2) {
    stop("plsi.log.auto.constraint() requires at least 2 exposures in X.name (the sigma ",
         "reparameterization needs p - 1 >= 1 free coordinates).")
  }
  if (!is.null(seed)) set.seed(seed)
  data <- as.data.frame(data)
  missing_cols <- setdiff(c(Y.name, X.name, Z.name), colnames(data))
  if (length(missing_cols) > 0) {
    stop("Column(s) not found in data: ", paste(missing_cols, collapse = ", "))
  }
  if (anyNA(data[, c(Y.name, X.name, Z.name)])) {
    stop("Missing values found in Y/X/Z columns; please handle NAs before calling ",
         "plsi.log.auto.constraint().")
  }
  if (!requireNamespace("mgcv", quietly = TRUE)) {
    stop("Package 'mgcv' is required for plsi.log.auto.constraint(). Please install it.")
  }

  # Validate the count outcome ---------------------------------------------
  y_raw <- data[, Y.name]
  if (!is.numeric(y_raw) || any(y_raw < 0) || any(y_raw != round(y_raw))) {
    stop("Y.name must be non-negative integer counts for plsi.log.auto.constraint(). ",
         "Found negative, non-integer, or non-numeric values.")
  }
  y <- as.numeric(y_raw)
  if (stats::var(y) == 0) {
    stop("Y.name has no variation (all values equal); cannot fit a count model.")
  }

  fam_obj <- if (family == "nb") mgcv::nb() else stats::poisson()

  # Anchor exposure: largest |correlation| with the count outcome (see
  # plsi.lr.auto.constraint() for why |r| rather than signed r).
  cor_linear <- stats::cor(data.frame(y = y, data[, X.name]))
  cor_linear <- as.data.frame(cor_linear[-1, 1])
  colnames(cor_linear)[1] <- "cor_with_Y"
  cor_linear$exposure <- X.name
  cor_linear <- cor_linear[order(abs(cor_linear$cor_with_Y), decreasing = TRUE), ]
  X.name <- cor_linear$exposure

  x <- as.matrix(data[, X.name]); z_df <- data[, Z.name, drop = FALSE]
  n <- nrow(data); x_length <- length(X.name); z_length <- length(Z.name)

  # Initial linear direction: MASS::glm.nb() for family = "nb", Poisson glm
  # otherwise (as in plsi.log.auto()). Only used for starting values and the
  # neg_penalty scale, not for inference. glm.nb() can fail to converge on
  # near-Poisson data, so fall back to the Poisson glm in that case.
  m0 <- NULL
  if (family == "nb") {
    m0 <- tryCatch(suppressWarnings(MASS::glm.nb(y ~ x + as.matrix(z_df))),
                   error = function(e) NULL)
  }
  if (is.null(m0)) {
    m0 <- stats::glm(y ~ x + as.matrix(z_df), family = stats::poisson())
  }

  # beta(sigma) reparameterization for the ||beta|| = 1, beta_1 > 0 constraint
  # (identical to plsi.lr.auto.constraint()).
  beta_from_sigma <- function(sigma) {
    ss <- sum(sigma^2)
    if (ss >= 1 - 1e-8) return(NULL)
    c(sqrt(1 - ss), sigma)
  }

  initial_table <- as.data.frame(matrix(NA, nrow = (1 + initial.random.num), ncol = x_length))
  colnames(initial_table) <- X.name
  rownames(initial_table) <- c("linear", paste0("random_", 1:initial.random.num))
  lin_coef <- m0$coefficients[2:(1 + x_length)]
  # glm() can return NA coefficients (aliasing);
  # fall back to a random start for the "linear" row if so.
  if (any(!is.finite(lin_coef)) || all(lin_coef == 0)) {
    warning("Linear count-model starting direction is not finite (possible collinearity); ",
            "replacing it with a random start.", call. = FALSE)
    lin_coef <- stats::runif(x_length, -1, 1)
  }
  initial_table[1, ] <- lin_coef
  initial_table[2:(1 + initial.random.num), ] <- stats::runif(initial.random.num * x_length, -1, 1)
  for (i in 1:nrow(initial_table)) {
    initial_table[i, ] <- initial_table[i, ] * sign(initial_table[i, 1]) / sqrt(sum(initial_table[i, ]^2))
  }
  initial_sigma_table <- initial_table[, -1, drop = FALSE]
  colnames(initial_sigma_table) <- paste0("sigma_", seq_len(x_length - 1))

  # REML-penalized smooth log link, refit for each candidate direction beta.
  fit_gam <- function(beta_0) {
    u0 <- as.vector(x %*% as.vector(beta_0))
    dat_gam <- cbind(data.frame(y = y, u0 = u0), z_df)
    form <- stats::as.formula(paste0("y ~ s(u0, bs='", bs, "', k=", k, ") + ",
                                     paste(Z.name, collapse = " + ")))
    mgcv::gam(form, data = dat_gam, family = fam_obj, method = "REML")
  }

  # Finite penalty for infeasible sigma / failed fits (see
  # plsi.lr.auto.constraint() for why this must not be -Inf).
  ref_loglik <- as.numeric(stats::logLik(m0))
  neg_penalty <- ref_loglik - 1e6 * (abs(ref_loglik) + 1)

  fn <- function(sigma) {
    utils::flush.console()
    beta_0 <- beta_from_sigma(sigma)
    if (is.null(beta_0)) return(neg_penalty)
    ll <- tryCatch(as.numeric(stats::logLik(fit_gam(beta_0))),
                   error = function(e) {
                     warning("fit_gam() failed for a candidate beta and was scored as a large ",
                             "penalty: ", conditionMessage(e), call. = FALSE)
                     neg_penalty
                   })
    if (!is.finite(ll)) neg_penalty else ll
  }

  # One free coordinate (p = 2): Brent on (-1, 1); otherwise Nelder-Mead.
  run_optim <- function(sigma_0) {
    if (length(sigma_0) == 1) {
      stats::optim(par = sigma_0, fn = fn, method = "Brent",
                   lower = -1 + 1e-6, upper = 1 - 1e-6,
                   control = list(fnscale = -1))
    } else {
      stats::optim(par = sigma_0, fn = fn, method = "Nelder-Mead",
                   control = list(fnscale = -1, maxit = maxit))
    }
  }

  fit_results <- vector("list", nrow(initial_sigma_table))
  crit_table <- data.frame(logLik = NA_real_, AIC = NA_real_, BIC = NA_real_)[rep(1, nrow(initial_sigma_table)), ]
  rownames(crit_table) <- rownames(initial_sigma_table)

  for (i in 1:nrow(initial_sigma_table)) {
    sigma_0 <- as.numeric(initial_sigma_table[i, ])
    m_temp <- run_optim(sigma_0)
    fit_results[[i]] <- m_temp
    g <- fit_gam(beta_from_sigma(m_temp$par))
    crit_table$logLik[i] <- as.numeric(stats::logLik(g))
    crit_table$AIC[i] <- stats::AIC(g)
    crit_table$BIC[i] <- stats::BIC(g)
  }

  optim_convergence <- data.frame(
    convergence = vapply(fit_results, function(r) as.integer(r$convergence), integer(1)),
    fn.evals = vapply(fit_results, function(r) if (is.na(r$counts[1])) NA_real_ else as.numeric(r$counts[1]), numeric(1)),
    logLik = vapply(fit_results, function(r) as.numeric(r$value), numeric(1)),
    row.names = rownames(initial_sigma_table))

  order_ob <- order(crit_table$BIC)[1]
  m_selected <- fit_results[[order_ob]]
  if (m_selected$convergence != 0) {
    warning("The selected starting value ('", rownames(crit_table)[order_ob], "') did not converge ",
            "within maxit = ", maxit, " Nelder-Mead iterations (optim convergence code ",
            m_selected$convergence, "). The single-index estimates may not be at the optimum and ",
            "their standard errors may be unreliable; consider increasing maxit. See ",
            "optim.convergence in the output.", call. = FALSE)
  }
  sigma_est <- m_selected$par
  beta_est <- beta_from_sigma(sigma_est)

  # Standard errors for beta ---------------------------------------------------
  # The search above is done in the sigma chart anchored on X.name[1]. If the
  # estimated beta_1 is small (the anchor is chosen by *marginal* |r|, which
  # need not match its weight in the fitted index), sigma_est sits close to
  # the ||sigma|| = 1 boundary and optimHess()'s finite-difference steps can
  # cross it, hitting neg_penalty and producing a finite but meaningless
  # Hessian (eigenvalues on the order of the penalty / step^2), which the
  # eigenvalue floor would then turn into near-zero SEs.
  #
  # To avoid this, the Hessian is computed in a *local* chart around beta_est
  # anchored on the exposure with the largest |beta_j| (with its sign held
  # fixed). That coordinate satisfies |beta_j| >= 1/sqrt(p), so the optimum is
  # always well inside the chart. Both charts parameterize the same unit
  # sphere and the delta-method covariance of beta depends only on the tangent
  # space at beta_est, so this changes nothing when the Hessian in the
  # original chart is already well behaved (and the chart is identical
  # whenever the largest |beta_j| is the anchor). beta_est itself, and the
  # beta_1 > 0 convention on X.name[1], are unchanged.
  j_star <- which.max(abs(beta_est))
  s_star <- sign(beta_est[j_star])
  tau_est <- beta_est[-j_star]
  beta_from_tau <- function(tau) {
    ss <- sum(tau^2)
    if (ss >= 1 - 1e-8) return(NULL)
    b <- numeric(x_length)
    b[j_star] <- s_star * sqrt(1 - ss)
    b[-j_star] <- tau
    b
  }
  hess_penalty_hits <- 0
  fn_tau <- function(tau) {
    beta_0 <- beta_from_tau(tau)
    if (is.null(beta_0)) {
      hess_penalty_hits <<- hess_penalty_hits + 1
      return(neg_penalty)
    }
    ll <- tryCatch(as.numeric(stats::logLik(fit_gam(beta_0))), error = function(e) NA_real_)
    if (!is.finite(ll)) {
      hess_penalty_hits <<- hess_penalty_hits + 1
      return(neg_penalty)
    }
    ll
  }

  hess <- tryCatch(
    stats::optimHess(par = tau_est, fn = fn_tau, control = list(fnscale = -1)),
    error = function(e) {
      warning("optimHess() could not compute a finite-difference Hessian at the selected optimum ",
              "(", conditionMessage(e), "). Falling back to a conservative (deliberately ",
              "large-variance) covariance for beta; consider increasing initial.random.num.",
              call. = FALSE)
      NULL
    }
  )
  if (!is.null(hess) && hess_penalty_hits > 0) {
    warning("The finite-difference Hessian hit ", hess_penalty_hits, " infeasible or failed ",
            "gam() evaluation(s), so it is not reliable. Falling back to a conservative ",
            "(deliberately large-variance) covariance for beta.", call. = FALSE)
    hess <- NULL
  }

  if (is.null(hess)) {
    cov_tau <- diag(1e8, length(tau_est))  # deliberately large fallback variance
  } else {
    info_matrix <- -hess
    eig_decomp <- eigen(info_matrix, symmetric = TRUE)
    eig_vals <- eig_decomp$values
    if (any(eig_vals <= 0)) {
      warning("Observed information matrix is not positive definite at the selected optimum ",
              "(min eigenvalue = ", signif(min(eig_vals), 3), "). Standard errors are computed from ",
              "an eigenvalue-floored version of the information matrix and should be treated as ",
              "conservative approximations; consider increasing initial.random.num for a more ",
              "reliable optimum.")
      # See plsi.lr.auto.constraint(): floor tiny/negative eigenvalues (caused
      # by the inner REML refit's small discontinuities) relative to the
      # largest one, so the affected SEs are conservative rather than NaN.
      floor_val <- max(abs(eig_vals)) * 1e-6
      eig_vals <- pmax(eig_vals, floor_val)
      info_matrix <- eig_decomp$vectors %*% diag(eig_vals, nrow = length(eig_vals)) %*% t(eig_decomp$vectors)
    }
    cov_tau <- suppressWarnings(MASS::ginv(info_matrix))
  }

  # Delta method from the local chart to beta: J = d(beta)/d(tau), with
  #   d(beta_jstar)/d(tau) = -tau / beta_jstar   (row j_star)
  #   d(beta_-jstar)/d(tau) = I                  (all other rows)
  # beta_jstar is bounded away from 0 (|beta_jstar| >= 1/sqrt(p)), so this
  # Jacobian is always well conditioned.
  jac <- matrix(0, nrow = x_length, ncol = x_length - 1)
  jac[-j_star, ] <- diag(x_length - 1)
  jac[j_star, ] <- -tau_est / beta_est[j_star]
  cov_beta <- jac %*% cov_tau %*% t(jac)
  beta_sigma <- suppressWarnings(sqrt(diag(cov_beta)))

  final_model <- fit_gam(beta_est)
  s <- summary(final_model)

  # Wald tests use the normal reference distribution (z), matching
  # plsi.log.auto() and summary.gam() under a Poisson/NB family.
  si.coefficient <- data.frame(Estimate = as.numeric(beta_est),
                               `Std.Error` = as.numeric(beta_sigma),
                               row.names = X.name, check.names = FALSE)
  si.coefficient$`z value` <- si.coefficient$Estimate / si.coefficient$`Std.Error`
  si.coefficient$`Pr(>|z|)` <- ifelse(2 * stats::pnorm(-abs(si.coefficient$`z value`)) < 0.0001, "<.0001",
                                      format(round(2 * stats::pnorm(-abs(si.coefficient$`z value`)), 4), nsmall = 4))
  si.coefficient$Lower.95CI <- si.coefficient$Estimate + stats::qnorm(0.025) * si.coefficient$`Std.Error`
  si.coefficient$Upper.95CI <- si.coefficient$Estimate + stats::qnorm(0.975) * si.coefficient$`Std.Error`
  si.coefficient$`Contribution proportion` <- format(round((si.coefficient$Estimate)^2, 3), nsmall = 3)
  si.coefficient <- si.coefficient[order(si.coefficient$Estimate, decreasing = TRUE), ]

  anchor.exposure <- X.name[1]
  anchor.cor <- cor_linear$cor_with_Y[cor_linear$exposure == anchor.exposure]
  anchor.note <- paste0(
    "First/anchor exposure (single-index coefficient constrained > 0 for identifiability, ",
    "||beta|| = 1): '", anchor.exposure, "', selected as the candidate exposure with the ",
    "largest absolute correlation with ", Y.name, " (r = ", signif(anchor.cor, 3), ")."
  )
  attr(si.coefficient, "note") <- anchor.note

  confounder.coefficient <- as.data.frame(s$p.table[Z.name, , drop = FALSE])
  confounder.coefficient$Lower.95CI <- confounder.coefficient$Estimate + stats::qnorm(0.025) * confounder.coefficient$`Std. Error`
  confounder.coefficient$Upper.95CI <- confounder.coefficient$Estimate + stats::qnorm(0.975) * confounder.coefficient$`Std. Error`
  confounder.coefficient$`Rate Ratio` <- exp(confounder.coefficient$Estimate)
  confounder.coefficient$RR.Lower.95CI <- exp(confounder.coefficient$Lower.95CI)
  confounder.coefficient$RR.Upper.95CI <- exp(confounder.coefficient$Upper.95CI)
  confounder.coefficient$`Pr(>|z|)` <- ifelse(confounder.coefficient$`Pr(>|z|)` < 0.0001, "<.0001",
                                              format(round(confounder.coefficient$`Pr(>|z|)`, 4), nsmall = 4))

  single_index_estimated <- as.vector(x %*% beta_est)

  # Confounder-free single-index model for prediction (offset approach; see
  # plsi.log.auto() for why this is fit against observed y).
  confounder_coefs <- stats::coef(final_model)[Z.name]
  confounder_offset <- as.vector(as.matrix(z_df) %*% confounder_coefs)

  dat_si_fun <- data.frame(single_index_estimated = single_index_estimated,
                           y = y,
                           confounder_offset = confounder_offset)
  form_si <- stats::as.formula(paste0("y ~ offset(confounder_offset) + s(single_index_estimated, bs='",
                                      bs, "', k=", k, ")"))
  m2 <- mgcv::gam(form_si, data = dat_si_fun, family = fam_obj, method = "REML")

  ord <- order(single_index_estimated)
  newdat <- data.frame(single_index_estimated = single_index_estimated[ord],
                       confounder_offset = 0)
  pred2 <- stats::predict(m2, newdata = newdat, type = "link", se.fit = TRUE)
  si.fun <- data.frame(single_index_estimated = single_index_estimated[ord],
                       fit = pred2$fit,
                       se  = pred2$se.fit)
  si.fun$lwr <- si.fun$fit - stats::qnorm(0.975) * si.fun$se
  si.fun$upr <- si.fun$fit + stats::qnorm(0.975) * si.fun$se
  si.fun$count.fit <- exp(si.fun$fit)
  si.fun$count.lwr <- exp(si.fun$lwr)
  si.fun$count.upr <- exp(si.fun$upr)

  return(list(original.data = list(y = y, x = x, z = z_df),
              original.par = list(family = family, k = k, bs = bs,
                                  initial.random.num = initial.random.num, seed = seed,
                                  maxit = maxit),
              si.coefficient = si.coefficient,          # carries attr(., "note") on anchor-exposure selection
              anchor.note = anchor.note,
              model.statistics = crit_table[order_ob, ],
              initial.table = initial_table,
              initial.sigma.table = initial_sigma_table,
              optim.convergence = optim_convergence,
              confounder.coefficient = confounder.coefficient,
              si.fun.edf = s$edf[1],
              si.fun = si.fun,
              si.fun.model = m2,              # confounder-free, log scale
              full.model = final_model))
}