rm(list = ls())
library(survival)
library(HQM) #used for loading the pbc2 data

"Biweight"<-function(x, ...) { ifelse(abs(x) < 1, (15/16) * (  1 - x^2)^2, ifelse(abs(x) > 1, 0, 0)) }
"Epanechnikov"<-    function(x, ... ) { ifelse(abs(x)<1,  0.75 * (1-x^2) , ifelse(abs(x)>1, 0,0))  }

"Gaussian"<- function(x, ...){ dnorm(x) }

"IntBiweight"<- function(x){ ifelse(x< -1, 0,ifelse(x> 1, 1, .5 + ((15*x)+ 3*(x^5))/16 - (5*x^3)/8)) }
"IntEpanechnikov"<-function(x) { ifelse(x< -2.236068, 0, ifelse(x> 2.236068, 1, .5- (2.236068*x*(x^2-15))/100)) }

"IntRectangular"<-  function(x){ ifelse(x< -1, 0, ifelse(x>1, 1, (x+1)/2))}
"IntTriangular"<- function(x){  ifelse(x< -1, 0, ifelse(x>1, 1, x+ ((x^2)*sign(x) - 1)/2)) }
"IntGaussian"<-  function(x){ pnorm(x) }

"Bandwidth"<-function(sample)
  #Rule of thumb bandwidth plug-in rules
  #input: sample data
  #output: a vector of two components:
  #Component 1 is the density bandwidth rule of Silverman (Density estimation for statistics and data analysis, Chapman and Hall, London, 1986), eq. 3.3.1
  #Component 2 is the distribution function bandwidth rule of Bowman, Halland Prvan (1988) (Bandwidth selection for the smoothing of distribution functions, Biometrika, Vol 85, page 801)
{
  n<-length(sample)
  StD<-sd(sample)
  IqR<-diff(quantile(sample, c(.25, .75)))
  tmp<-min(StD, IqR/1.34)
  c( 0.9*tmp*n^(-1/5) , 3.572* tmp * n^(-1/3)) #first entry is density simple plug in rule, second entry is distribution function simple plug in rule
}

"Epanechnikov"<-
  function(x){
    #Epanechikov kernel, square root of 5 is taken as 2.236068 (to save comp. time).
    ifelse(abs(x)< 2.236068, (3/4)*(1-((x^2)/5 ))/2.236068, ifelse(abs(x)>2.236068, 0,0))
  }


"IntEpanechnikov"<-
  function(x){ #Integrated Epanechikov kernel, square root of 5 is taken as 2.236068 (to save comp. time).
    ifelse(x< -2.236068, 0, ifelse(x> 2.236068, 1, .5- (2.236068*x*(x^2-15))/100))
  }


"HazardRateEst"<- function(xin, xout, kfun, h, ci)
{
  n<- length(xin)
  n1<-1:n
  xin.use<-sort(xin)
  ci<-ci[order(xin)]
  arg1<-(sapply(xout, "-", xin.use))/h
  arg2<-kfun(arg1)* ci/(n-n1+1)
  arg3<-colSums(arg2)/h
  arg3
}


"HStar"<-function(xin, cens, xout)
  #Kaplan-Meier Estimate, slightly modified so that it avoids division by zero.
  #In: sample and censoring indicators, i.e. xin and cens correspond to $(Z_i, \Lambda_i)$ in manuscript
  #(in fact the first two lines of the routine ensure that this will be the case anyway). The xout argument denotes the points where the estimate will be evaluated
  #the most common option is to use xin in the place of xout
  #Out: 1-\hat H^\star(X_i) (as given on the manuscript)
{
  ordxin <-  sort(xin) #order input data=Xi
  ordcens <- cens[order(xin)] #order censoring =Di
  n <- length(xin)
  nn <- length(xout)
  vKuse <- sapply(1:nn, function (i, ordxin, xout) max(c(max(which(ordxin <= xout[i]))
                                                         , min(which(xout[i] <= ordxin)))), ordxin, xout)
  vK <- vKuse - 1
  Huse<-sapply(vK, function (i,n,ordcens) ((n-1:i+1)/(n-1:i+2))^(1-ordcens[1:i]), n, ordcens)
  result<-sapply(Huse, prod)
  result[order(xin)]<-result
  result
}



"ikde"<-function(xin, xout, h, kfun, ci)
{ ###SORT VERSION!!!!
  n<- length(xin)
  ci <- ci[order(xin)]
  xin<-sort(xin)
  arg1<-(sapply(xout, "-", xin))/h

  test21<- survfit(Surv(xin, ci) ~ 1)
  denom.km<- summary(test21, times = xin)$surv

  arg2<-kfun(arg1) * ci/denom.km # HStar(xin,  ci, xin) #(fit1$surv+0.01) #
  arg3<-colSums(arg2)
  arg3/n #
}

"HazardRatioEst"<- function(xin, xin0, xout0, kfun, ikfun, h,   ci, ci0)
{
  n<- length(xin0)
  b0<-   Bandwidth(xin0)[2] # NEED TO REPLACE WITH OPTIMAL RULE
  hatF0x<- ikde(xin0, xout0, b0, ikfun, ci0)
  hatF0xi<-  ikde(xin0 , sort(xin), b0, ikfun, ci0)  # ci0 or ci?


  arg1<-   (sapply(hatF0x, "-", hatF0xi))/h

  test2<- survfit(Surv(xin, ci) ~ 1)
  hatG1xi.km<- summary(test2, times = xin)$surv

  test2<- survfit(Surv(xin0, ci0) ~ 1)
  hatF1xi.km<- summary(test2, times = xin)$surv

  arg2<-kfun(arg1)*  ci  *  (( 1-hatF0xi) /  (hatG1xi.km *hatF1xi.km ) )  #((1- hatF1xi )*(1 - hatG1xi))
  arg3<-colSums(arg2)/( h*n)
  arg3
}


 
pbc2.use<-pbc2

pbc2.use2<- matrix(0, nrow=312, ncol=20)
pbc2.use2 <- as.data.frame(pbc2.use2)
 colnames(pbc2.use2)<-names(pbc2)
 for(i in 1:312)
 {
   test1<-pbc2[which(pbc2[, 1]==i), ]
   pbc2.use2[i,]<-test1[which.max(test1[, 7]), ]
 }

pbc2.placebo<-pbc2.use2[which(pbc2.use2[,4]==1),]
time.placebo<-pbc2.placebo[,2]


pbc2.Dpenic<-pbc2.use2[which(pbc2.use2[,4]==2), ]
time.Dpenic<-pbc2.Dpenic[,2]

xout1<-time.placebo
xout2<-time.Dpenic
xout<-max(xout1, xout2)
xout<-seq(0.1, xout, length=50)
cen = pbc2.use2[which(pbc2.use2[,4]==1), ]$status2

cen0 = pbc2.use2[which(pbc2.use2[,4]==2), ]$status2

 

 
xout2<-seq(0.1, 12, length=50)

arg2.placebo<-HazardRateEst(time.placebo, xout2, Epanechnikov, 1.2, cen)
arg2.Dpenic<-HazardRateEst(time.Dpenic, xout2, Epanechnikov, 1.2, cen0)

 
plot(xout2, arg2.placebo, type="l", xlab="time (years)", ylim=c(0, 0.3), ylab="Baseline hazards: Placebo vs D-penicil", col=4, lwd=2)
lines(xout2, arg2.Dpenic, col=2, lwd=2)


