[R] numerical approximation of survival function

[Carlo Giovanni Camarda]

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]]