R/plot.tvcure.R

In tvcure: Additive Cure Survival Model with Time-Varying Covariates

#' Plot visual information related to a tvcure object.
#' @description Visualization of the estimated additive terms and of the reference (cumulative) hazard function in a tvcure object.
#'
#' @usage \method{plot}{tvcure}(x, ngrid=300, ci.level=.95, pages=0, select=NULL,
#'                              fill=TRUE, pointwise=TRUE, mar=c(4,5,1,1),
#'                              xlim0=NULL, ylim0=NULL,
#'                              xlim1=NULL, ylim1=NULL, xlim2=NULL, ylim2=NULL,
#'                              equal.ylims=TRUE,...)
#'
#' @param x a \code{\link{tvcure.object}}.
#' @param ngrid (optional) number of points used to plot the fitted additive terms. (Default: 300).
#' @param ci.level (optional) nominal level for the plotted pointwise credible intervals. (Default: 0.95).
#' @param pages The number of pages over which to spread the output. For example, if pages=1 then all terms will be plotted on one page with the layout performed automatically. Set to 0 to have the routine leave all graphics settings as they are. (Default 0).
#' @param select Allows the plot for a single model term to be selected for printing. e.g. if you just want the plot for the second smooth term set select=2. The plot of the reference hazard \eqn{f_0(t)} and cumulative hazard \eqn{F_0(t)} functions is provided when select=0. When select=-1, only the reference hazard is plotted.(Default: NULL).
#' @param fill Logical indicating whether credible regions should be greyed out. (Default: TRUE).
#' @param pointwise Logical indicating whether only pointwise credible intervals should be plotted for additive terms. When FALSE, simultaneous credible regions are also provided. (Default: TRUE).
#' @param mar A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot. (Default: c(4,5,1,1)).
#' @param xlim0 Vector of length 2 specifying x-axis limits when plotting the estimated reference hazard function \eqn{\exp(\beta_0)f_0(t)}.
#' @param ylim0 Vector of length 2 specifying y-axis limits when plotting the estimated reference hazard function \eqn{\exp(\beta_0)f_0(t)}.
#' @param xlim1 Vector of length 2 specifying (common) x-axis limits when plotting the fitted additive term(s) in the long-term survival submodel. (Default: NULL).
#' @param ylim1 Vector of length 2 specifying (common) y-axis limits when plotting the fitted additive term(s) in the long-term survival submodel. (Default: NULL).
#' @param xlim2 Vector of length 2 specifying (common) x-axis limits when plotting the fitted additive term(s) in the short-term survival submodel. (Default: NULL).
#' @param ylim2 Vector of length 2 specifying (common) y-axis limits when plotting the fitted additive term(s) in the short-term survival submodel. (Default: NULL).
#' @param equal.ylims logical indicating if the same y-limits must be used when plotting the fitted additive terms from the same submodel. It can be overriden by non-NULL values for \code{ylim1} or \code{ylim2}. (Default: TRUE).
#' @param ... additional generic plotting arguments.
#'
#' @details Plot of the fitted additive terms, as well as of the reference hazard \eqn{f_0(t)} and cumulative hazard \eqn{F_0(t)} functions of the fitted tvcure model in \code{x}.
#'
#' @return In addition to the plots, an invisible list generated by the \code{\link{additive.tvcure}} function is returned.
#'
#' @author Philippe Lambert \email{p.lambert@uliege.be} based on the plot.gam function in mgcv for the \code{pages} argument.
#' @references Lambert, P. and Kreyenfeld, M. (2025).
#' Time-varying exogenous covariates with frequently changing values in double additive cure survival model: an application to fertility.
#' \emph{Journal of the Royal Statistical Society, Series A}. <doi:10.1093/jrsssa/qnaf035>
#'
#' @examples
#' \donttest{
#' require(tvcure)
#' ## Simulated data generation
#' beta = c(beta0=.4, beta1=-.2, beta2=.15) ; gam = c(gam1=.2, gam2=.2)
#' data = simulateTVcureData(n=500, seed=123, beta=beta, gam=gam,
#'                           RC.dist="exponential",mu.cens=550)$rawdata
#' ## TVcure model fitting
#' tau.0 = 2.7 ; lambda1.0 = c(40,15) ; lambda2.0 = c(25,70) ## Optional
#' model = tvcure(~z1+z2+s(x1)+s(x2), ~z3+z4+s(x3)+s(x4), data=data,
#'                tau.0=tau.0, lambda1.0=lambda1.0, lambda2.0=lambda2.0)
#' plot(model,pages=1)
#' }
#'
#' @seealso \code{\link{tvcure}}, \code{\link{tvcure.object}}, \code{\link{print.tvcure}}
#'
#' @export
#'
plot.tvcure = function(x, ngrid=300, ci.level=.95, pages=0, select=NULL,
                       fill=TRUE, pointwise=TRUE, mar=c(4,5,1,1),
                       xlim0=NULL, ylim0=NULL, xlim1=NULL, ylim1=NULL, xlim2=NULL, ylim2=NULL,
                       equal.ylims=TRUE,...){
    oldpar <- par(no.readonly = TRUE)
    on.exit(par(oldpar))
    ##
    obj = x
    add.grid = TRUE
    ##
    ## Compute additive term + envelope
    ## --------------------------------
    fhat = additive.tvcure(obj,ngrid=ngrid,ci.level=ci.level)
    ##
    ## Plot organization (extracted from mgcv::plot.gam)
    ## -------------------------------------------------
    if (is.null(select)){
        n.plots = 2 + fhat$J1 + fhat$J2 ## Number of plots: f0, F0 + additive terms
    } else {
        n.plots = length(select)
        if (all(select == 0)) n.plots = 2
    }
    if (pages > n.plots) pages = n.plots
    if (pages < 0) pages = 0
    if (pages != 0) {
        ppp = n.plots%/%pages
        if (n.plots%%pages != 0) {
            ppp = ppp + 1
            while (ppp * (pages - 1) >= n.plots) pages = pages - 1
        }
        c = r = trunc(sqrt(ppp))
        if (c < 1) r = c = 1
        if (c * r < ppp) c = c + 1
        if (c * r < ppp) r = r + 1
        par(mfrow = c(r, c))
    }
    else {
        ppp = 1
    }
    if ((pages == 0 && prod(par("mfcol")) < n.plots && dev.interactive()) ||
        pages > 1 && dev.interactive())
        ask = TRUE
    else ask = FALSE
    if (!is.null(select) && (all(select!=0))) {
        ask = FALSE
    }
    if (ask) {
        oask = devAskNewPage(TRUE)
        on.exit(devAskNewPage(oask))
    }
    ## Plot baseline f0 & F0
    ## ---------------------
    if (is.null(xlim0)) xlim0 = attr(fhat$f0,"support")
    ## Plot only f0 if select=-1
    if (!is.null(select) && (all(select==-1))){
        par(mar=mar)
        beta0 = x$fit$beta[1,1]
        curve(exp(beta0)*fhat$f0(x),
              xlim=xlim0,ylim=ylim0,
              xlab="time",ylab=bquote(e^{beta[0]}*~f[0](t)), type="n",...)
        grid(lwd=.5,lty=1)
        curve(exp(beta0)*fhat$f0(x),lwd=1.5,
              xlim=attr(fhat$f0,"support"),ylim=ylim0, add=TRUE,
              xlab="time",ylab=bquote(e^{beta[0]}*~f[0](t)), ...)
    }
    ## Plot both if select=0 or NULL
    if (is.null(select) || (all(select==0))){
        par(mar=mar)
        beta0 = x$fit$beta[1,1]
        curve(exp(beta0)*fhat$f0(x),
              xlim=xlim0,ylim=ylim0,
              xlab="time",ylab=bquote(e^{beta[0]}*~f[0](t)), type="n",...)
        grid(lwd=.5,lty=1)
        curve(exp(beta0)*fhat$f0(x),lwd=1.5,
              xlim=attr(fhat$f0,"support"),ylim=ylim0, add=TRUE,
              xlab="time",ylab=bquote(e^{beta[0]}*~f[0](t)), ...)
        ## with(fhat, curve(f0,
        ##                  xlim=attr(f0,"support"),xlab="time",ylab=bquote(e^{beta[0]}*~f[0](t))))
        with(fhat, curve(F0, xlim=xlim0, type="n",las=1,
                         xlab="time",ylab=bquote(F[0](t))))
        grid(lwd=.5,lty=1,ny=0)
        abline(h=seq(0,1,by=.1),col="grey",lwd=.5)
        with(fhat, curve(F0, xlim=attr(F0,"support"), add=TRUE,lwd=1.5,
                         xlab="time",ylab=bquote(F[0](t))))
    }
    ## Plot Additive terms
    ## -------------------
    plotAdd = function(x,y,y2=NULL,col=1,colfill=c("#CCCCCC80","#E5E5E580"),las=1,...) {
        matplot(x, y,type="n",col=col,las=las,
                xlim=xlims,ylim=ylims,xlab=xlab,ylab=ylab,
                lwd=c(1.5,1,1),lty=c(1,2,2),...)
        grid(lwd=.5,lty=1)
        if (!fill){
            if (!is.null(y2)){
                matplot(x,y2,type="l",add=TRUE,col="grey",las=las,
                        xlim=xlims,ylim=ylims,xlab=xlab,ylab=ylab,
                        lwd=c(1.5,1,1),lty=c(1,3,3),...)
            }
            matplot(x,y,type="l",add=TRUE,col=col,las=las,
                    xlim=xlims,ylim=ylims,xlab=xlab,ylab=ylab,
                    lwd=c(1.5,1,1),lty=c(1,2,2),...)
        } else {
            if (!is.null(y2)){
                plotRegion(x,y2,add=TRUE,col=col,colfill=colfill[2],las=las,
                           xlim=xlims,ylim=ylims,
                           xlab=xlab,ylab=ylab,lwd=2,...)
            }
            plotRegion(x,y,add=TRUE,col=col,colfill=colfill[1],las=las,
                       xlim=xlims,ylim=ylims,
                       xlab=xlab,ylab=ylab,lwd=2,...)
        }
    }
    if (fhat$J1 > 0){
        xlims = ylims = NULL
        if (equal.ylims){
            if (pointwise){
                ylims =  range(lapply(fhat$f1.grid, function(x) range(x$y.mat)))
            } else {
                ylims =  range(lapply(fhat$f1.grid, function(x) range(x$y.mat2)))
            }
        }
        if (!is.null(ylim1)) ylims = ylim1
        if (!is.null(xlim1)) xlims = xlim1
        for (j in 1:fhat$J1){
            if ((is.null(select)) || (any(select == j))){
                par(mar=mar)
                xlab = names(fhat$f1.grid)[j]
                ylab = bquote('f'[.(j)]*(.(xlab)))
                if (!pointwise) {
                    with(fhat$f1.grid[[j]], plotAdd(x,y.mat,y.mat2,...))
                } else {
                    with(fhat$f1.grid[[j]], plotAdd(x,y.mat,...))
                }
                if (obj$regr1$has.ref[j]){
                    xt = obj$regr1$ref.values[j]
                    rug(xt,ticksize=-.015,lwd=1)
                    rug(xt,ticksize=.015,lwd=1)
                }
            }
        }
    }
    if (fhat$J2 > 0){
        xlims = ylims = NULL
        if (equal.ylims){
            if (pointwise){
                ylims =  range(lapply(fhat$f2.grid, function(x) range(x$y.mat)))
            } else {
                ylims =  range(lapply(fhat$f2.grid, function(x) range(x$y.mat2)))
            }
        }
        if (!is.null(ylim2)) ylims = ylim2
        if (!is.null(xlim2)) xlims = xlim2
        for (j in 1:fhat$J2){
            if ((is.null(select)) || (any(select == fhat$J1+j))){
                par(mar=mar)
                xlab = names(fhat$f2.grid)[j]
                ylab = bquote(tilde('f')[.(j)]*(.(xlab)))
                if (!pointwise) {
                    with(fhat$f2.grid[[j]], plotAdd(x,y.mat,y.mat2,...))
                } else {
                    with(fhat$f2.grid[[j]], plotAdd(x,y.mat,...))
                }
                if (obj$regr2$has.ref[j]){
                    xt = obj$regr2$ref.values[j]
                    rug(xt,ticksize=-.015,lwd=1)
                    rug(xt,ticksize=.015,lwd=1)
                }
            }
        }
    }
    return(invisible(fhat))
}