TransHazRateEst {NPHazardRate} | R Documentation |
Implements the transformated kernel hazard rate estimator of Bagkavos (2008). The estimate is expected to have less bias compared to the ordinary kernel estimate HazardRateEst
. The estimate results by first transforming the data to a sample from the exponential distribution through the integrated hazard rate function, estimated by iHazardRateEst
and uses the result as input to the classical kernel hazard rate estimate HazardRateEst
. An inverse transform turn the estimate to a hazard rate estimate of the original sample. See section "Details" below.
TransHazRateEst(xin, xout, kfun, ikfun, h1, h2, ci)
xin |
A vector of data points. Missing values not allowed. |
xout |
A vector of points at which the hazard rate function will be estimated. |
kfun |
Kernel function to use. Supported kernels: Epanechnikov, Biweight, Gaussian, Rectangular, Triangular, HigherOrder. |
ikfun |
An integrated kernel function to use. Supported kernels: Epanechnikov, Biweight, Gaussian, Rectangular, Triangular, HigherOrder. |
h1 |
A scalar, pilot bandwidth. |
h2 |
A scalar, transformed kernel bandwidth. |
ci |
A vector of censoring indicators: 1's indicate uncensored observations, 0's correspond to censored obs. |
The transformed kernel hazard rate estimate of Bagkavos (2008) is given by
\hat \lambda_t(x;h_1, h_2) = \sum_{i=1}^n \frac{K_{h_2}\left \{ (\hat \Lambda(x;h_1 ) - \hat \Lambda(X_{(i)};h_1 ) ) \right \}\delta_{(i)}}{n-i+1}\hat \lambda(x;h_1 ).
The estimate uses the classical kernel hazard rate estimate \lambda(x; h_1)
implemented in HazardRateEst
and its integrated version
\hat \Lambda(x; h_1) = \sum_{i=1}^n \frac{k\left \{(x-X_{(i)})h_1^{-1}\right \}\delta_{(i)}}{n-i+1}
where
k(x) = \int_{-\infty}^x K(y)\,dy
implemented in iHazardRateEst
. The pilot bandwidth h_1
is determined by an optimal bandwidth rule such as PlugInBand
.
TO DO: Insert a rule for the adaptive bandwidth h_2
.
A vector with the values of the function at the designated points xout.
Bagkavos (2008), Transformations in hazard rate estimation, J. Nonparam. Statist., 20, 721-738
VarBandHazEst, HazardRateEst, PlugInBand
x<-seq(0, 5,length=100) #design points where the estimate will be calculated
plot(x, HazardRate(x, "weibull", .6, 1), type="l",
xlab = "x", ylab="Hazard rate") #plot true hazard rate function
SampleSize <- 100
mat<-matrix(nrow=SampleSize, ncol=20)
for(i in 1:20)
{ #Calculate the average of 20 estimates and draw on the screen
ti<- rweibull(SampleSize, .6, 1) #draw a random sample from the actual distribution
ui<-rexp(SampleSize, .05) #draw a random sample from the censoring distribution
cat("\n AMOUNT OF CENSORING: ", length(which(ti>ui))/length(ti)*100, "\n")
x1<-pmin(ti,ui) #this is the observed sample
cen<-rep.int(1, SampleSize) #censoring indicators
cen[which(ti>ui)]<-0 #censored values correspond to zero
h2<-DefVarBandRule(ti, cen) #Deafult Band. Rule - Weibull Reference
huse1<- PlugInBand(x1, x, cen, Biweight) #
mat[,i]<-TransHazRateEst(x1,x,Epanechnikov,IntEpanechnikov,huse1,h2,cen)
}
lines(x, rowMeans(mat) , lty=2) #draw the average transformed estimate