if (FALSE) {
## 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 warning), 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
## bandwidth 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