R/nomogram.R

In QHScrnomo: Construct Nomograms for Competing Risks Regression Models

##' Construct a Nomogram for a Competing Risks Regression Model
##'
##' Draws a partial nomogram from a \code{\link[QHScrnomo]{crr.fit}} object that can be used to manually obtain predicted values from from a competing risks regression model.
##'
##' @param fit A model fit by \code{\link[QHScrnomo]{crr.fit}}
##' @param failtime A vector of time points to display failure probability axes for.
##' @param ci Should the failure probability be displayed? Defaults to \code{TRUE}. If \code{FALSE}, the event-free probability is displayed.
##' @param ... Settings of variables to use in constructing axes.
##' If \code{datadist} was in effect, the default is to use \code{pretty(total range, nint)}
##' for continuous variables, and the class levels for discrete ones.
##' For \code{legend.nomabbrev}, \code{...} specifies optional parameters to pass to the legend.
##' Common ones are \code{bty = "n"} to suppress drawing the box. You may want to
##' specify a non-proportionally spaced font (e.g., courier) number if
##' abbreviations are more than one letter long. This will make the abbreviation
##' definitions line up (e.g., specify font = 2, the default for courier).
##' Ignored for print.
##' @param adj.to If \code{datadist} was not defined for all predictors, it is required to define adjustment settings for the undefined ones (e.g. \code{adj.to=list(age=50, sex="female")}).
##' @param lp Set to \code{FALSE} to suppress creation of an axis for scoring \eqn{X\beta}{X beta}. Defaults to \code{TRUE}.
##' @param lp.at If \code{lp=TRUE}, specifies a vector of settings of \eqn{X\beta}{X beta}. Default is to use \code{pretty(range of linear predictors, nint)}.
##' @param lplabel label for linear predictor axis. Default is \code{"Linear Predictor"}.
##' @param fun.at Function values to label on axis. Default \code{fun}
##'   evaluated at \code{lp.at}. If more than one \code{fun} was specified,
##'   using a vector for \code{fun.at} will cause all functions to be evaluated
##'   at the same argument values. To use different values, specify a list of
##'   vectors for \code{fun.at}, with elements corresponding to the different
##'   functions (lists of vectors also applies to \code{fun.lp.at} and
##'   \code{fun.side}).
##' @param fun.lp.at If you want to evaluate one of the functions at a
##'   different set of linear predictor values than may have been used in
##'   constructing the linear predictor axis, specify a vector or list of
##'   vectors of linear predictor values at which to evaluate the function.
##'   This is especially useful for discrete functions. The presence of this
##'   attribute also does away with the need for \code{nomogram} to compute
##'   numerical approximations of the inverse of the function.  It also allows
##'   the user-supplied function to return \code{factor} objects, which is
##'   useful when e.g. a single tick mark position actually represents a range.
##'   If the \code{fun.lp.at} parameter is present, the \code{fun.at} vector
##'   for that function is ignored.
##' @param funlabel Label for \code{fun} axis. If more than one function was
##'   given but \code{funlabel} is of length one, it will be duplicated as needed. If
##'   \code{fun} is a list of functions for which you specified names (see the
##'   final example below), these names will be used as labels.
##' @param fun.side A vector or list of vectors of \code{side} parameters for
##'   the \code{axis} function for labeling function values. Values may be 1 to
##'   position a tick mark label below the axis (the default), or 3 for above
##'   the axis. If for example an axis has 5 tick mark labels and the second
##'   and third will run into each other, specify \code{fun.side=c(1,1,3,1,1)}
##'   (assuming only one function is specified as \code{fun}).
##' @param interact When a continuous variable interacts with a discrete one,
##'   axes are constructed so that the continuous variable moves within the
##'   axis, and separate axes represent levels of interacting factors. For
##'   interactions between two continuous variables, all but the axis variable
##'   must have discrete levels defined in \code{interact}. For discrete
##'   interacting factors, you may specify levels to use in constructing the
##'   multiple axes. For continuous interacting factors, you must do this.
##'   Examples: \code{interact=list(age=seq(10,70,by=10),
##'   treat=c("A","B","D"))}.
##' @param intercept For models such as the ordinal logistic model with
##'   multiple intercepts, specifies which one to use in evaluating the linear
##'   predictor.
##' @param conf.int Confidence levels to display for each scoring. Default is
##'   \code{FALSE} to display no confidence limits. Setting \code{conf.int} to
##'   \code{TRUE} is the same as setting it to \code{c(0.7, 0.9)}, with the
##'   line segment between the 0.7 and 0.9 levels shaded using gray scale.
##' @param col.conf Colors corresponding to \code{conf.int}. Use fractions for
##'   gray scale (for UNIX S-PLUS).
##' @param conf.space A 2-element vector with the vertical range within which
##'   to draw confidence bars, in units of 1=spacing between main bars. Four
##'   heights are used within this range (8 for the linear predictor if more
##'   than 16 unique values were evaluated), cycling them among separate
##'   confidence intervals to reduce overlapping.
##' @param conf.lp Default is \code{"representative"} to group all linear
##'   predictors evaluated into deciles, and to show, for the linear predictor
##'   confidence intervals, only the mean linear predictor within the deciles
##'   along with the median standard error within the deciles. Set
##'   \code{conf.lp="none"} to suppress confidence limits for the linear
##'   predictors, and to \code{"all"} to show all confidence limits.
##' @param est.all To plot axes for only the subset of variables named in
##'   \code{\dots{}}, set \code{est.all=FALSE}. Note: This option only works
##'   when zero has a special meaning for the variables that are omitted from
##'   the graph.
##' @param abbrev Set to \code{TRUE} to use the \code{abbreviate} function to
##'   abbreviate levels of categorical factors, both for labeling tick marks
##'   and for axis titles. If you only want to abbreviate certain predictor
##'   variables, set \code{abbrev} to a vector of character strings containing
##'   their names.
##' @param minlength Applies if \code{abbrev=TRUE}. Is the minimum
##'   abbreviation length passed to the \code{abbreviate} function. If you set
##'   \code{minlength=1}, the letters of the alphabet are used to label tick
##'   marks for categorical predictors, and all letters are drawn no matter how
##'   close together they are. For labeling axes (interaction settings),
##'   \code{minlength=1} causes \code{minlength=4} to be used.
##' @param maxscale Default maximum point score is 100.
##' @param nint Number of intervals to label for axes representing continuous variables. See \code{\link{pretty}}.
##' @param label.every Specify \code{label.every=i} to label on every \code{i}th tick mark.
##' @param force.label Set to \code{TRUE} to force every tick mark intended to
##'   be labeled to have a label plotted (whether the labels run into each
##'   other or not)
##' @param xfrac Fraction of horizontal plot to set aside for axis titles
##' @param cex.axis Character size for tick mark labels
##' @param cex.var Character size for axis titles (variable names)
##' @param col.grid If \code{col.grid=1}, no gray scale is used, but an
##'   ordinary line is drawn. If \code{0<col.grid<1}, a \code{col} (gray
##'   scale) of \code{col.grid} is used to draw vertical reference lines for
##'   major axis divisions and \code{col.grid/2} for minor divisions. The
##'   default is \code{col.grid=FALSE}, i.e., reference lines are omitted.
##'   Specifying \code{col.grid=TRUE} is the same as specifying a gray scale
##'   level of \code{col.grid=.2} (5 for Windows S-PLUS).
##' @param vnames By default, variable labels are used to label axes. Set
##'   \code{vnames="names"} to instead use variable names.
##' @param varname.label In constructing axis titles for interactions, the
##'   default is to add \code{"(interacting.varname=level)"} on the right.
##'   Specify \code{varname.label=FALSE} to instead use \code{"(level)"}.
##' @param varname.label.sep If \code{varname.label=TRUE}, you can change the
##'   separator to something other than \code{=} by specifying this parameter.
##' @param ia.space When multiple axes are draw for levels of interacting
##'   factors, the default is to group combinations related to a main effect.
##'   This is done by spacing the axes for the second to last of these within a
##'   group only 0.7 (by default) of the way down as compared with normal space
##'   of 1 unit.
##' @param tck See \code{tck} under \code{par}.
##' @param lmgp Spacing between numeric axis labels and axis (see \code{par} for \code{mgp})
##' @param omit Vector of character strings containing names of variables for
##'   which to suppress drawing axes. Default is to show all variables.
##' @param naxes Maximum number of axes to allow on one plot. If the nomogram
##'   requires more than one "page", the "Points" axis will be repeated at the
##'   top of each page when necessary.
##' @param points.label A character string giving the axis label for the points scale
##' @param total.points.label A character string giving the axis label for the total points scale
##' @param total.sep.page Set to \code{TRUE} to force the total points and later axes to be placed on a separate page
##' @param total.fun A user-provided function that will be executed before the
##'   total points axis is drawn. Default is not to execute a function. This
##'   is useful e.g. when \code{total.sep.page=TRUE} and you wish to use
##'   \code{locator} to find the coordinates for positioning an abbreviation
##'   legend before it's too late and a new page is started (i.e.,
##'   \code{total.fun=function()print(locator(1))}).
##' @param verbose Set to \code{TRUE} to get printed output detailing how tick
##'   marks are chosen and labeled for function axes. This is useful in seeing
##'   how certain linear predictor values cannot be solved for using inverse
##'   linear interpolation on the (requested linear predictor values, function
##'   values at these lp values). When this happens you will see \code{NA}s in
##'   the \code{verbose} output, and the corresponding tick marks will not
##'   appear in the nomogram.
##' @param total.min Setting the minimal value in the total point axis on the nomogram.
##' @param total.max Setting the maximal value in the total point axis.
##' @param mikeomit The predictor variables specified by their names here will
##'   not be shown in the nomogram. The predicted outcome based on this
##'   reduced nomogram would be the same as if users were using the full
##'   version of the nomogram by entering the some values for the predictors
##'   remaining in the reduced nomogram but adjusted values for the hiden
##'   predictors so that 0 points will be achieved from these hiden predictor
##'   variables in the full nomogram.
##' @return A list of class \code{"nomogram"} that contains information used in
##'   plotting the axes. Please see \code{\link[rms]{nomogram}} for details.
##' @author Changhong Yu, Michael Kattan, Ph.D \cr Department of Quantitative
##'   Health Sciences\cr Cleveland Clinic\cr
##' @author Frank Harrell\cr
##' Department of Biostatistics\cr
##' Vanderbilt University\cr
##' \email{f.harrell@@vanderbilt.edu}
##'
##' @export
##' @seealso \code{\link[rms]{nomogram}} \code{\link[QHScrnomo]{nomogram.mk6}} \code{\link[QHScrnomo]{crr.fit}}
##' @references Banks J: Nomograms. Encylopedia of Statistical Sciences, Vol 6.
##'   Editors: S Kotz and NL Johnson.  New York: Wiley; 1985.
##'
##' Lubsen J, Pool J, van der Does, E: A practical device for the application
##'   of a diagnostic or prognostic function.  Meth. Inform. Med. 17:127--129;
##'   1978.
##'
##' Wikipedia: Nomogram, \url{https://en.wikipedia.org/wiki/Nomogram}.
##'
##' Michael W. Kattan, Glenn Heller and Murray F. Brennan (2003). A
##'   competing-risks nomogram \cr for sarcoma-specific death following local
##'   recurrence. Statistics in Medicine. \code{Stat Med}. 2003;22:3515-3525.
##' @examples
##' data(prostate.dat)
##' dd <- datadist(prostate.dat)
##' options(datadist = "dd")
##' prostate.f <- cph(Surv(TIME_EVENT,EVENT_DOD == 1) ~ TX  + rcs(PSA,3) +
##'            BX_GLSN_CAT +  CLIN_STG + rcs(AGE,3) +
##'            RACE_AA, data = prostate.dat,
##'            x = TRUE, y= TRUE, surv=TRUE,time.inc = 144)
##' prostate.crr <- crr.fit(prostate.f,cencode = 0,failcode = 1)
##' ## make a CRR nomogram
##' nomogram.crr(prostate.crr,failtime = 120,lp=FALSE,
##' funlabel = "Predicted 10-year cumulative incidence")
##'
##' @keywords models regression hplot
##'
nomogram.crr <-
  function(
      fit,
      failtime,
      ci = TRUE,
      ...,
      adj.to,
      lp = TRUE,
      lp.at,
      lplabel = "Linear Predictor",
      fun.at,
      fun.lp.at,
      funlabel = "Predicted Value",
      fun.side,
      interact = NULL,
      intercept = 1,
      conf.int = FALSE,
      col.conf = c(1, 12),
      conf.space = c(0.08, 0.2),
      conf.lp = c("representative", "all", "none"),
      est.all = TRUE,
      abbrev = FALSE,
      minlength = 4,
      maxscale = 100,
      nint = 10,
      label.every = 1,
      force.label = FALSE,
      xfrac = 0.35,
      cex.axis = 0.85,
      cex.var = 1,
      col.grid = FALSE,
      vnames = c("labels", "names"),
      varname.label = TRUE,
      varname.label.sep = "=",
      ia.space = 0.7,
      tck = -0.009,
      lmgp = 0.4,
      omit = NULL,
      naxes,
      points.label = "Points",
      total.points.label = "Total Points",
      total.sep.page = FALSE,
      total.fun,
      verbose = FALSE,
      total.min,
      total.max,
      mikeomit = NULL
    ) {

    # Check for fit
    if(missing(fit))
      stop("Please supply a model fit from crr.fit")

    # Check if the object was fit from crr.fit
    if(!inherits(fit, "cmprsk"))
      stop("The fit is not a 'cmprsk' object fit from crr.fit")

    # Check for a time point
    if(missing(failtime))
      stop("Please specify one or more time points for computing the event probability.")

    # Get the number of axes to draw
    nt <- length(failtime)

    # Extract the rms::cph object
    cph.f <- fit$cph.f

    # Assign the crr.fit coefficients to the rms::cph object
    nm <- names(cph.f$coefficients)
    cph.f$coefficients <- fit$coef
    names(cph.f$coefficients) <- nm

    # Set up list to hold transformation functions
    func <- mapply(list, rep(NA, nt))

    # Set a string for transformation function depending on failure vs. survival probability
    ci_string <- ""
    if(!ci)
      ci_string <- "1-"

    # Iterate the requested time points
    for(i in seq_len(nt)) {

      # Create and parse the string template of the transformation function for this time point
      expr <- parse(text = paste0("function(x)", ci_string, "nomo2.crr(x + cph.f$center, fit, time = ", failtime[i], ")"))

      # Set the function into the list
      func[[i]] <- eval(expr)

    }

    # Pass arguments to workhorse
    nomogram.mk6(
      cph.f,
      ...,
      adj.to = adj.to,
      lp = lp,
      lp.at = lp.at,
      lplabel = lplabel,
      fun = func,
      fun.at = fun.at,
      fun.lp.at = fun.lp.at,
      funlabel = funlabel,
      fun.side = fun.side,
      interact = interact,
      intercept = intercept,
      conf.int = conf.int,
      col.conf = col.conf,
      conf.space = conf.space,
      conf.lp = conf.lp,
      est.all = est.all,
      abbrev = abbrev,
      minlength = minlength,
      maxscale = maxscale,
      nint = nint,
      label.every = label.every,
      force.label = force.label,
      xfrac = xfrac,
      cex.axis = cex.axis,
      cex.var = cex.var,
      col.grid = col.grid,
      vnames = vnames,
      varname.label = varname.label,
      varname.label.sep = varname.label.sep,
      ia.space = ia.space,
      tck = tck,
      lmgp = lmgp,
      naxes = naxes,
      points.label = points.label,
      total.points.label = total.points.label,
      total.sep.page = total.sep.page,
      total.fun = total.fun,
      verbose = verbose,
      total.min = total.min,
      total.max = total.max,
      mikeomit = mikeomit
    )

  }

##' Construct a Nomogram for a Regression Model
##'
##' Draws a partial nomogram that can be used to manually obtain predicted values. This is a modified version of \code{\link[rms]{nomogram}}.
##'
##' @param fit A regression model that was created with \code{library(rms)}
##' in effect, and (usually) with \code{options(datadist = "object.name")} in
##' effect.
##' @param ... Settings of variables to use in constructing axes.
##' If \code{datadist} was in effect, the default is to use
##' \code{pretty(total range, nint)} for continuous variables, and the class
##' levels for discrete ones.
##' For \code{legend.nomabbrev}, \code{\dots} specifies optional parameters to
##' pass to \code{legend}.  Common ones are \code{bty = "n"} to suppress
##' drawing the box.  You may want to specify a non-proportionally spaced font
##' (e.g., courier) number if abbreviations are more than one letter long.
##' This will make the abbreviation definitions line up (e.g., specify
##' \code{font = 2}, the default for courier).  Ignored for \code{print}.
##' @param adj.to   If you didn't define \code{datadist} for all predictors,
##' you will have to define adjustment settings for the undefined ones, e.g.
##' \code{adj.to= list(age = 50, sex = "female")}.
##' @param lp  Set to \code{FALSE} to suppress creation of an axis for scoring
##' \eqn{X\beta}{X beta}.
##' @param lp.at If \code{lp=TRUE}, \code{lp.at} may specify a vector of
##' settings of \eqn{X\beta}{X beta}.
##' Default is to use \code{pretty(range of linear predictors, nint)}.
##' @param lplabel label for linear predictor axis.
##' Default is \code{"Linear Predictor"}.
##' @param fun on another axis.  If more than one transformation is plotted,
##' put them in a list, e.g. \code{list(function(x) x/2, function(x) 2*x)}.
##' Any function values equal to \code{NA} will be ignored.
##' @param fun.at function values to label on axis.  Default \code{fun}
##' evaluated at \code{lp.at}.   If more than one \code{fun} was specified,
##' using a vector for \code{fun.at} will cause all functions to be evaluated
##' at the same argument values.  To use different values, specify a list of
##' vectors for \code{fun.at}, with elements corresponding to the different
##' functions (lists of vectors also applies to \code{fun.lp.at} and
##' \code{fun.side}).
##' @param fun.lp.at If you want to
##' evaluate one of the functions at a different set of linear predictor
##' values than may have been used in constructing the linear predictor axis,
##' specify a vector or list of vectors
##' of linear predictor values at which to evaluate the function.  This is
##' especially useful for discrete functions.  The presence of this attribute
##' also does away with the need for \code{nomogram} to compute numerical
##' approximations of the inverse of the function.  It also allows the
##' user-supplied function to return \code{factor} objects, which is useful
##' when e.g. a single tick mark position actually represents a range.
##' If the \code{fun.lp.at} parameter is present, the \code{fun.at}
##' vector for that function is ignored.
##' @param funlabel label for \code{fun} axis.  If more than one function was
##' given but funlabel is of length one, it will be duplicated as needed.
##' If \code{fun} is a list of functions for which you specified names (see the
##' final example below), these names will be used as labels.
##' @param fun.side a vector or list of vectors of \code{side} parameters for
##' the \code{axis} function for labeling function values.
##' Values may be 1 to position a tick mark label below the axis (the default),
##' or 3 for above the axis.  If for example an axis has 5 tick mark labels
##' and the second and third will run into each other, specify
##' \code{fun.side=c(1,1,3,1,1)} (assuming only one function is specified as
##' \code{fun}).
##' @param interact When a continuous variable interacts with a discrete one,
##' axes are constructed so that the continuous variable moves within the axis,
##' and separate axes represent levels of interacting factors.  For
##' interactions between two continuous variables, all but the axis variable
##' must have discrete levels defined in \code{interact}.
##' For discrete interacting factors, you may specify levels to use in
##' constructing the multiple axes.  For continuous interacting factors,
##' you must do this.  Examples: \code{interact = list(age = seq(10,70,by=10),
##' treat = c("A","B","D"))}.
##' @param intercept for models such as the ordinal logistic model with
##' multiple intercepts, specifies which one to use in evaluating the linear
##' predictor.
##' @param conf.int confidence levels to display for each scoring.
##' Default is \code{FALSE} to display no confidence limits.
##' Setting \code{conf.int} to \code{TRUE} is the same as
##' setting it to \code{c(0.7, 0.9)},
##' with the line segment between the 0.7 and 0.9 levels shaded using
##' gray scale.
##' @param col.conf colors corresponding to \code{conf.int}.
##' Use fractions for gray scale(for UNIX S-PLUS).
##' @param conf.space a 2-element vector with the vertical range within which
##' to draw confidence bars, in units of 1=spacing between main bars.  Four
##' heights are used within this range (8 for the linear predictor if more than
##' 16 unique values were evaluated), cycling them among separate confidence
##' intervals to reduce overlapping
##' @param conf.lp default is \code{"representative"} to group all linear
##' predictors evaluated into deciles, and to show, for the linear predictor
##' confidence intervals, only the mean linear predictor within the deciles
##' along with the median standard error within the deciles.  Set
##' \code{conf.lp = "none"} to suppress confidence limits for the linear
##' predictors, and to \code{"all"} to show all confidence limits.
##' @param est.all To plot axes for only the subset of variables named in
##'  \code{\dots}, set \code{est.all = FALSE}.  Note: This option only works
##' when zero has a special meaning for the variables that are omitted from the
##' graph
##' @param abbrev Set to \code{TRUE} to use the \code{\link{abbreviate}}
##' function to abbreviate levels of categorical factors, both for labeling
##' tick marks and for axis titles. If you only want to abbreviate certain
##' predictor variables, set \code{abbrev} to a vector of character strings
##' containing their names.
##' @param minlength \code{\link{abbreviate}} function. If you set
##' \code{minlength = 1}, the letters of the alphabet are used to label tick
##' marks for categorical predictors, and all letters are drawn no matter how
##' close together they are.  For labeling axes (interaction settings),
##' \code{minlength = 1} causes \code{minlength = 4} to be used
##' @param maxscale default maximum point score is 100
##' @param nint number of intervals to label for axes representing continuous
##' variables.
##' See \code{\link{pretty}}.
##' @param label.every Specify \code{label.every = i} to label on every
##' \code{i}th tick mark
##' @param force.label set to \code{TRUE} to force every tick mark intended to
##' be labeled to have a label plotted (whether the labels run into each other
##' or not)
##' @param xfrac fraction of horizontal plot to set aside for axis titles
##' @param cex.axis character size for tick mark labels
##' @param cex.var character size for axis titles (variable names)
##' @param col.grid If left unspecified, no vertical reference lines are drawn.
##' Specify a vector of length one (to use the same color for both minor and
##' major reference lines) or two (corresponding to the color for the major and
##' minor divisions, respectively) containing colors, to cause vertical
##' reference lines to the top points scale to be drawn.  For R, a good choice
##' is \code{col.grid = gray(c(0.8, 0.95))}.
##' @param vnames By default, variable labels are used to label axes.  Set
##' \code{vnames = "names"} to instead use variable names.
##' @param varname.label In constructing axis titles for interactions, the
##' default is to add \code{(interacting.varname = level)} on the right.
##' Specify \code{varname.label = FALSE} to instead use \code{"(level)"}.
##' @param varname.label.sep If \code{varname.label = TRUE}, you can change the
##' separator to something other than
##' \code{=} by specifying this parameter.
##' @param ia.space When multiple axes are draw for levels of interacting
##' factors, the default is to group combinations related to a main effect.
##' This is done by spacing the axes for the second to last of these
##' within a group only
##' 0.7 (by default) of the way down as compared with normal space of 1 unit.
##' @param tck see \code{tck} under \code{\link{par}}
##' @param tcl length of tick marks in nomogram
##' @param lmgp spacing between numeric axis labels and axis
##' (see \code{\link{par}} for \code{mgp})
##' @param omit vector of character strings containing names of variables for
##' which to suppress drawing axes.  Default is to show all variables.
##' @param naxes  maximum number of axes to allow on one plot.  If the nomogram
##' requires more than one \dQuote{page}, the \dQuote{Points} axis will be
##' repeated at the top of each page when necessary.
##' @param points.label a character string giving the axis label for the points
##' scale
##' @param total.points.label  a character string giving the axis label for
##' the total points scale
##' @param total.sep.page set to \code{TRUE} to force the total points and
##' later axes to be placed on a separate page
##' @param total.fun  a user-provided function that will be executed before the
##' total points axis is drawn.  Default is not toe xecute a function.
##' This is useful e.g. when \code{total.sep.page = TRUE} and you wish to use
##' \code{locator} to find the coordinates for positioning an abbreviation
##' legend before it's too late and a new page is started (i.e.,
##' \code{total.fun = function() print(locator(1))}).
##' @param verbose set to \code{TRUE} to get printed output detailing how tick
##' marks are chosen and labeled for function axes.  This is useful in seeing
##' how certain linear predictor values cannot be solved for using inverse
##' linear interpolation on the (requested linear predictor values, function
##' values at these lp values).  When this happens you will see \code{NA}s in
##' the verbose output, and the corresponding tick marks will not appear in
##' the nomogram.
##' @param cap.labels logical: should the factor labels have their first
##' letter capitalized?
##' @param total.min the minimum point for the total point axis
##' @param total.max the maxmum point for the total point axis
##' @param survtime specified survival time for the predicted survival
##' probability
##' @param mikeomit a modified version of \code{omit}
##'
##' @return A \code{\link[rms]{nomogram}}
##' @export
##' @note Used internally for \code{\link[QHScrnomo]{nomogram.crr}}. See \code{\link[rms]{nomogram}}.
##' @references
##' Banks J: Nomograms. Encylopedia of Statistical Sciences, Vol 6.
##' Editors: S Kotz and NL Johnson.  New York: Wiley; 1985.
##'
##' Lubsen J, Pool J, van der Does, E: A practical device for the application
##' of a diagnostic or prognostic function.  Meth. Inform. Med. 17:127--129;
##' 1978.
##'
##' Wikipedia: Nomogram, \url{https://en.wikipedia.org/wiki/Nomogram}.
##'
##' @seealso \code{\link[QHScrnomo]{nomogram.crr}} \code{\link{rms}} \code{\link[rms]{nomogram}} \code{\link{plot.summary.rms}} \code{\link{axis}} \code{\link{pretty}} \code{\link{approx}}
##'
##' @keywords graphics
##'
##' @examples
##' data(prostate.dat)
##' dd <- datadist(prostate.dat)
##' options(datadist = "dd")
##' prostate.f <- cph(Surv(TIME_EVENT,EVENT_DOD == 1) ~ TX  + rcs(PSA,3) +
##'            BX_GLSN_CAT +  CLIN_STG + rcs(AGE,3) +
##'            RACE_AA, data = prostate.dat,
##'            x = TRUE, y= TRUE, surv=TRUE,time.inc = 144)
##' ## make a cph nomogram
##' nomogram.mk6(prostate.f, survtime=120, lp=FALSE,
##' funlabel = "Predicted 10-year cumulative incidence")
##'
nomogram.mk6 <-
  function(
      fit,
      ...,
      adj.to,
      lp = TRUE,
      lp.at,
      lplabel = "Linear Predictor",
      fun,
      fun.at,
      fun.lp.at,
      funlabel = "Predicted Value",
      fun.side,
      interact = NULL,
      intercept = 1,
      conf.int = FALSE,
      col.conf = c(1, 12),
      conf.space = c(0.08, 0.2),
      conf.lp = c("representative", "all", "none"),
      est.all = TRUE,
      abbrev = FALSE,
      minlength = 4,
      maxscale = 100,
      nint = 10,
      label.every = 1,
      force.label = FALSE,
      xfrac = 0.35,
      cex.axis = 0.85,
      cex.var = 1,
      col.grid = NULL,
      vnames = c("labels", "names"),
      varname.label = TRUE,
      varname.label.sep = "=",
      ia.space = 0.7,
      tck = NA,
      tcl = -0.25,
      lmgp = 0.4,
      omit = NULL,
      naxes,
      points.label = "Points",
      total.points.label = "Total Points",
      total.sep.page = FALSE,
      total.fun,
      verbose = FALSE,
      cap.labels = FALSE,
      total.min,
      total.max,
      survtime,
      mikeomit = NULL
  ) {

    if (!missing(survtime)) {
      # lp <- FALSE
      surv <- Survival(fit)
      fun <- function(x) surv(survtime, x)
      assign("surv", surv)
      assign("survtime", survtime)
    }
    if (any(stats::na.omit(match(attr(fit, "class"), "lrm"))) &
        missing(fun)) {
      # lp <- FALSE
      fun <- function(x) 1 / (1 + exp(-x))
    }
    ols <- any(stats::na.omit(match(attr(fit, "class"), "ols")))
    # if (!missing(total.min) && !missing(total.max) && lp && !ols)
    # stop("Total points control combined with lp is not setup yet.",
    #       "Pick one.")

    conf.lp <- match.arg(conf.lp)
    vnames <- match.arg(vnames)
    abb <- (is.logical(abbrev) && abbrev) || is.character(abbrev)
    if (is.logical(conf.int) && conf.int) {
      conf.int <- c(0.7, 0.9)
    }
    if (!is.logical(conf.int) && (length(conf.int) != length(col.conf))) {
      stop("conf.int and col.conf must have same length")
    }
    oldpar <- Hmisc::oPar()
    mgp <- oldpar$mgp
    mar <- oldpar$mar
    graphics::par(mgp = c(mgp[1], lmgp, mgp[3]), mar = c(
      mar[1], 1.1, mar[3],
      mar[4]
    ))
    on.exit(Hmisc::setParNro(oldpar))
    tck2 <- tck / 2
    tcl2 <- tcl / 2
    tck3 <- tck / 3
    tcl3 <- tcl / 3
    se <- FALSE
    if (any(conf.int > 0)) {
      se <- TRUE
      zcrit <- stats::qnorm((conf.int + 1) / 2)
      bar <- function(x, y, zcrit, se, col.conf, nlev = 4) {
        y <- rep(seq(y[1], y[2], length = nlev),length.out = length(x))
        for (j in seq(1,length(x))) {
          xj <- x[j]
          yj <- y[j]
          W <- c(0, zcrit) * se[j]
          for (i in seq(1,length(zcrit))) {
            graphics::segments(
              xj - W[i + 1], yj, xj - W[i], yj,
              col = col.conf[i], lwd = 1
            )
            graphics::segments(
              xj + W[i + 1], yj, xj + W[i], yj,
              col = col.conf[i], lwd = 1
            )
          }
        }
      }
    }
    nfun <- if (missing(fun)) {
      0
    } else if (is.list(fun)) {
      length(fun)
    } else {
      1
    }
    if (nfun > 1 && length(funlabel) == 1) {
      funlabel <- rep(funlabel, nfun)
    }
    if (nfun > 0 && is.list(fun) && length(names(fun))) {
      funlabel <- names(fun)
    }
    if (!missing(fun.at)) {
      if (!is.list(fun.at)) {
        fun.at <- rep(list(fun.at), nfun)
      }
    }
    if (!missing(fun.lp.at)) {
      if (!is.list(fun.lp.at)) {
        fun.lp.at <- rep(list(fun.lp.at), nfun)
      }
    }
    if (!missing(fun.side)) {
      if (!is.list(fun.side)) {
        fun.side <- rep(list(fun.side), nfun)
      }
      if (any(!(unlist(fun.side) %in% c(1, 3)))) {
        stop("fun.side must contain only the numbers 1 and 3")
      }
    }
    at <- fit$Design
    if (!length(at)) {
      at <- getOldDesign(fit)
    }
    assume <- at$assume.code
    if (any(assume == 10)) {
      stop("does not currently work with matrix factors in model")
    }
    name <- at$name
    names(assume) <- name
    parms <- at$parms
    label <- if (vnames == "labels") {
      at$label
    } else {
      name
    }
    if (any(d <- duplicated(name))) {
      stop("duplicated variable names:", name[d])
    }
    label <- name
    if (vnames == "labels") {
      label <- at$label
      if (any(d <- duplicated(label))) {
        stop("duplicated variable labels:", label[d])
      }
    }
    ia <- at$interactions
    factors <- list(...)
    nf <- length(factors)
    which <- if (est.all) {
      seq(1,length(assume))[assume != 8]
    } else {
      seq(1,length(assume))[assume != 8 & assume != 9]
    }
    if (nf > 0) {
      jw <- charmatch(names(factors), name, 0)
      if (any(jw == 0)) {
        stop(
          "factor name(s) not in the design:",
          names(factors)[jw == 0])
      }
      if (!est.all) {
        which <- jw
      }
    }
    Limval <- Getlim(at, allow.null = TRUE, need.all = FALSE)
    values <- Limval$values
    lims <- Limval$limits[c(6, 2, 7), , drop = FALSE]
    lims <- unclass(lims)
    for (i in seq(1,length(lims))) {
      if (is.factor(lims[[i]])) {
        lims[[i]] <- as.character(lims[[i]])
      }
    }
    attr(lims, "class") <- "data.frame"
    ucat <- rep(FALSE, length(assume))
    names(ucat) <- name
    for (i in seq(1,length(assume))[assume != 5 & assume < 8]) {
      ucat[i] <- !is.null(V <- values[[name[i]]])
      if (ucat[i]) {
        parms[[name[i]]] <- V
      }
    }
    discrete <- assume == 5 | assume == 8 | ucat
    names(discrete) <- name
    nrp <- Hmisc::num.intercepts(fit)
    Intercept <- if (nrp > 0) {
      fit$coefficients[intercept]
    } else if (!is.null(fit$center)) {
      -fit$center
    } else {
      0
    }
    intercept.offset <- if (nrp < 2) {
      0
    } else {
      fit$coefficients[intercept] - fit$coefficients[1]
    }
    settings <- list()
    for (i in which[assume[which] < 9]) {
      ni <- name[i]
      z <- factors[[ni]]
      lz <- length(z)
      if (lz < 2) {
        settings[[ni]] <- value.chk(
          at, i, NA, -nint, Limval,
          type.range = "full"
        )
      } else if (lz > 0 && any(is.na(z))) {
        stop("may not specify NA as a variable value")
      }
      if (lz == 1) {
        lims[2, i] <- z
      } else if (lz > 1) {
        settings[[ni]] <- z
        if (is.null(lims[[ni]]) || is.na(lims[2, ni])) {
          lims[[ni]] <- c(NA, z[1], NA)
          stop(
            "adjustment values for ", ni,
            " not defined in datadist;",
            "taken to be first value specified (",
            z[1], ")",
          )
        }
      }
    }
    adj <- lims[2, , drop = FALSE]
    if (!missing(adj.to)) {
      for (nn in names(adj.to)) adj[[nn]] <- adj.to[[nn]]
    }
    isna <- vapply(adj, is.na, FALSE)
    if (any(isna)) {
      stop(
        "adjustment values not defined here or with datadist for",
        name[assume != 9][isna]
      )
    }
    num.lines <- 0
    space.used <- 0
    entities <- 0
    set <- list()
    nset <- character(0)
    iset <- 0
    start <- len <- NULL
    end <- 0
    main.effects <- which[assume[which] < 8]
    if (any(assume == 9)) {
      main.effects <- main.effects[
        order(10 * discrete[main.effects] + (
          name[main.effects] %in% names(interact)))]
    }
    rel <- related.predictors(at)
    already.done <- structure(rep(FALSE, length(name)), names = name)
    for (i in main.effects) {
      nam <- name[i]
      ## changhong add
      if (already.done[nam] || (nam %in% omit) || (nam %in% mikeomit)) {
        next
      }
      r <- if (length(rel[[nam]])) {
        sort(rel[[nam]])
      } else {
        NULL
      }
      if (length(r) == 0) {
        num.lines <- num.lines + 1
        space.used <- space.used + 1
        entities <- entities + 1
        x <- list()
        x[[nam]] <- settings[[nam]]
        iset <- iset + 1
        attr(x, "info") <- list(
          predictor = nam, effect.name = nam,
          type = "main"
        )
        set[[iset]] <- x
        nset <- c(nset, label[i])
        start <- c(start, end + 1)
        n <- length(settings[[nam]])
        len <- c(len, n)
        end <- end + n
        NULL
      } else {
        namo <- name[r]
        s <- !(name[r] %in% names(interact))
        if (any(s)) {
          if (is.null(interact)) {
            interact <- list()
          }
          for (j in r[s]) {
            nj <- name[j]
            if (discrete[j]) {
              interact[[nj]] <- parms[[nj]]
            }
            NULL
          }
          s <- !(name[r] %in% names(interact))
        }
        if (any(s)) {
          stop(
            "factors not defined in interact=list(...):",
            name[r[s]]
          )
        }
        combo <- expand.grid(interact[namo])
        oldClass(combo) <- NULL
        acombo <- combo
        if (abb) {
          for (n in if (is.character(abbrev)) {
            abbrev
          } else {
            names(acombo)
          }) {
            if (discrete[n]) {
              acombo[[n]] <- abbreviate(
                parms[[n]],
                minlength = if (minlength ==1) {
                  4
                } else {
                  minlength
                })[combo[[n]]]
            }
          }
        }
        for (n in names(combo)) {
          if (is.factor(combo[[n]])) {
            combo[[n]] <- as.character(combo[[n]])
            acombo[[n]] <- as.character(acombo[[n]])
            NULL
          }
        }
        entities <- entities + 1
        already.done[namo] <- TRUE
        for (k in seq(1,length(combo[[1]]))) {
          num.lines <- num.lines + 1
          space.used <- space.used + if (k == 1) {
            1
          } else {
            ia.space
          }
          x <- list()
          x[[nam]] <- settings[[nam]]
          for (nm in namo) x[[nm]] <- combo[[nm]][k]
          iset <- iset + 1
          set.name <- paste(nam, " (", sep = "")
          for (j in seq(1,length(acombo))) {
            set.name <- paste(set.name, if (varname.label) {
              paste(namo[j], varname.label.sep, sep = "")
            } else {
              ""
            }, format(acombo[[j]][k]), sep = "")
            if (j < length(acombo)) {
              set.name <- paste(set.name, " ", sep = "")
            }
          }
          set.name <- paste(set.name, ")", sep = "")
          ia.names <- NULL
          for (j in r) {
            ia.names <- c(ia.names, name[interactions.containing(
              at,
              j
            )])
          }
          ia.names <- unique(ia.names)
          attr(x, "info") <- list(predictor = nam, effect.name = c(
            nam,
            namo[assume[namo] != 8], ia.names
          ),
          type = if (k == 1) {
            "first"
          } else {
            "continuation"
          })
          set[[iset]] <- x
          nset <- c(nset, set.name)
          start <- c(start, end + 1)
          n <- length(settings[[nam]])
          len <- c(len, n)
          end <- end + n
          NULL
        }
        NULL
      }
    }
    xadj <- unclass(Design.levels(adj, at))
    for (k in seq(1,length(xadj))) xadj[[k]] <- rep(xadj[[k]], sum(len))
    j <- 0
    for (S in set) {
      j <- j + 1
      ns <- names(S)
      nam <- names(S)
      for (k in seq(1,length(nam))) {
        xadj[[nam[k]]][start[j]:(start[j] + len[j] - 1)] <- S[[k]]
        NULL
      }
    }
    xadj <- structure(
      xadj, class = "data.frame",
      row.names = as.character(seq(1,sum(len))))
    xx <- predictDesign(
      fit,
      newdata = xadj, type = "terms",
      center.terms = FALSE, se.fit = FALSE, kint = intercept
    )
    if (any(is.infinite(xx))) {
      stop(
        "variable limits and transformations are such that an",
        "infinite axis value has resulted.\nRe-run specifying",
        "your own limits to variables.")
    }
    if (se) {
      xse <- predictDesign(
        fit,
        newdata = xadj, se.fit = TRUE,
        kint = intercept
      )
    }
    R <- matrix(NA, nrow = 2, ncol = length(main.effects), dimnames = list(
      NULL,
      name[main.effects]
    ))
    R[1, ] <- 1e+30
    R[2, ] <- -1e+30
    for (i in seq(1,num.lines)) {
      is <- start[i]
      ie <- is + len[i] - 1
      s <- set[[i]]
      setinfo <- attr(s, "info")
      nam <- setinfo$effect.name
      xt <- xx[is:ie, nam]
      if (length(nam) > 1) {
        xt <- apply(xt, 1, sum)
      }
      set[[i]]$Xbeta <- xt
      r <- range(xt)
      pname <- setinfo$predictor
      R[1, pname] <- min(R[1, pname], r[1])
      R[2, pname] <- max(R[2, pname], r[2])
      if (se) {
        set[[i]]$Xbeta.whole <- xse$linear.predictors[is:ie]
        set[[i]]$se.fit <- xse$se.fit[is:ie]
        NULL
      }
      NULL
    }
    R <- R[, R[1, ] < 1e+30, drop = FALSE]
    sc <- maxscale / max(R[2, ] - R[1, ])
    xl <- -xfrac * maxscale
    Intercept <- Intercept + sum(R[1, ])
    if (missing(naxes)) {
      naxes <- if (total.sep.page) {
        max(space.used + 1, nfun + lp + 1)
      } else {
        space.used + 1 + nfun + lp + 1
      }
    }
    Format <- function(x) {
      f <- character(l <- length(x))
      for (i in seq(1,l)) f[i] <- format(x[i])
      f
    }
    newpage <- function(
    naxes, xl, maxscale, cex.var, nint, space.used,
    col.grid, cex.axis, tck, tck2, tcl, tcl2,
    label.every, force.label, points = TRUE,
    points.label = "Points", usr) {
      y <- naxes - 1
      plot(
        0, 0,
        xlim = c(xl, maxscale), ylim = c(0, y), type = "n",
        axes = FALSE, xlab = "", ylab = ""
      )
      if (!missing(usr)) {
        graphics::par(usr = usr)
      }
      if (!points) {
        return(y + 1)
      }
      ax <- c(0, maxscale)
      graphics::text(xl, y, points.label, adj = 0, cex = cex.var)
      x <- pretty(ax, n = nint)
      dif <- x[2] - x[1]
      x2 <- seq((x[1] + x[2]) / 2, max(x), by = dif)
      x2 <- sort(c(x2 - dif / 4, x2, x2 + dif / 4))
      if (length(col.grid)) {
        graphics::segments(
          x, y, x, y - space.used,
          col = col.grid[1],
          lwd = 1
        )
        graphics::segments(
          x2, y, x2, y - space.used,
          col = col.grid[-1],
          lwd = 1
        )
      }
      axisf(
        3,
        at = x, pos = y, cex = cex.axis, tck = tck,
        tcl = tcl, label.every = label.every,
        force.label = force.label, padj = 0
      )
      axisf(
        3,
        at = x2, labels = FALSE, pos = y, tck = tck2,
        tcl = tcl2, cex = cex.axis
      )
      y
    }
    y <- newpage(
      naxes, xl, maxscale, cex.var, nint, space.used,
      col.grid, cex.axis, tck, tck2, tcl, tcl2,
      label.every = label.every,
      force.label = force.label, points.label = points.label
    )
    i <- 0
    Abbrev <- list()
    for (S in set) {
      i <- i + 1
      setinfo <- attr(S, "info")
      type <- setinfo$type
      y <- y - (if (type == "continuation") {
        ia.space
      } else {
        1
      })
      if (y < -0.05) {
        y <- newpage(
          naxes, xl, maxscale, cex.var, nint,
          space.used, col.grid, cex.axis, tck, tck2, tcl,
          tcl2,
          label.every = label.every, force.label = force.label,
          points.label = points.label
        ) - (if (type == "continuation") {
          ia.space
        } else {
          1
        })
      }
      graphics::text(
        xl, y, paste(
          Hmisc::strgraphwrap(
            nset[[i]],
            abs(xl), cex = cex.var),
          collapse = "\n"
        ), adj = 0, cex = cex.var)
      x <- S[[1]]
      nam <- names(S)[1]
      fx <- if (is.character(x)) {
        x
      } else {
        Hmisc::sedit(Format(x), " ", "")
      }
      if (abb && discrete[nam] &&
          (is.logical(abbrev) || nam %in% abbrev)) {
        old.text <- fx
        fx <- if (abb && minlength == 1) {
          letters[seq(1,length(fx))]
        } else {
          abbreviate(fx, minlength = minlength)
        }
        Abbrev[[nam]] <- list(abbrev = fx, full = old.text)
      }
      j <- match(nam, name, 0)
      #if (any(j == 0)) stop("program logic error 1")
      is <- start[i]
      ie <- is + len[i] - 1
      xt <- (S$Xbeta - R[1, nam]) * sc
      set[[i]]$points <- xt
      r <- rle(xt)
      if (any(r$length > 1)) {
        is <- 1
        for (j in r$length) {
          ie <- is + j - 1
          if (j > 1) {
            fx[ie] <- if (discrete[nam] || ie < length(xt)) {
              paste(fx[is], "-", fx[ie], sep = "")
            } else {
              paste(fx[is], "+", sep = "")
            }
            fx[is:(ie - 1)] <- ""
            xt[is:(ie - 1)] <- NA
          }
          is <- ie + 1
        }
        fx <- fx[!is.na(xt)]
        xt <- xt[!is.na(xt)]
      }
      side <- c(1, 3)
      padj <- c(1, 0)
      new.mgp <- vector(mode = "list", 2)
      new.mgp[[2]] <- c(0, lmgp, 0)
      new.mgp[[1]] <- new.mgp[[2]] - c(0, 0.6, 0)
      ch <- if (length(xt) > 2) {
        c(FALSE, FALSE, diff(diff(xt) > 0) != 0)
      } else {
        rep(FALSE, length(xt))
      }
      if (discrete[nam] && length(xt) > 1) {
        j <- order(xt)
        graphics::lines(range(xt), rep(y, 2))
        for (k in c(1,2)) {
          is <- j[seq(k, length(j), by = 2)]
          new.labs <- if (cap.labels) {
            Hmisc::capitalize(fx[is])
          } else {
            fx[is]
          }
          axisf(
            side[k],
            at = xt[is], labels = new.labs,
            pos = y, cex = cex.axis, tck = tck, tcl = tcl,
            force.label = force.label ||
              (abb && minlength == 1 && (
                is.logical(abbrev) || nam %in% abbrev)),
            disc = TRUE, mgp = new.mgp[[k]], padj = padj[k]
          )
          if (se) {
            bar(
              xt[is], if (k == 1) {
                y - conf.space - 0.32
              } else {
                y + conf.space + 0.32
              }, zcrit, sc * S$se.fit[is],
              col.conf
            )
          }
        }
      } else if (!any(ch)) {
        axisf(
          1,
          at = xt, labels = fx, pos = y, cex = cex.axis,
          tck = tck, tcl = tcl, mgp = new.mgp[[1]],
          label.every = label.every, force.label = force.label,
          disc = discrete[nam], padj = padj[1]
        )
        if (se) {
          bar(
            xt, y + conf.space, zcrit, sc * S$se.fit,
            col.conf
          )
        }
      } else {
        graphics::lines(range(xt), rep(y, 2))
        j <- seq(1,length(ch))[ch]
        if (max(j) < length(ch)) {
          j <- c(j, length(ch) + 1)
        }
        flag <- 1
        is <- 1
        for (k in j) {
          ie <- k - 1
          axisf(
            side[flag],
            at = xt[is:ie], labels = fx[is:ie],
            pos = y, cex = cex.axis, tck = tck, tcl = tcl,
            label.every = label.every, force.label = force.label,
            mgp = new.mgp[[flag]], disc = discrete[nam],
            padj = padj[flag]
          )
          if (se) {
            bar(
              xt[is:ie], if (side[flag] == 1) {
                y - conf.space - 0.32
              } else {
                y + conf.space + 0.32
              }, zcrit, sc * S$se.fit[is:ie],
              col.conf
            )
          }
          flag <- if (flag == 2) {
            1
          } else {
            2
          }
          is <- ie + 1
        }
      }
    }
    # browser()
    if (missing(lp.at)) {
      xb <- fit$linear.predictors
      if (!length(xb)) {
        xb <- fit$fitted.values
      }
      if (!length(xb)) {
        xb <- fit$fitted
      }
      if (!length(xb)) {
        stop(
          "lp.at not given and fit did not store",
          "linear.predictors or fitted.values")
      }
      if (nrp > 1) {
        xb <- xb + intercept.offset
      }
      lp.at <- pretty(range(xb), n = nint)
    }
    sum.max <- if (entities == 1) {
      maxscale
    } else {
      max(maxscale, sc * max(lp.at - Intercept))
    }
    x <- pretty(c(0, sum.max), n = nint)
    new.max <- max(x)
    xl.old <- xl
    xl <- -xfrac * new.max
    u <- graphics::par()$usr
    if (!missing(total.fun)) {
      total.fun()
    }
    usr <- c(xl * u[1] / xl.old, new.max * u[2] / maxscale, u[3:4])
    graphics::par(usr = usr)
    x.double <- seq(x[1], new.max, by = (x[2] - x[1]) / 5)
    y <- y - 1
    if (y < -0.05 || total.sep.page) {
      y <- newpage(
        naxes, xl, maxscale, cex.var, nint, space.used,
        col.grid, cex.axis, tck, tck2, tcl, tcl2,
        label.every = label.every,
        force.label = force.label, points = FALSE, usr = usr
      ) -
        1
    }
    graphics::text(xl, y, total.points.label, adj = 0, cex = cex.var)
    # changhong edit here
    if (missing(total.max) || missing(total.min)) {
      axisf(
        1,
        at = x, pos = y, cex = cex.axis, tck = tck, tcl = tcl,
        label.every = label.every, force.label = force.label,
        mgp = c(0, lmgp - 0.6, 0), padj = 1
      )
      axisf(
        1,
        at = x.double, labels = FALSE, pos = y, tck = tck2,
        tcl = tcl2, cex = cex.axis
      )
    } else {
      axisf(
        1,
        at = x, pos = y,
        labels = round(seq(
          total.min,total.max,
          by = (total.max - total.min) / (length(x) - 1)
        )),
        cex = cex.axis, tck = tck, tcl = tcl,
        label.every = label.every, force.label = force.label,
        mgp = c(0, lmgp - 0.6, 0), padj = 1
      )
      axisf(
        1,
        at = x.double, labels = FALSE, pos = y, tck = tck2,
        tcl = tcl2, cex = cex.axis
      )
    }
    iset <- iset + 1
    nset <- c(nset, "total.points")
    set[[iset]] <- list(x = x, y = y)
    # browser()
    if (lp) {
      x2 <- seq(lp.at[1], max(lp.at), by = (lp.at[2] - lp.at[1]) / 2)
      scaled.x <- (lp.at - Intercept) * sc
      scaled.x2 <- (x2 - Intercept) * sc
      y <- y - 1
      if (y < -0.05) {
        y <- newpage(
          naxes, xl, maxscale, cex.var, nint,
          space.used, col.grid, cex.axis, tck, tck2, tcl,
          tcl2,
          label.every = label.every, force.label = force.label,
          points = FALSE, usr = usr
        ) - 1
      }
      graphics::text(xl, y, lplabel, adj = 0, cex = cex.var)
      # changhong edit here on 05/12/2010
      if (missing(total.max) || missing(total.min)) {
        axisf(
          1,
          at = scaled.x, labels = Format(lp.at), pos = y,
          cex = cex.axis, tck = tck, tcl = tcl,
          label.every = label.every,
          force.label = force.label, mgp = c(
            0, lmgp - 0.6,
            0
          ), padj = 1
        )
      } else {
        axisf(
          1,
          at = ((scaled.x - total.min) * (
            max(x) - min(x))) / (total.max - total.min),
          labels = Format(lp.at), pos = y,
          cex = cex.axis, tck = tck, tcl = tcl,
          label.every = label.every,
          force.label = force.label, mgp = c(
            0, lmgp - 0.6,
            0
          ), padj = 1
        )
      }
      axisf(
        1,
        at = scaled.x2, labels = FALSE, tck = tck2,
        tcl = tcl2, pos = y, cex = cex.axis
      )
      iset <- iset + 1
      nset <- c(nset, "lp")
      set[[iset]] <- list(x = scaled.x, y = y, x.real = lp.at)
      if (se && conf.lp != "none") {
        xxb <- NULL
        xse <- NULL
        for (S in set) {
          xxb <- c(xxb, S$Xbeta.whole)
          xse <- c(xse, S$se.fit)
        }
        i <- order(xxb)
        if (length(xxb) < 16 | conf.lp == "representative") {
          nlev <- 4
          w <- 1
        } else {
          nlev <- 8
          w <- 2
        }
        if (conf.lp == "representative") {
          deciles <- Hmisc::cut2(xxb[i], g = 10)
          mean.xxb <- tapply(xxb[i], deciles, mean)
          median.se <- tapply(xse[i], deciles, stats::median)
          bar((mean.xxb - Intercept) * sc, y + c(
            conf.space[1],
            conf.space[1] + w * diff(conf.space)
          ), zcrit,
          sc * median.se, col.conf,
          nlev = nlev
          )
        } else {
          bar((xxb[i] - Intercept) * sc, y + c(
            conf.space[1],
            conf.space[1] + w * diff(conf.space)
          ), zcrit,
          sc * xse[i], col.conf,
          nlev = nlev
          )
        }
      }
    }
    # browser()
    if (nfun > 0) {
      if (!is.list(fun)) {
        fun <- list(fun)
      }
      i <- 0
      for (func in fun) {
        i <- i + 1
        if (!missing(fun.lp.at)) {
          xseq <- fun.lp.at[[i]]
          fat <- func(xseq)
          w <- xseq
        } else {
          if (missing(fun.at)) {
            fat <- pretty(func(range(lp.at)), n = nint)
          } else {
            fat <- fun.at[[i]]
          }
          if (verbose) {
            cat(
              "Function", i,
              "values at which to place tick marks:\n")
            print(fat)
          }
          xseq <- seq(min(lp.at), max(lp.at), length = 1000)
          fu <- func(xseq)
          s <- !is.na(fu)
          w <- stats::approx(fu[s], xseq[s], fat, ties = mean)$y
          if (verbose) {
            cat("Estimated inverse function values (lp):\n")
            print(w)
          }
        }
        s <- !(is.na(w) | is.na(fat))
        w <- w[s]
        fat <- fat[s]
        fat.orig <- fat
        fat <- if (is.category(fat)) {
          as.character(fat)
        } else {
          Format(fat)
        }
        scaled <- (w - Intercept) * sc
        y <- y - 1
        if (y < -0.05) {
          y <- newpage(
            naxes, xl, maxscale, cex.var, nint,
            space.used, col.grid, cex.axis, tck, tck2,
            tcl, tcl2,
            label.every = label.every, force.label = force.label,
            points = FALSE, usr = usr
          ) - 1
        }
        graphics::text(xl, y, funlabel[i], adj = 0, cex = cex.var)
        sides <- if (missing(fun.side)) {
          rep(1, length(fat))
        } else {
          (fun.side[[i]])[s]
        }
        if (length(sides) != length(fat)) {
          stop(
            "fun.side vector not same length as",
            "fun.at or fun.lp.at")
        }

        # changhong edit here
        for (jj in seq(1,length(fat))) {
          if (missing(total.max) || missing(total.min)) {
            graphics::axis(
              sides[jj],
              at = scaled[jj], label = fat[jj],
              pos = y, cex.axis = cex.axis, tck = tck, tcl = tcl,
              mgp = if (sides[jj] == 1) {
                c(0, lmgp - 0.6, 0)
              } else {
                c(0, lmgp, 0)
              }, padj = if (sides[jj] == 1) {
                1
              } else {
                0
              }
            )
            graphics::lines(range(scaled), rep(y, 2))
          } else {
            graphics::axis(
              sides[jj],
              at = ((
                scaled[jj] - total.min) * (
                  max(x) - min(x))) /
                (total.max - total.min),
              label = fat[jj],
              pos = y, cex.axis = cex.axis,
              tck = tck, tcl = tcl,
              mgp = if (sides[jj] == 1) {
                c(0, lmgp - 0.6, 0)
              } else {
                c(0, lmgp, 0)
              }, padj = if (sides[jj] == 1) {
                1
              } else {
                0
              }
            )
            graphics::lines(((
              range(scaled) - total.min) *
                (max(x) - min(x))) /
                (total.max - total.min), rep(y, 2))
          }
        }

        iset <- iset + 1
        nset <- c(nset, funlabel[i])
        set[[iset]] <- list(x = scaled, y = y, x.real = fat.orig)
      }
    }
    names(set) <- nset
    set$abbrev <- Abbrev
    oldClass(set) <- "nomogram"
    invisible(set)

  }

nomo2.crr <-
  function(x, f.crr, time) {
    assign("f.crr", f.crr)
    assign("time", time)
    vapply(x, function(x) pred2.crr(f.crr, x, time), 0.5)
  }

is.category <- function (x)
  length(attr(x, "levels")) > 0 && mode(x) == "numeric"

getOldDesign <- function(fit) {
  at <- attr(fit$terms,'Design')
  if(is.null(at))
    stop('fit was not created by a Design library fitting function')
  at
}

Design.levels <- function(df, at) {
  ac <- at$assume.code
  for (nn in names(df)) {
    j <- match(nn, at$name, 0)
    if (j > 0) {
      if ((ac[j] == 5 | ac[j] == 8) & length(lev <- at$parms[[nn]])) {
        df[[nn]] <- factor(df[[nn]], lev)
      }
    }
  }
  df
}

value.chk <- rms:::value.chk
axisf <- rms:::axisf