Summary of Splus Estimation Steps for Kernel-Smoothed 
Hazard Estimator, illustrated with Bone Marrow 
Transplant Data from Sec. 1.3.
============================================== 10/24/03

Bone Marrow Transplant data already imported as data-frame 
"BMT" (in problem (B) of HWset2 -- see HW2 solutions) with 
variables  

> names(BMT)
[1] "Gp"    "Dtime" "Rtime" "Dind"  "Rind" 

and with patients in 3 groups 

> table(BMT$Gp)
  1  2  3          1=ALL, 2=AML lo-risk, 3 = High-risk
 38 54 45

### We are interested in disease-free survival time, which is
###   Rtime , and consider the range of observed times
> summary(BMT$Rtime[(BMT$Rind==1 | BMT$Dind==1) & BMT$Gp==1])
 Min. 1st Qu. Median  Mean 3rd Qu. Max. 
    1   108.5    182 248.6   391.8  662    ### 24 values

### We use Epanechnikov kernel (6.2.2) with modification
###   (6.2.7) at small and large times, as given in the Klein and
###   Moeschberger book. Note that in Splus we want to code the
##    kernel in such a way that the t-input can be a vector.

> abfcn <- function(q) {
       aux <- (1+q)^4*(19-18*q+3*q^2)
       matrix(c(64*(2-4*q+6*q^2-3*q^3),240*(1-q)^2)/aux, ncol=2) }

> abfcn(0.5)
         [,1]     [,2] 
[1,] 1.322997 1.102498
> abfcn(0.62)
         [,1]      [,2] 
[1,] 1.148371 0.5595052


> Kernmat <- function(t, tv, b, tD) { 
### Modified Epanechnkov kernel function, made into matrix
    kmat <- pmax(1-(outer(t,tv,"-")/b)^2,0)*(0.75/b)
    loind <- t <= b
    hiind  <- t >= tD-b
    qL <- t[loind]/b
    qH <- (tD-t[hiind])/b 
    ablo <- abfcn(qL)
    abhi <- abfcn(qH)
    if(sum(loind)) kmat[loind,] <- kmat[loind,] * (outer(ablo[,1],
        rep(1,length(tv)),"*") + outer(ablo[,2]/b, rep(1,length(tv)),
        "*")*outer(t[loind],tv,"-"))
    if(sum(hiind)) kmat[hiind,] <- kmat[hiind,] * (outer(abhi[,1], 
        rep(1,length(tv)),"*") - outer(abhi[,2]/b, rep(1,length(tv)),
        "*")*outer(t[hiind],tv,"-"))
    kmat }

### Now the kernel estimator at t, if tD is selected as
###   rightmost interval point and  b  as bandwidth, is
###   c(Kernmat(t,tv,b,tD) %*% hv)  if  tv  and  
###   respectively  hv  are the obs-death points and 
###   (Nelson-Aalen) hazard increments. 

> BMTfit <- survfit(Surv(BMT$Rtime[BMT$Gp==1], 
       pmax(BMT$Rind,BMT$Dind)[BMT$Gp==1]))
### 23 distinct death-times in this fitted object
> tv <- BMTfit$time[BMTfit$n.ev > 0]
  hv <- BMTfit$n.ev[BMTfit$n.ev > 0]/BMTfit$n.risk[BMTfit$n.ev > 0]

> round(c(Kernmat(c(150,50,600),tv, 100, 662 ) %*% hv),6)
[1] 0.002570 0.001515 0.001324       ### Agrees with book, pp.168-9

> round(Kernmat(c(150,50,600),tv, 100, 662 ),6)
         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8] 
[1,] 0.000000 0.000731 0.003168 0.004428 0.005913 0.006113 0.006239 0.006300
[2,] 0.010619 0.009485 0.007481 0.006046 0.003866 0.003517 0.003288 0.003175
[3,] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
         [,9]    [,10]    [,11]    [,12]    [,13]  [,14] [,15] [,16] [,17] 
[1,] 0.006912 0.007169 0.007137 0.006177 0.006048 0.0027     0     0     0
[2,] 0.001911 0.001274 0.000000 0.000000 0.000000 0.0000     0     0     0
[3,] 0.000000 0.000000 0.000000 0.000000 0.000000 0.0000     0     0     0
     [,18] [,19] [,20]    [,21]    [,22]    [,23] 
[1,]     0     0     0 0.000000 0.000000 0.000000
[2,]     0     0     0 0.000000 0.000000 0.000000
[3,]     0     0     0 0.002492 0.008918 0.006904