# NOT RUN {
# }
# NOT RUN {
## Analysis with the example dataset (pscheck)
## Cumulative hazard function (chf) estimation
# Presmoothed estimate of chf with bootstrap bandwidth (fixing the randomization seed
# makes further comparisons easier)
set.seed(1)
Hboot1 <- presmooth(t, delta, pscheck, estimand = "H", bw.selec = "bootstrap")
# As above, but: 1) specifying the points where the estimate is computed (note the war-
# ning), 2) specifying the search grid for the bandwidth, and 3) saving the bootstrap
# MISE
set.seed(1)
Hboot2 <- presmooth(t, delta, pscheck, estimand = "H", bw.selec = "bootstrap", x.est =
seq(0, 1, by = 0.02), grid.bw = seq(0.55, 0.7, by = 0.01), control =
control.presmooth(save.mise = TRUE))
# A plot of the MISE, showing the bootstrap bandwidth
with(Hboot2,{
plot(grid.bw$grid.bw, mise, xlab = "Bandwidth", ylab = "MISE", main =
expression(paste("Bootstrap bandwidth, ", b[boot])), type = "l")
points(bandwidth, mise[grid.bw$grid.bw == bandwidth], pch = 46, cex = 5)
segments(bandwidth, 0, bandwidth, mise[grid.bw$grid.bw == bandwidth], lty = "dotted")
text(bandwidth, min(mise), bquote(paste(" ", b[boot] == .(bandwidth))), adj = c(0, 0))
})
# A plot of the presmoothed chf compared with Nelson-Aalen estimate and the true curve.
# The point (0, 0) must be added.
plot(c(0, Hboot2$x.est), c(0, Hboot2$estimate), xlab = "t", ylab = "Cumulative hazard",
type = "s")
Hna <- presmooth(t, delta, pscheck, estimand = "H", bw.selec = "fixed", fixed.bw = 0)
lines(c(0, Hna$x.est), c(0, Hna$estimate), type = "s", col = "red")
curve(x^3, add = TRUE, col = "grey", lty = "dotted")
legend("topleft", legend = c("Presmoothed Nelson-Aalen", "Nelson-Aalen", "True"), col =
c("black", "red", "grey"), lty = c("solid", "solid", "dotted"))
# An alternative way of obtaining the Nelson-Aalen estimate
Hna <- presmooth(t, delta, pscheck, "H", presmoothing = FALSE)
## Hazard function (hf) estimation
# (An example where right-boundary correction is successful) Presmoothed estimate of hf
# 1) with plug-in bandwidth with and without right-boundary correction, 2) specifying the
# grid for presmoothing bandwidth selection and 3) specifying the support of the weight
# function
hpi1 <- presmooth(t, delta, pscheck, estimand = "h", bw.selec = "plug-in", x.est =
seq(0, max(pscheck$t), by = 0.02), grid.bw.pil = seq(range(pscheck$t)[1],
range(pscheck)[2], by = 0.01), control = control.presmooth(q.weight = c(0.25, 0.75)))
hpi2 <- presmooth(t, delta, pscheck, estimand = "h", bw.selec = "plug-in", bound = "right",
x.est = seq(0, max(pscheck$t), by = 0.02), grid.bw.pil = seq(range(pscheck$t)[1],
range(pscheck$t)[2], by = 0.01), control = control.presmooth(q.weight = c(0.25, 0.75)))
plot(hpi1$x.est, hpi1$estimate, xlab = "t", ylab = "Hazard", ylim = c(0,
max(pmax(hpi1$estimate, hpi2$estimate))), type = "l")
lines(hpi2$x.est, hpi2$estimate, col = 2)
legend("bottomright", legend = c("none", "right"), title = "Boundary effect correction",
col = 1:2, lty = 1)
# Estimation of hf using a bootstrap bandwidth. The values chosen for the grid.bw argument
# are based on the result of preliminary trials (Warning: it may take a while ...)
set.seed(1)
hboot <- presmooth(t, delta, pscheck, estimand = "h", bw.selec = "bootstrap", bound = "right",
x.est = seq(0, max(pscheck$t), by = 0.02), grid.bw.pil = seq(range(pscheck$t)[1],
range(pscheck)[2], by = 0.01), grid.bw = list(seq(0.4, 0.8, by = 0.05),
seq(0.6, 0.9, by = 0.005)), control = control.presmooth(q.weight = c(0.25, 0.75),
save.mise = TRUE))
# A plot of the bootstrap MISE, showing the bootstrap bandwidth
with(hboot,{
contour(grid.bw$grid.bw.1, grid.bw$grid.bw.2, mise, levels = seq(min(mise), max(mise),
length.out = 20), xlab = "Presmoothing bandwidth", ylab = "Smoothing bandwidth",
main = "Bootstrap MISE")
points(bandwidth[1], bandwidth[2], pch = 46, cex = 6)
text(bandwidth[1], bandwidth[2], bquote(paste(" ", b[boot], symbol("="), symbol("("),
.(bandwidth[1]), symbol(","), .(bandwidth[2]), symbol(")"))), adj = c(1, 0))
})
# Compare with the hf estimate obtained with plug-in bandwidth and the true curve
plot(hpi2$x.est, hpi2$estimate, xlab = "t", ylab = "Hazard", ylim = c(0,
max(pmax(hpi2$estimate, hboot$estimate))), type = "l")
lines(hboot$x.est, hboot$estimate, col = 2)
curve(3*x^2, add = TRUE, col = "grey", lty = "dotted")
legend("bottomright", legend = c("Plug-in bandwidth", "Bootstrap bandwidth", "True"),
col = c("black", "red", "grey"), lty = c("solid", "solid", "dotted"))
## Density function (df) estimation
# Presmoothed estimate of df with plug-in and bootstrap bandwidths (with default options)
# and comparison with the true df
dpi <- presmooth(t, delta, pscheck, estimand = "f", bw.selec = "plug-in")
# The bootstrap presmoothing bandwidth is on the right boundary of the default grid (which
# in fact is the upper bound for the bandwidth: the range of the observed times)
set.seed(1)
dboot <- presmooth(t, delta, pscheck, estimand = "f", bw.selec = "bootstrap")
# For this example, the estimates with either plugin or bootstrap banwidth are very similar
plot(dpi$x.est, dpi$estimate, xlab = "t", ylab = "Density", ylim = c(0,
max(pmax(dpi$estimate, dboot$estimate))), type = "l", col = "blue", lty = 2)
lines(dboot$x.est, dboot$estimate, col = "red", lty = 4)
curve(3*x^2*exp(-x^3), add = TRUE, lty = 1)
legend("bottomright", legend = c("Plug-in bandwidth", "Bootstrap bandwidth", "True"),
col = c("blue", "red", "black"), lty = c(2, 4, 1))
# }
Run the code above in your browser using DataLab