lines.CopulaCenR {CopulaCenR}

Plotting for CopulaCenR fits

Description

Plotting for CopulaCenR fits from ic_spTran_copula, rc_spCox_copula, ic_par_copula and rc_par_copula.

Usage

## S3 method for class 'CopulaCenR'
lines(
  x,
  y,
  class = "joint",
  newdata,
  evalPoints = 50,
  evalTimes1 = NULL,
  evalTimes2 = NULL,
  plot_margin = 1,
  cond_time = NULL,
  cond_margin = 2,
  plotly_object = NULL,
  ...
)

Arguments

x

an object of ic_spTran_copula, rc_spCox_copula, ic_par_copula, rc_par_copula

y

new data frame with colname names id, ind and covariate

class

one of "joint", "conditional" or "marginal"

newdata

new data frame (ignored if y is included)

evalPoints

number of time points to be evaluated; default is 50

evalTimes1

a vector of times for margin 1 to be evaluated; default is NULL; will override evalPoints if non-NULL

evalTimes2

a vector of times for margin 2 to be evaluated

plot_margin

for class = "marginal" only; indicator of which margin to plot (either 1 or 2); default is 1 for margin 1

cond_time

for class = "conditional" only; the time by which event has occurred in the margin indicated by cond_margin; must be smaller than the largest observed time

cond_margin

for class = "conditional" only; indicator of the margin where event has occurred (either 1 or 2); default is 2 for margin 2

plotly_object

only for class = "joint", an object of plot.CopulaCenR

...

further arguments

Details

y must be a data frame with columns id (subject id), ind (1,2 for two margins) and covariates.
The argument class determines the plot: "joint" for joint survival probabilities, "conditional" for conditional probabilities and "marginal" for marginal probabilities.

The function evaluates on a series of time points (given by evalPoints or evalTimes; evalTimes will override evalPoints). By default, the time points are automatically selected by specifying the number of points (evalPoints = 50). Users can also provide the specific time points through evalTimes1 and evalTimes2 for the two margins, respectively. When class = "conditional", only evalTimes1 is needed and the evaluation times are actually evalTimes1 plus cond_time.

If class = "conditional", one needs to specify the margin that has the event (by cond_margin) and time when the event has occurred (by cond_time). For example, if cond_margin = 2 and cond_time = 5, then the function produces the conditional survival probability (after time 5) in margin 1 given that margin 2 has got an event by time 5. This measurement is useful for predicting the second event given the first event has occurred. See the example for details.

If class = "marginal", one needs to specify which margin to plot through the argument plot_margin. See the example for details.

If class = "joint", one needs to include a plot_ly object (from plot.CopulaCenR with class = "joint") through the argument plotly_object. See the example for details.

Value

a 3D joint survival distribution plot if class = "joint"; a 2D survival distribution plot if class = "marginal" or "conditional".

Examples

data(AREDS)
# fit a Copula2-Sieve model
copula2_sp <- ic_spTran_copula(data = AREDS, copula = "Copula2",
              l = 0, u = 15, m = 3, r = 3,
              var_list = c("ENROLLAGE","rs2284665","SevScaleBL"))
newdata = data.frame(id = rep(1:3, each=2), ind = rep(c(1,2),3),
                     SevScaleBL = rep(3,6), ENROLLAGE = rep(60,6),
                     rs2284665 = c(0,0,1,1,2,2))
# Plot marginal survival probabilities
plot(x = copula2_sp, class = "marginal",
     newdata = newdata[newdata$id==1,],
     plot_margin = 1, ylim = c(0.6,1),
     ylab = "Marginal Survival Probability")
lines(x = copula2_sp, class = "marginal",
      newdata = newdata[newdata$id==2,],
      plot_margin = 1, lty = 2)
legend("bottomleft", c("id: 1","id: 2"), lty = c(1,2))

# Plot conditional survival probabilities
plot(x = copula2_sp, class = "conditional",
     newdata = newdata[newdata$id==1,],
     cond_margin = 2, cond_time = 5, ylim = c(0.25,1),
     xlab = "years", ylab = "Conditional Survival Probability")
lines(x = copula2_sp, class = "conditional",
      newdata = newdata[newdata$id==2,],
     cond_margin = 2, cond_time = 5, lty = 2)
legend("bottomleft", c("id: 1","id: 2"), lty = c(1,2))

# Plot joint survival probabilities
plot3d <- plot(x = copula2_sp, class = "joint",
               newdata = newdata[newdata$id==1,])
plot3d <- lines(x = copula2_sp, class = "joint",
                newdata = newdata[newdata$id==2,], plotly_object = plot3d)