# NOT RUN {
## Some artificial data
set.seed(123)
n <- 50
x <- runif(n, -2, 2) ## Covariate values
y <- rweibull(n, shape = .5*(x + 4)) ## True lifetimes
c <- rexp(n) ## Censoring values
p <- exp(2*x)/(1 + exp(2*x)) ## Probability of being susceptible
u <- runif(n)
t <- ifelse(u < p, pmin(y, c), c) ## Observed times
d <- ifelse(u < p, ifelse(y < c, 1, 0), 0) ## Uncensoring indicator
data <- data.frame(x = x, t = t, d = d)
## Latency estimates for covariate value 0.5...
x0 <- .5
## ... (a) with global bandwidths 0.5, 1, 2.
## By default, estimates are computed at the time values of 't'
S1 <- latency(x, t, d, data, x0 = x0, h = c(.5, 1, 2), local = FALSE)
plot(S1$testim, S1$S$h0.5$x0.5, type = "s", xlab = "Time", ylab =
"Latency", ylim = c(0, 1))
lines(S1$testim, S1$S$h1$x0.5, type = "s", lty = 2)
lines(S1$testim, S1$S$h2$x0.5, type = "s", lty = 3)
## The true latency curve is plotted as reference
lines(S1$testim, pweibull(S1$testim, shape = .5*(x0 + 4), lower.tail =
FALSE), col = 2)
legend("topright", c(paste("Estimate, ", c("h = 0.5", "h = 1",
"h = 2")), "True"), lty = c(1:3, 1), col = c(rep(1, 3), 2))
## As before, but with estimates computed at times 0.1, 0.2,..., 1
S2 <- latency(x, t, d, data, x0 = x0, h = c(.5, 1, 2), local = FALSE,
testimate = .1*(1:10))
## ... (b) with local bandwidth 2.
S3 <- latency(x, t, d, data, x0 = x0, h = 2, local = TRUE)
#### Note that with only one covariate value the results with
#### 'local = FALSE' and 'local = TRUE' coincide, but the output formats
#### differ slightly. Compare with
S3 <- latency(x, t, d, data, x0 = x0, h = 2, local = FALSE)
## ... (c) with local bootstrap bandwidth
b <- latencyhboot(x, t, d, data, x0 = x0)
S4 <- latency(x, t, d, data, x0 = x0, h = b$h)
## ... (d) when the bandwidth is not specified, the bootstrap bandwidth
#### selector given by the 'latencyhboot' function is used by default.
#### The computation of 95% confidence intervals based on 1999 bootstrap
#### resamples is also illustrated
S5 <- latency(x, t, d, data, x0 = x0, conflevel = .95, bootpars =
controlpars(B = 1999))
plot(S5$testim, S5$S$x0, type = "s", xlab = "Time", ylab = "Latency",
ylim = c(0, 1))
lines(S5$testim, S5$conf$x0$lower, type = "s", lty = 2)
lines(S5$testim, S5$conf$x0$upper, type = "s", lty = 2)
lines(S5$testim, pweibull(S5$testim, shape = .5*(x0 + 4), lower.tail =
FALSE), col = 2)
legend("topright", c("Estimate", "95% CI limits", "True"), lty = c(1,
2, 1), col = c(1, 1, 2))
# }
# NOT RUN {
## Example with the dataset 'bmt' of the 'KMsurv' package
## to study the survival of the uncured patients aged 25 and 40
data("bmt", package = "KMsurv")
x0 <- c(25, 40)
S <- latency(z1, t2, d3, bmt, x0 = x0, conflevel = .95)
## Plot of predicted latency curves and confidence intervals
plot(S$testim, S$S$x25, type = "s", xlab = "Time (days)",
ylab = "Latency", ylim = c(0,1))
lines(S$testim, S$conf$x25$lower, type = "s", lty = 2)
lines(S$testim, S$conf$x25$upper, type = "s", lty = 2)
lines(S$testim, S$S$x40, type = "s", lty = 1, col = 2)
lines(S$testim, S$conf$x40$lower, type = "s", lty = 2, col = 2)
lines(S$testim, S$conf$x40$upper, type = "s", lty = 2, col = 2)
legend("topright", c("Age 25: Estimate", "Age 25: 95% CI limits",
"Age 40: Estimate","Age 40: 95% CI limits"), lty = 1:2,
col = c(1, 1, 2, 2))
# }
Run the code above in your browser using DataLab