Carlo Giovanni Camarda
2016-Nov-18 16:30 UTC
[R] numerical approximation of survival function
Dear R-users, my question concerns numerical approximation and, somehow, survival analysis. Let?s assume we have a density function and we aim to numerically compute the hazard, which is, in theory, the ratio between density and survival. In the following example, I take a pdf from a log-normal and proceed as I would have known only this function. I numerically approximate the survival function by splinefun() and integrate().Finally I compute the hazard. It?s easy to see that (1) true hazard is different from the numerically approximated one and (2) discrepancy depends upon the approximation carried out to compute the survival. Is there any way to overcome or attenuate this issue? Or do I just miss something obvious? In my real study I do not deal with analytically defined pdf, so I would like to be sure that approximation is done at its best. Thanks in advance for your help. Giancarlo Camarda ## x-values delta <- 0.1 x <- seq(0, 200, by=delta) m <- length(x) ## true pdf fxT <- dlnorm(x, meanlog=3.5, sdlog=0.5) fxT <- fxT/sum(fxT) ## true survival SxT <- 1 - plnorm(x, meanlog=3.5, sdlog=0.5) ## true hazard hxT <- fxT/SxT ## numerical approximation of the pdf fx <- spline(x, fxT, n=m)$y ## numerical approximation of the survival Sx <- numeric(m) fxfun <- splinefun(x, fx/delta) for(j in 1:m){ Sx[j] <- 1 - integrate(fxfun, x[1], x[j])$value } ## hazard: pdf divided by survival hx <- fx/Sx ## both numerical approx hx1 <- fx/SxT ## only the survival numerically approx hx2 <- fxT/Sx ## only the pdf numerically approx ## plotting hazards plot(x, hxT) lines(x, hx, col=2, lwd=2) lines(x, hx1, col=3, lwd=2) lines(x, hx2, col=4, lty=2, lwd=2) ## it seems that we have got something ## at the survival approx level: plot(x, Sx/SxT) [[alternative HTML version deleted]]