R/plot.p3state.R

In p3state.msm: Analyzing Survival Data from an Illness-Death Model

#' Plot method for p3state objects
#'
#' Produces various plots related to the \code{p3state} model, including transition probabilities,
#' marginal distributions, and bivariate distribution functions.
#'
#' @param x An object of class \code{p3state}.
#' @param plot.trans Character or logical; specifies which transition probabilities to plot.
#'   Options include \code{"P11"}, \code{"P12"}, \code{"P22"}, \code{"P23"}, or \code{"all"}.
#'   Defaults to \code{NULL} (no plot).
#' @param plot.marginal Logical; whether to plot marginal distributions. Defaults to \code{NULL}.
#' @param plot.bivariate Logical; whether to plot the bivariate distribution function. Defaults to \code{NULL}.
#' @param time1 Numeric; start time for the plot. Defaults to 0 if missing.
#' @param time2 Numeric; end time for the plot. Defaults to the maximum observed time.
#' @param xlab Character; label for the x-axis. Defaults to "Time" if missing.
#' @param ylab Character; label for the y-axis. Optional.
#' @param zlab Character; label for the z-axis in 3D plots. Optional.
#' @param col Color(s) used for plots. Optional.
#' @param col.biv Color specification for bivariate contour plots. Optional.
#' @param ... Additional graphical parameters passed to plotting functions.
#'
#' @return Invisibly returns \code{NULL}. Plots are drawn as side effects.
#'
#' @examples
#' \dontrun{
#' # Suponha que 'fit' é um objeto p3state ajustado
#' plot(fit, plot.trans = "P11")
#' plot(fit, plot.marginal = TRUE)
#' plot(fit, plot.bivariate = TRUE)
#' }
#'
#' @export
#' 
plot.p3state <- function(
  x,
  plot.trans=NULL,
  plot.marginal=NULL,
  plot.bivariate=NULL,
  time1,
  time2,
  xlab,
  ylab,
  zlab,
  col,
  col.biv=NULL,
  ...
) {
  if ( missing(x) ) stop("Argument 'x' is missing with no default");
  if ( !inherits(x, "p3state") ) stop("'x' must be of class 'p3state'");
  mydata <- x$datafr;
  xlab.aux <- FALSE;
  ylab.aux <- FALSE;
  pmar <- FALSE;
  pbiv <- FALSE;
  if ( missing(xlab) ) {
    xlab <- "Time";
    xlab.aux <- TRUE;
  }
  if ( missing(ylab) ) ylab.aux <- TRUE;
  if ( missing(plot.trans) ) plot.trans <- FALSE;
  if ( missing(plot.marginal) ) plot.marginal <- FALSE;
  if (missing( plot.bivariate) ) plot.bivariate <- FALSE;
  if ( missing(col.biv) ) col.biv <- FALSE;
  if ( missing(time1) ) time1 <- 0;
  if ( missing(time2) ) time2 <- max(mydata[, 1]);
  if (time1 < 0 | time2 < 0) stop("'time1' and 'time2' must be positive");
  if (time1 > time2) stop("Argument 'time1' cannot be greater then 'time2'");
  if (sum(mydata[, 2]) == 0) stop("This is not a multi-state model");
  if (sum(mydata[, 2] == 0 & mydata[, 5] == 1) > 0) {ntrans <- 3;}
  else {ntrans <- 2;}
  if (plot.bivariate == TRUE) pbiv <- TRUE;
  if (plot.marginal == TRUE) pmar <- TRUE;
  if (ntrans == 3) {
    plot.marginal <- FALSE;
    plot.bivariate <- FALSE;
  }
  if (plot.trans == "P11" | plot.trans == "all") {
    x <- which((mydata[, 2] == 1 | mydata[, 5] == 1) & mydata[, 1] > time1);
    y1 <- sort(mydata[x, 1]);
    y2 <- unique(y1);
    y <- c(time1, y2);
    mydata2 <- matrix( data=0, ncol=2, nrow=length(y) );
    for ( k in 1:length(y) ) {
      mydata2[k, 1] <- y[k];
    }
    for ( k in 1:length(y) ) {
      mydata2[k, 2] <- pLIDA(mydata, time1, mydata2[k, 1], tp="p11");
    }
    if (ylab.aux == TRUE) ylab <- paste("Estimated prob. p11(", time1, ",t)");
    matplot(mydata2[,1], mydata2[,2], xlab=xlab, ylab=ylab, type="s", ...);
    oask <- devAskNewPage(TRUE);
    on.exit( devAskNewPage(oask) );
  }
  if (plot.trans == "P12" | plot.trans == "all") {
    x <- which(mydata[, 1] > time1 & mydata[, 5] == 1);
    y1 <- c(mydata[x, 1] + 1e-09, mydata[x, 4] + 1e-09);
    y2 <- unique(y1);
    y <- sort(y2);
    y <- c(time1, y);
    mydata2 <- matrix( data=0, ncol=2, nrow=length(y) );
    for ( k in 1:length(y) ) {mydata2[k, 1] <- y[k];}
    for ( k in 1:length(y) ) {
      mydata2[k,2] <- pLIDA(mydata, time1, mydata2[k,1], tp="p12");
    }
    if (ylab.aux == TRUE) ylab <- paste("Estimated prob. p12(", time1, ",t)");
    matplot(mydata2[,1], mydata2[,2], xlab=xlab, ylab=ylab, type="s", ...);
    oask <- devAskNewPage(TRUE);
    on.exit( devAskNewPage(oask) );
  }
  if (plot.trans == "P22" | plot.trans == "P23" | plot.trans == "all") {
    x <- which(mydata[,5] == 1 & mydata[,1] <= time1);
    q2 <- c(mydata[x,1], mydata[x,4]);
    q3 <- q2[q2 >= time1];
    y1 <- sort(q3);
    y2 <- unique(y1);
    y <- c(time1, y2);
    mydata2 <- matrix( data=0, ncol=2, nrow=length(y) );
    for ( k in 1:length(y) ) {mydata2[k, 1] <- y[k];}
    for ( k in 1:length(y) ) {
      mydata2[k,2] <- pLIDA(mydata, time1, mydata2[k,1], tp="p22");
    }
    if (ylab.aux == TRUE) ylab <- paste("Estimated prob. p22(", time1, ",t)");
    if (plot.trans == "P22" | plot.trans == "all") 
      matplot(mydata2[,1], mydata2[,2], xlab=xlab, ylab=ylab, type="s", ...);
    oask <- devAskNewPage(TRUE);
    on.exit( devAskNewPage(oask) );
    if (ylab.aux == TRUE) ylab <- paste("Estimated prob. p23(", time1, ",t)");
    if (plot.trans == "P23" | plot.trans == "all") 
      matplot(mydata2[,1], 1 - mydata2[,2], xlab=xlab, ylab=ylab, type="s", ...);
    oask <- devAskNewPage(TRUE);
    on.exit( devAskNewPage(oask) );
  }
  if (plot.marginal == "TRUE") {
    x <- which(mydata[,5] == 1);
    y1 <- mydata[x,4];
    y2 <- unique(y1);
    y3 <- sort(y2);
    y <- c(0, y3);
    mydata2 <- matrix( data=0, ncol=2, nrow=length(y) );
    for ( k in 1:length(y) ) {mydata2[k, 1] <- y[k];}
    for ( k in 1:length(y) ) {
      mydata2[k,2] <- Biv(mydata, max(mydata[,1]), mydata2[k,1]);
    }
    oask <- devAskNewPage(TRUE);
    on.exit( devAskNewPage(oask) );
    if (ylab.aux == TRUE) ylab <- "Estimated marginal dist. F2(t)";
    matplot(mydata2[,1], mydata2[,2], xlab=xlab, ylab=ylab, type="s", ...);
    oask <- devAskNewPage(TRUE);
    on.exit( devAskNewPage(oask) );
  }
  if (plot.bivariate == "TRUE") {
    x <- which(mydata[,5] == 1);
    x1 <- unique( sort(mydata[x,1]) );
    y1 <- unique( sort(mydata[x,3]) );
    mydata2 <- matrix( data=0, nrow=length(x1), ncol=length(y1) );
    for ( k in 1:length(x1) ) {
      for ( j in 1:length(y1) ) {mydata2[k, j] <- Biv(mydata, x1[k], y1[j]);}
    }
    z <- mydata2;
    if (xlab.aux == TRUE) xlab <- "Time in state 1";
    if (ylab.aux == TRUE) ylab <- "Time in state 2";
    if ( missing(zlab) ) zlab <- "F(t1,t2)";
    if ( missing(col) ) col <- "lightblue";
    persp(x1, y1, z, theta=30, phi=30, ticktype="detailed", expand=0.5,
      xlab=xlab, ylab=ylab, zlab=zlab, col=col, shade=0.2, ...);
    oask <- devAskNewPage(TRUE);
    on.exit( devAskNewPage(oask) );
    x1 <- seq(0, max(mydata[mydata[,5] == 1, 1]) * 1.25, length.out=30);
    y1 <- seq(0, max(mydata[mydata[,5] == 1, 3]) * 1.25, length.out=30);
    data <- expand.grid(x1, y1);
    z <- seq(0, 1, length.out=900);
    for (k in 1:900) {z[k] <- Biv(mydata, data[k, 1], data[k, 2]);}
    z1 <- matrix(z, 30);
    if (col.biv == FALSE) {
      bw <- colours()[350-3*0:19];
      filled.contour(x1, y1, z1, xlab="Time in state 1",
        ylab="Time in state 2", col=bw, ...);
    } else {
      filled.contour(x1, y1, z1, xlab="Time in state 1",
        ylab="Time in state 2", ...);
    }
  }
  if (ntrans == 3 & pmar == TRUE) 
    cat("The plot for the marginal distribution of the second time cannot 
      be given for the illness-death model\n");
  if (ntrans == 3 & pbiv == TRUE) 
    cat("The plot for the bivariate distribution function cannot be given
      for the illness-death model\n");
} # plot.p3state