Nothing
R/Uhaz.R
In npsurv: Nonparametric Survival Analysis
Defines functions
uh.initial
plotgraduh
plotsurvuh
plotdenuh
plotchazuh
plothazuh
plot.uh
plot.Uhaz
denuh
survuh
chazuh
hazuh
print.uh
uh
collapse
simplify
maxgrad
dlogLik
logLikuh
grad2
grad1
grad0
backtrack
updatemass
Uhaz
Documented in
chazuh
denuh
hazuh
logLikuh
plotchazuh
plotdenuh
plotgraduh
plothazuh
plotsurvuh
plot.uh
plot.Uhaz
print.uh
survuh
uh
Uhaz
############################################
# Estimation of a U-shaped Hazard Function #
############################################
##'U-shaped Hazard Function Estimation
##'
##'
##'\code{Uhaz} computes the nonparametric maximum likelihood esimate (NPMLE) of
##'a U-shaped hazard function from exact or interval-censored data, or a mix of
##'the two types of data.
##'
##'
##'If \code{data} is a vector, it contains only exact observations, with
##'weights given in \code{w}.
##'
##'If \code{data} is a matrix with two columns, it contains interval-censored
##'observations, with the two columns storing their left and right end-points,
##'respectively. If the left and right end-points are equal, then the
##'observation is exact. Weights are provided by \code{w}.
##'
##'If \code{data} is a matrix with three columns, it contains interval-censored
##'observations, with the first two columns storing their left and right
##'end-points, respectively. The weight of each observation is the third-column
##'value multiplied by the corresponding weight value in \code{w}.
##'
##'The algorithm used for the computing the NPMLE of a hazard function under
##'the U-shape restriction is is proposed by Wang and Fani (2015). Such a
##'hazard function is given by
##'
##'A U-shaped hazard function is given by
##'
##'\deqn{h(t) = \alpha + \sum_{j = 1}^k \nu_j(\tau_j - t)_+^p + \sum_{j = 1}^{m} \mu_j (t-\eta_j)_+^p,}{ h(t) = alpha + sum_{j=1}^k nu_j (tau_j - t)_+^p
##' + sum_{j=1}^m mu_j (t - eta_j)_+^p,}
##'
##'where \eqn{\alpha,\nu_j,\mu_j \ge 0}{alpha, nu_j, mu_j \ge 0},
##'\eqn{\tau_1 < \cdots < \tau_k \le \eta_1 < \cdots < \eta_m,}{tau_1 < ... < tau_k <= eta_1 < ... < eta_m,} and \eqn{p \ge 0}{p >= 0} is the the spline degree which
##'determines the smoothness of the U-shaped hazard. As p increases, the family
##'of hazard functions becomes increasingly smoother, but at the time, smaller.
##'When p = 0, the hazard function is U-shaped, as studied by Bray et al.
##'(1967). When p = 1, the hazard function is convex, as studied by Jankowski
##'and Wellner (2009a,b).
##'
##'Note that \code{deg} (i.e., p in the above mathematical display) can take on
##'any nonnegative real value.
##'
##'@aliases Uhaz Uhaz.object
##'@param data vector or matrix, or an object of class \code{icendata}.
##'@param w weights or multiplicities of the observations.
##'@param deg nonnegative real number for spline degree (i.e., p in the formula
##'below).
##'@param maxit maximum number of iterations.
##'@param tol tolerance level for stopping the algorithm. It is used as the
##'threshold on the increase of the log-likelihood after each iteration.
##'@param verb verbosity level for printing intermediate results in each
##'iteration.
##'@return
##'
##'An object of class \code{Uhaz}, which is a list with components:
##'
##'\item{convergence}{= \code{TRUE}, converged successfully;
##'
##'= \code{FALSE}, maximum number of iterations reached.}
##'
##'\item{grad}{gradient values at the knots.}
##'
##'\item{numiter}{number of iterations used.}
##'
##'\item{ll}{log-likelihood value of the NPMLE \code{h}.}
##'
##'\item{h}{NPMLE of the U-shaped hazard function, an object of class
##'\code{uh}.}
##'@author Yong Wang <yongwang@@auckland.ac.nz>
##'@seealso \code{\link{icendata}}, \code{\link{nzmort}}.
##'@references
##'
##'Bray, T. A., Crawford, G. B., and Proschan, F. (1967). \emph{Maximum
##'Likelihood Estimation of a U-shaped Failure Rate Function}. Defense
##'Technical Information Center.
##'
##'Jankowski, H. K. and Wellner, J. A. (2009a). Computation of nonparametric
##'convex hazard estimators via profile methods. \emph{Journal of Nonparametric
##'Statistics}, \bold{21}, 505-518.
##'
##'Jankowski, H. K. and Wellner, J. A. (2009b). Nonparametric estimation of a
##'convex bathtub-shaped hazard function. \emph{Bernoulli}, \bold{15},
##'1010-1035.
##'
##'Wang, Y. and Fani, S. (2018). Nonparametric maximum likelihood
##' computation of a U-shaped hazard function. \emph{Statistics and
##' Computing}, \bold{28}, 187-200.
##'
##'@keywords function
##'@examples
##'
##'## Interval-censored observations
##'data(ap)
##'(r = Uhaz(ap, deg=0))
##'plot(r, ylim=c(0,.3), col=1)
##'for(i in 1:6) plot(Uhaz(ap, deg=i/2), add=TRUE, col=i+1)
##'legend(15, 0.01, paste0("deg = ", 0:6/2), lwd=2, col=1:7, xjust=1, yjust=0)
##'
##'## Exact observations
##'data(nzmort)
##'x = with(nzmort, nzmort[ethnic=="maori",])[,1:2] # Maori mortality
##'(h0 = Uhaz(x[,1]+0.5, x[,2], deg=0)$h) # U-shaped hazard
##'(h1 = Uhaz(x[,1]+0.5, x[,2], deg=1)$h) # convex hazard
##'(h2 <- Uhaz(x[,1]+0.5, x[,2], deg=2)$h) # smooth U-shaped hazard
##'
##'plot(h0, pch=2) # plot hazard functions
##'plot(h1, add=TRUE, col="green3", pch=1)
##'plot(h2, add=TRUE, col="red3", pch=19)
##'
##'age = 0:max(x[,1]) # plot densities
##'count = integer(length(age))
##'count[x[,"age"]+1] = x[,"deaths"]
##'barplot(count/sum(count), space=0, col="lightgrey")
##'axis(1, pos=NA, at=0:10*10)
##'plot(h0, fn="d", add=TRUE, pch=2)
##'plot(h1, fn="d", add=TRUE, col="green3", pch=1)
##'plot(h2, fn="d", add=TRUE, col="red3", pch=19)
##'
##'plot(h0, fn="s", pch=2) # plot survival functions
##'plot(h1, fn="s", add=TRUE, col="green3", pch=1)
##'plot(h2, fn="s", add=TRUE, col="red3", pch=19)
##'
##'## Exact and right-censored observations
##'data(gastric)
##'plot(h0<-Uhaz(gastric, deg=0)$h) # plot hazard functions
##'plot(h1<-Uhaz(gastric, deg=1)$h, add=TRUE, col="green3")
##'plot(h2<-Uhaz(gastric, deg=2)$h, add=TRUE, col="red3")
##'
##'plot(npsurv(gastric), fn="s", col="grey") # plot survival functions
##'plot(h0, fn="s", add=TRUE)
##'plot(h1, fn="s", add=TRUE, col="green3")
##'plot(h2, fn="s", add=TRUE, col="red3")
##'
##'@export Uhaz
Uhaz = function(data, w=1, deg=1, maxit=100, tol=1e-6, verb=0) {
x = icendata(data, w)
h = uh.initial(x, deg)
attr(h, "ll") = logLikuh(h, x)
expdH = NULL
bc = TRUE # boundary change
convergence = 1
for(i in 1:maxit){
h.old = h
if(nrow(x$o) > 0) expdH = exp(chazuh(x$o[,1],h) - chazuh(x$o[,2],h))
maxima = maxgrad(h, x, expdH, bc=bc)
np1 = maxima$np1
np2 = maxima$np2
h = uh(h$alpha, c(h$tau, np1), c(h$nu, double(length(np1))),
c(h$eta, np2), c(h$mu, double(length(np2))),
h$upper, h$deg, collapse=TRUE)
r = updatemass(h, x, expdH, tol=tol)
h = r$h
if(h$deg == 0) {h = simplify(h); attr(h, "ll") = logLikuh(h, x)}
if(verb>0) { cat("##### Iteration", i, "#####\n")
cat("Log-likelihood: ", signif(attr(h,"ll"), 6), "\n")
if(verb>1) cat("Gradient values: ", signif(dlogLik(h, x), 6), "\n")
if(verb>2) {cat("hazard function:\n"); print(h)} }
if(r$convergence == 1) bc = FALSE # backtracking failed.
else if(attr(h, "ll") <= attr(h.old, "ll") + tol) {convergence = 0; break}
}
r = list(convergence=convergence, grad=dlogLik(h, x), numiter=i,
ll=attr(h, "ll"), h=h)
class(r) = "Uhaz"
r
}
# Update masses
updatemass = function(h, x, expdH=NULL, tol=1e-10) {
tau = h$tau
k = length(tau)
j2 = h$eta != h$upper
eta = h$eta = h$eta[j2]
h$mu = h$mu[j2]
m = length(eta)
p = h$deg
D1 = D2 = NULL
t1 = x$t[x$i1]
n1 = length(t1)
if(n1 > 0) {
tau.r = rep(tau, rep.int(n1,k))
dim(tau.r) = c(n1, k)
if(p > 0) tau.t = pmax(tau.r - t1, 0)
if(m > 0) {
eta.r = rep(eta, rep.int(n1,m))
dim(eta.r) = c(n1, m)
if(p > 0) t.eta = pmax(t1 - eta.r, 0)
}
D1 = switch(as.character(p),
"0" = cbind(1, tau.r >= t1, if(m>0) t1 >= eta.r else NULL) / hazuh(t1, h),
"1" = cbind(1, tau.t, if(m>0) t.eta else NULL) / hazuh(t1, h),
"2" = cbind(1, tau.t * tau.t,
if(m>0) t.eta * t.eta else NULL) / hazuh(t1, h),
cbind(1, tau.t^p, if(m>0) t.eta^p else NULL) / hazuh(t1, h) )
}
n2 = nrow(x$o)
if(n2 > 0) {
if(is.null(expdH)) expdH = exp(chazuh(x$o[,1],h) - chazuh(x$o[,2],h))
delta = sqrt(expdH) / (1 - expdH)
tau.r1 = rep(tau, rep.int(n2,k))
dim(tau.r1) = c(n2, k)
tau.x1 = pmax(tau.r1 - x$o[,1], 0)
tau.x2 = pmax(tau.r1 - x$o[,2], 0)
xd0 = x$o[,1] - x$o[,2]
xd1 = switch(as.character(p),
"0" = tau.x2 - tau.x1,
"1" = .5 * (tau.x2 * tau.x2 - tau.x1 * tau.x1),
"2" = (tau.x2 * tau.x2 * tau.x2 - tau.x1 * tau.x1 * tau.x1) / 3,
(tau.x2^(p+1) - tau.x1^(p+1)) / (p+1) )
if(m > 0) {
eta.r2 = rep(eta, rep.int(n2,m))
dim(eta.r2) = c(n2, m)
x1.eta = pmax(x$o[,1] - eta.r2, 0)
x2.eta = pmax(x$o[,2] - eta.r2, 0)
xd2 = switch(as.character(p),
"0" = x1.eta - x2.eta,
"1" = .5 * (x1.eta * x1.eta - x2.eta * x2.eta),
"2" = (x1.eta * x1.eta * x1.eta - x2.eta * x2.eta * x2.eta) / 3,
(x1.eta^(p+1) - x2.eta^(p+1)) / (p+1) )
}
else xd2 = NULL
D2 = cbind(xd0, xd1, xd2) * delta
D2[delta == 0] = 0
}
D = rbind(D1, D2) * sqrt(c(x$wt[x$i1], x$wo))
H = crossprod(D) # Choleski decomposition
v = sqrt(diag(H))
jv = v != 0
Hv = H[jv,jv,drop=FALSE] / tcrossprod(v[jv])
diag(Hv) = diag(Hv) + 1e-10
Rv = chol(Hv)
gv = dlogLik(h, x, expdH, interior=TRUE)[jv] / v[jv]
plus = forwardsolve(Rv, gv, upper.tri=TRUE, transpose=TRUE)
par = c(h$alpha, h$nu, h$mu[j2])[jv] * v[jv]
w.new = double(length(v))
w.new[jv] = nnls(Rv, Rv %*% par + plus)$x / v[jv]
alpha = w.new[1]
nu = if(k > 0) w.new[2:(k+1)] else numeric()
mu = if(m > 0) w.new[(k+2):(k+m+1)] else numeric()
newh = uh(alpha=alpha, tau=h$tau, nu=nu, eta=h$eta, mu=mu,
upper=x$upper, h$deg, collapse=FALSE)
if(h$deg == 0) b = backtrack(h, newh, x, expdH, alpha=0)
else b = backtrack(h, newh, x, expdH)
newh = b$h2
j1 = newh$nu != 0
j2 = newh$mu != 0
h2 = uh(newh$alpha, newh$tau[j1], newh$nu[j1], newh$eta[j2], newh$mu[j2],
upper=h$upper, h$deg, collapse=FALSE)
h = collapse(h2, x, tol=pmax(tol,1e-10))
list(h=h, convergence=b$convergence)
}
# Backtracking line search. h and h2 must have the same knots
backtrack = function(h, h2, x, expdH, tol=1e-10, alpha=0.33){
j = h$eta != h$upper
h2$eta = h$eta = h$eta[j]
h$mu = h$mu[j]
h2$mu = h2$mu[j]
ll.h = logLikuh(h, x)
d = c(h2$alpha - h$alpha, h2$nu - h$nu, h2$mu - h$mu)
g = alpha * sum(dlogLik(h, x, expdH) * d)
convergence = 0
r = 1
repeat {
hr = uh((1-r) * h$alpha + r * h2$alpha,
h$tau, (1-r) * h$nu + r * h2$nu,
h$eta, (1-r)*h$mu + r * h2$mu, upper=h$upper, deg=h$deg)
ll.hr = logLikuh(hr, x)
if (ll.hr >= ll.h + r * g) {convergence =0; break}
r = 0.5 * r
if(r < tol) {r = 0; hr = h; ll.h2 = ll.h; convergence = 1; break}
}
attr(h2, "ll") = ll.hr
list(h2=hr, r=r, convergence=convergence)
}
# Zeroth gradient function
grad0 = function(h, x, expdH=NULL) {
if(length(x$t) > 0)
d0 = sum(x$wt[x$i1] / hazuh(x$t[x$i1], h)) - sum(x$wt * x$t)
else d0 = 0
if(nrow(x$o) > 0) {
if(is.null(expdH)) expdH = exp(chazuh(x$o[,1],h) - chazuh(x$o[,2],h))
Delta = expdH / (1 - expdH)
xd = x$o[,1] - x$o[,2]
xd[Delta == 0] = 0
d0 = d0 - sum(x$wo * (x$o[,1] + xd * Delta))
}
d0
}
# First gradient function
grad1 = function(tau, h, x, expdH=NULL, order=0) {
g = vector("list", length(order))
names(g) = paste("d", order, sep="")
g[1:length(g)] = 0
if(length(tau) == 0) return(NULL)
m = length(tau)
n1 = length(x$t)
p = h$deg
if(n1 > 0) { # for exact observations
tau.r1 = rep(tau, rep.int(n1,m))
dim(tau.r1) = c(n1,m)
ind = tau.r1 >= x$t
tau.t = pmax(tau.r1 - x$t, 0)
if(0 %in% order) {
g$d0 = switch(as.character(p),
"0" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
ind[x$i1,,drop=FALSE]) - crossprod(x$wt, tau.r1 - tau.t))[1,],
"1" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
tau.t[x$i1,,drop=FALSE]) -
0.5 * crossprod(x$wt, tau.r1 * tau.r1 - tau.t * tau.t))[1,],
"2" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
tau.t[x$i1,,drop=FALSE]^2) -
crossprod(x$wt, tau.r1*tau.r1*tau.r1 - tau.t*tau.t*tau.t) / 3)[1,],
(crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
tau.t[x$i1,,drop=FALSE]^p) -
crossprod(x$wt, tau.r1^(p+1) - tau.t^(p+1))[1,] / (p+1))[1,] )
}
if(1 %in% order) {
g$d1 = switch(as.character(p),
"0" = double(m),
"1" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
ind[x$i1,,drop=FALSE]) - crossprod(x$wt, tau.r1 - tau.t))[1,],
"2" = (2 * crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
tau.t[x$i1,,drop=FALSE]) -
crossprod(x$wt, tau.r1 * tau.r1 - tau.t * tau.t))[1,],
{ tau.t.pm1 = tau.t[x$i1,,drop=FALSE]^(p-1)
if(p < 1) tau.t.pm1[tau.t.pm1 == Inf] = 0
(p * crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h), tau.t.pm1) -
crossprod(x$wt, tau.r1^p - tau.t^p))[1,]
} )
}
if(2 %in% order) {
g$d2 = switch(as.character(p),
"0" = double(m),
"1" = - sum(x$wt) + crossprod(x$wt, ind)[1,],
"2" = (2 * crossprod(x$wt[x$i1], ind[x$i1,,drop=FALSE]) -
crossprod(x$wt, tau.r1 - tau.t))[1,],
{ tau.t1.pm2 = tau.t[x$i1,,drop=FALSE]^(p-2)
if(p < 2) tau.t1.pm2[tau.t1.pm2 == Inf] = 0
tau.t.pm1 = tau.t^(p-1)
if(p < 1) tau.t.pm1[tau.t.pm1 == Inf] = 0
(p * (p-1) * crossprod(x$wt[x$i1], tau.t1.pm2) -
p * crossprod(x$wt, tau.r1^(p-1) - tau.t.pm1))[1,]
} )
}
}
n2 = nrow(x$o)
if(n2 > 0) { # for interval-censored observations
if(is.null(expdH)) expdH = exp(chazuh(x$o[,1],h) - chazuh(x$o[,2],h))
Delta = expdH / (1 - expdH)
tau.r2 = rep(tau, rep.int(n2,m))
dim(tau.r2) = c(n2,m)
tau.x1 = pmax(tau.r2 - x$o[,1], 0)
tau.x2 = pmax(tau.r2 - x$o[,2], 0)
ind1 = tau.r2 >= x$o[,1]
ind2 = tau.r2 >= x$o[,2]
if(0 %in% order) {
xd0 = switch(as.character(p),
"0" = {tau.r2.p1 = tau.r2; tau.x2 - (tau.x1.p1 <- tau.x1)},
"1" = {tau.r2.p1 = tau.r2 * tau.r2
tau.x2 * tau.x2 - (tau.x1.p1 <- tau.x1 * tau.x1)},
"2" = {tau.r2.p1 = tau.r2 * tau.r2 * tau.r2
tau.x2 * tau.x2 * tau.x2 -
(tau.x1.p1 <- tau.x1 * tau.x1 * tau.x1)},
{ tau.r2.p1 = tau.r2^(p+1)
tau.x2^(p+1) - (tau.x1.p1 <- tau.x1^(p+1))} )
xd0[Delta == 0] = 0
g$d0 = g$d0 -
crossprod(x$wo, (tau.r2.p1 - tau.x1.p1 + xd0 * Delta))[1,] / (p+1)
}
if(1 %in% order) {
xd1 = switch(as.character(p),
"0" = {tau.r2.p = 1; ind2 - (tau.x1.p <- ind1)},
"1" = {tau.r2.p = tau.r2; tau.x2 - (tau.x1.p <- tau.x1)},
"2" = {tau.r2.p = tau.r2 * tau.r2;
tau.x2 * tau.x2 - (tau.x1.p <- tau.x1 * tau.x1)},
{tau.r2.p = tau.r2^p; tau.x2^p - (tau.x1.p <- tau.x1^p)} )
xd1[Delta == 0] = 0
g$d1 = g$d1 - crossprod(x$wo, tau.r2.p - tau.x1.p + xd1 * Delta)[1,]
}
if(2 %in% order) {
g$d2 = switch(as.character(p),
"0" = double(m),
"1" = g$d2 - sum(x$wo) +
crossprod(x$wo, ind1 - (ind2 - ind1) * Delta)[1,],
"2" = { xd2 = tau.x2 - tau.x1
xd2[Delta == 0] = 0
g$d2 -
2 * crossprod(x$wo, tau.r2 - tau.x1 + xd2 * Delta)[1,] },
{ tau.x1.pm1 = tau.x1^(p-1)
if(p < 1) tau.x1.pm1[tau.x1.pm1 == Inf] = 0
xdp = tau.x2^(p-1) - tau.x1^(p-1)
xdp[Delta == 0] = 0
g$d2 - p * crossprod(x$wo,
tau.r2^(p-1) - tau.x1^(p-1) + xdp * Delta)[1,]
} )
}
}
g
}
# Second gradient function
grad2 = function(eta, h, x, expdH=NULL, order=0) {
g = vector("list", length(order))
names(g) = paste("d", order, sep="")
g[1:length(g)] = 0
if(length(eta) == 0) return(NULL)
m = length(eta)
n1 = length(x$t)
p = h$deg
if(n1 > 0) {
eta.r1 = rep(eta, rep.int(n1,m))
dim(eta.r1) = c(n1,m)
t.eta = pmax(x$t - eta.r1, 0)
ind = x$t >= eta.r1
if(0 %in% order) {
g$d0 = switch(as.character(p),
"0" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
ind[x$i1,,drop=FALSE]) - crossprod(x$wt, t.eta))[1,],
"1" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
t.eta[x$i1,,drop=FALSE]) - 0.5 * crossprod(x$wt, t.eta^2))[1,],
"2" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
t.eta[x$i1,,drop=FALSE]^2) -
crossprod(x$wt, t.eta * t.eta * t.eta) / 3)[1,],
(crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
t.eta[x$i1,,drop=FALSE]^p) -
crossprod(x$wt, t.eta^(p+1)) / (p+1))[1,] )
}
if(1 %in% order) {
g$d1 = switch(as.character(p),
"0" = double(m),
"1" = (-crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
ind[x$i1,,drop=FALSE]) +
crossprod(x$wt, t.eta))[1,],
"2" = (-2 * crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
t.eta[x$i1,,drop=FALSE]) +
crossprod(x$wt, t.eta^2))[1,],
{ t.eta.pm1 = t.eta[x$i1,,drop=FALSE]^(p-1)
if(p < 1) t.eta.pm1[t.eta.pm1 == Inf] = 0
(-p * crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h), t.eta.pm1) +
crossprod(x$wt, t.eta^p))[1,]
} )
}
if(2 %in% order) {
g$d2 = switch(as.character(p),
"0" = double(m),
"1" = - crossprod(x$wt, ind)[1,],
"2" = - 2 * crossprod(x$wt, t.eta)[1,],
{ t.eta.pm1 = t.eta^(p-1)
if(p < 1) t.eta.pm1[t.eta.pm1 == Inf] = 0
- p * crossprod(x$wt, t.eta.pm1)[1,]
} )
}
}
n2 = nrow(x$o)
if(n2 > 0) {
if(is.null(expdH)) expdH = exp(chazuh(x$o[,1],h) - chazuh(x$o[,2],h))
Delta = expdH / (1 - expdH)
eta.r2 = rep(eta, rep.int(n2,m))
dim(eta.r2) = c(n2,m)
x1.eta = pmax(x$o[,1] - eta.r2, 0)
x2.eta = pmax(x$o[,2] - eta.r2, 0)
ind1 = x$o[,1] >= eta.r2
ind2 = x$o[,2] >= eta.r2
if(0 %in% order) {
xd0 = switch(as.character(p),
"0" = (x1.eta.p1 <- x1.eta) - x2.eta,
"1" = (x1.eta.p1 <- x1.eta * x1.eta) - x2.eta * x2.eta,
"2" = (x1.eta.p1 <- x1.eta * x1.eta * x1.eta) -
x2.eta * x2.eta * x2.eta,
(x1.eta.p1 <- x1.eta^(p+1)) - x2.eta^(p+1) )
xd0[Delta == 0] = 0
g$d0 = g$d0 - crossprod(x$wo, x1.eta.p1 + xd0 * Delta)[1,] / (p+1)
}
if(1 %in% order) {
xd1 = switch(as.character(p),
"0" = (x1.eta.p <- ind1) - ind2,
"1" = (x1.eta.p <- x1.eta) - x2.eta,
"2" = (x1.eta.p <- x1.eta * x1.eta) - x2.eta * x2.eta,
(x1.eta.p <- x1.eta^p) - x2.eta^p )
xd1[Delta == 0] = 0
g$d1 = g$d1 + crossprod(x$wo, x1.eta.p + xd1 * Delta)[1,]
}
if(2 %in% order) {
g$d2 = switch(as.character(p),
"0" = double(m),
"1" = g$d2 - crossprod(x$wo, ind1 + (ind1 - ind2) * Delta)[1,],
"2" = { xd2 = x1.eta - x2.eta
xd2[Delta == 0] = 0
g$d2 - 2 * crossprod(x$wo, x1.eta + xd2 * Delta)[1,] },
{ x1.eta.pm1 = x1.eta^(p-1)
x2.eta.pm1 = x2.eta^(p-1)
if(p < 1) {
x1.eta.pm1[x1.eta.pm1 == Inf] = 0
x2.eta.pm1[x2.eta.pm1 == Inf] = 0
}
xdp = x1.eta.pm1 - x2.eta.pm1
xdp[Delta == 0] = 0
g$d2 - p * crossprod(x$wo, x1.eta.pm1 + xdp * Delta)[1,]
} )
}
}
g
}
##'Computes the Log-likelihood Value of a U-shaped Hazard Function
##'
##'
##'\code{logLikuh} returns the log-likelihood value of a U-shaped hazard
##'function, given a data set.
##'
##'
##'@param h an object of class \code{uh}.
##'@param data numeric vector or matrix for exact or interval-censored
##'observations, or an object of class \code{icendata}.
##'@return
##'
##'Log-likelihood value evaluated at \code{h}, given \code{data}.
##'@author Yong Wang <yongwang@@auckland.ac.nz>
##'@seealso \code{\link{Uhaz}}, \code{\link{icendata}}, \code{\link{plot.uh}}
##'@references
##'
##'Wang, Y. and Fani, S. (2018). Nonparametric maximum likelihood
##' computation of a U-shaped hazard function. \emph{Statistics and
##' Computing}, \bold{28}, 187-200.
##'
##'@keywords function
##'@examples
##'
##'data(ap)
##'(h0 = uh(.2, NULL, NULL, NULL, NULL, 15, 1)) # Uniform hazard
##'plot(h0, ylim=c(0,.3))
##'logLikuh(h0, ap)
##'
##'r = Uhaz(ap, deg=2)
##'r$ll
##'logLikuh(r$h, ap)
##'plot(r$h, add=TRUE, col="red3")
##'
##'@export logLikuh
logLikuh = function(h, data) {
x = icendata(data)
if(length(x$t) > 0) ll = sum(x$wt[x$i1] * log(hazuh(x$t[x$i1], h))) -
sum(x$wt * chazuh(x$t, h))
else ll = 0
if(nrow(x$o) > 0) {
ll = ll + sum(x$wo *
(log(exp(-chazuh(x$o[,1], h)) - exp(-chazuh(x$o[,2], h)))))
}
ll
}
dlogLik = function(h, x, expdH=NULL, interior=FALSE) {
if(interior) eta = h$eta[h$eta != h$upper]
else eta = h$eta
m = length(eta)
d = c(grad0(h, x, expdH), grad1(h$tau, h, x, expdH)$d0,
if(m>0) grad2(eta, h, x, expdH)$d0 else NULL)
names(d) = c("alpha", paste("nu",1:length(h$tau),sep=""),
if(m>0) paste("mu",1:m,sep="") else NULL)
d
}
### Finds all local maxima of the gradient functions
# bc boundary change
maxgrad = function(h, x, expdH=NULL, maxit=100, grid=100, bc=TRUE) {
if(!is.icendata(x)) x = icendata(x, w=1)
u = sort(unique(c(x$u, h$tau, h$eta)))
p = h$deg
if(p > 0 && p < 0.1)
u = rep(u, rep.int(21, length(u))) + seq(-h$upper*1e-2, h$upper*1e-2, len=21)
if(p == 1) maxit = 1
tau1 = c(0, h$tau[h$tau > 0])
k = length(tau1) - 1
eta1 = c(h$eta[h$eta<h$upper], h$upper)
m = length(eta1) - 1
if(p < 1) grid = grid * log2(2/(p*p))
## grad1
if(p < 0.1) { # optimal points
if(bc) pt1 = u[u<=eta1[1]]
else pt1 = u[u<eta1[1]]
gpt1 = grad1(pt1, h, x, expdH)$d0
}
else {
if(p == 1) u1 = u[u <= eta1[1]]
else u1 = seq(0, eta1[1], len=grid)
if(!bc) u1 = u1[u1 < eta1[1]]
m1 = length(u1)
pt1 = gpt1 = numeric()
if(p == 1) jd = rep(TRUE, m1-1)
else {
g = grad1(u1, h, x, expdH, order=1)
jd = g$d1[-m1] > 0 & g$d1[-1] < 0
}
if(any(jd)) {
left = u1[-m1][jd]
right = u1[-1][jd]
pt1 = (left + right) * .5
for(i in 1:maxit) {
g = grad1(pt1, h, x, expdH, order=1:2)
left[g$d1>0] = pt1[g$d1>0]
right[g$d1<0] = pt1[g$d1<0]
pt1.old = pt1
pt1 = pt1 - g$d1 / g$d2
j = is.na(pt1) | pt1 < left | pt1 > right
pt1[j] = (left[j] + right[j]) * .5
if( max(abs(pt1 - pt1.old)) <= 1e-14 * h$upper) break
}
if(p == 1) pt1 = pt1[!j]
gpt1 = grad1(pt1, h, x, expdH, order=0)$d0
}
}
i = pt1 > 0 & pt1 <= tau1[k+1]
pt1i = pt1[i]
gpt1i = gpt1[i]
if(k > 0 && length(pt1i) > 0) {
grp = apply(outer(pt1i, tau1[-(k+1)], ">") & outer(pt1i, tau1[-1], "<="),
1, which.max)
r1 = aggregate(gpt1i, by=list(group=grp), which.max)
j = integer(k)
j[r1[,1]] = r1[,2]
j = j + c(0,cumsum(tabulate(grp, nbins=k))[-k])
np1 = pt1i[j]
gnp1 = gpt1i[j]
j0 = gnp1 > 0
np1 = np1[j0]
gnp1 = gnp1[j0]
}
else np1 = gnp1 = numeric()
## grad2
if(p < 0.1) {
if(bc) pt2 = u[u>=tau1[k+1] & u<=h$upper]
else pt2 = u[u>tau1[k+1] & u<h$upper]
gpt2 = grad2(pt2, h, x, expdH)$d0
}
else {
if(p == 1) u2 = u[u >= tau1[k+1]]
else u2 = seq(tau1[k+1], h$upper, len=grid)
if(!bc) u2 = u2[u2 > tau1[k+1]]
m2 = length(u2)
pt2 = gpt2 = numeric()
if(p == 1) jd = rep(TRUE, m2-1)
else {
g = grad2(u2, h, x, expdH, order=1)
jd = g$d1[-m2] > 0 & g$d1[-1] < 0
}
if(any(jd)) {
left = u2[-m2][jd]
right = u2[-1][jd]
pt2 = (left + right) * .5
for(i in 1:maxit) {
g = grad2(pt2, h, x, expdH, order=1:2)
left[g$d1>0] = pt2[g$d1>0]
right[g$d1<0] = pt2[g$d1<0]
pt2.old = pt2
pt2 = pt2 - g$d1 / g$d2
j = is.na(pt2) | pt2 < left | pt2 > right
pt2[j] = (left[j] + right[j]) * .5
if( max(abs(pt2 - pt2.old)) <= 1e-14 * h$upper) break
}
if(p == 1) pt2 = pt2[!j]
gpt2 = grad2(pt2, h, x, expdH, order=0)$d0
}
}
i = pt2 >= eta1[1] & pt2 < eta1[m+1]
pt2i = pt2[i]
gpt2i = gpt2[i]
if(m > 0 && length(pt2i) > 0) {
grp = apply(outer(pt2i, eta1[-(m+1)], ">=") & outer(pt2i, eta1[-1], "<"), 1,
which.max)
r2 = aggregate(gpt2i, by=list(group=grp), which.max)
j = integer(m)
j[r2[,1]] = r2[,2]
j = j + c(0, cumsum(tabulate(grp, nbins=m))[-m])
np2 = pt2i[j]
gnp2 = gpt2i[j]
j0 = gnp2 > 0
np2 = np2[j0]
gnp2 = gnp2[j0]
}
else np2 = gnp2 = numeric()
## grad1 and grad2
if(max(h$tau) != h$eta[1]) {
jj1 = pt1 >= tau1[k+1] & pt1 <= eta1[1]
if(p == 0) { uj1 = pt1[jj1]; gj1 = gpt1[jj1] }
else {
uj1 = c(tau1[k+1], if(bc) eta1[1] else NULL, pt1[jj1])
gj1 = c(grad1(c(tau1[k+1], if(bc) eta1[1] else NULL),
h, x, expdH)$d0, gpt1[jj1])
}
jmax = which.max(gj1)
np31 = uj1[jmax]
gnp31 = gj1[jmax]
jj2 = pt2 >= tau1[k+1] & pt2 <= eta1[1]
if(p == 0) { uj2 = pt2[jj2]; gj2 = gpt2[jj2] }
else {
uj2 = c(if(bc) tau1[k+1] else NULL, eta1[1], pt2[jj2])
gj2 = c(grad2(c(if(bc) tau1[k+1] else NULL, eta1[1]),
h, x, expdH)$d0, gpt2[jj2])
}
jmax = which.max(gj2)
np32 = uj2[jmax]
gnp32 = gj2[jmax]
if(gnp31 > gnp32) {np1 = c(np1, np31); gnp1 = c(gnp1, gnp31)}
else {np2 = c(np2, np32); gnp2 = c(gnp2, gnp32)}
}
list(np1=np1, gnp1=gnp1, np2=np2, gnp2=gnp2)
}
simplify = function(h) {
i1 = order(h$tau) # remove identical knots
tau = h$tau[i1]
nu = h$nu[i1]
i2 = order(h$eta)
eta = h$eta[i2]
mu = h$mu[i2]
if(h$deg != 0 || length(tau) == 0 || length(eta) == 0) return(h)
if(tau[length(tau)] != eta[1]) return(h)
if(nu[length(tau)] < mu[1]) {
if(eta[1] == 0) {
h$alpha = h$alpha + mu[1]
eta = eta[-1]
mu = mu[-1]
}
else {
h$alpha = h$alpha + nu[length(nu)]
mu[1] = mu[1] - nu[length(nu)]
tau = tau[-length(tau)]
nu = nu[-length(nu)]
}
}
else {
h$alpha = h$alpha + mu[1]
nu[length(nu)] = nu[length(nu)] - mu[1]
eta = eta[-1]
mu = mu[-1]
}
uh(h$alpha, tau, nu, eta, mu, h$upper, h$deg, collapse=FALSE)
}
# Collapse similar knots
collapse = function(h, x, tol=0){
ll = attr(h, "ll")
i1 = order(h$tau) # remove identical knots
tau = h$tau[i1]
ind1 = cumsum(!duplicated(c(tau)))
tau = unique(tau)
nu = aggregate(h$nu[i1], by=list(group=ind1), sum)[,2]
nu[tau == 0] = 0
i2 = order(h$eta)
eta = h$eta[i2]
ind2 = cumsum(!duplicated(eta))
eta = unique(eta)
mu = aggregate(h$mu[i2], by=list(group=ind2), sum)[,2]
mu[eta == h$upper] = 0
h = uh(h$alpha, tau, nu, eta, mu, h$upper, h$deg, collapse=FALSE)
if(tol > 0) {
if(is.null(ll)) ll = logLikuh(h, x)
# if(h$deg < 1) {attr(h, "ll") = ll; return(h)} # why?
h2 = h
break1 = break2 = FALSE
repeat {
if(!break1 && length(h2$nu) > 1) {
j = which.min(diff(h$tau))
h2$nu[j] = h2$nu[j] + h2$nu[j+1]
h2$nu = h2$nu[-(j+1)]
h2$tau[j] = (h2$tau[j] + h2$tau[j+1]) * .5
h2$tau = h2$tau[-(j+1)]
ll2 = logLikuh(h2, x)
if(ll2 + tol >= ll) {h = h2; ll = ll2; break1 = FALSE}
else {h2 = h; break1 = TRUE}
}
else break1 = TRUE
if(!break2 && length(h2$mu) > 1) {
j = which.min(diff(h2$eta))
h2$mu[j] = h2$mu[j] + h2$mu[j+1]
h2$mu = h2$mu[-(j+1)]
h2$eta[j] = (h2$eta[j] + h2$eta[j+1]) * .5
h2$eta = h2$eta[-(j+1)]
ll2 = logLikuh(h2, x)
if(ll2 + tol >= ll) {h = h2; ll = ll2; break2 = FALSE}
else {h2 = h; break2 = TRUE}
}
else break2 = TRUE
if(break1 && break2) break
}
attr(h, "ll") = ll
}
h
}
# deg - polynomial degree
##'U-shaped Hazard Function
##'
##'
##' Class \code{uh} can be used to store U-shaped hazard functions.
##' There are a couple of functions associated with the class.
##'
##'\code{uh} creates an object of class \code{uh}, which stores a U-shaped
##'hazard function.
##'
##'\code{print.uh} prints an object of class \code{uh}.
##'
##'
##'A U-shape hazard function, as generalized by Wang and Fani (2018), is given
##'by
##'
##'\deqn{h(t) = \alpha + \sum_{j = 1}^k \nu_j(\tau_j - t)_+^p + \sum_{j = 1}^{m} \mu_j (t-\eta_j)_+^p,}{ h(t) = alpha + sum_{j=1}^k nu_j (tau_j - t)_+^p
##' + sum_{j=1}^m mu_j (t - eta_j)_+^p,}
##'
##'where \eqn{\alpha,\nu_j,\mu_j \ge 0}{alpha, nu_j, mu_j \ge 0},
##'\eqn{\tau_1 < \cdots < \tau_k \le \eta_1 < \cdots < \eta_m,}{tau_1 < ... < tau_k <= eta_1 < ... < eta_m,} and \eqn{p \ge 0}{p >= 0} is the the spline degree which
##'determines the smoothness of the U-shaped hazard. As p increases, the family
##'of hazard functions becomes increasingly smoother, but at the same time,
##'smaller. When \eqn{p = 0}{p = 0}, the hazard function is U-shaped, as
##'studied by Bray et al. (1967). When \eqn{p = 1}{p = 1}, the hazard function
##'is convex, as studied by Jankowski and Wellner (2009a,b).
##'
##'\code{print.uh} prints an object of class \code{uh}. While \code{alpha},
##'\code{upper} and \code{deg} are printed as they are, \code{tau} and
##'\code{nu} are printed as a two-column matrix, and so are \code{eta} and
##'\code{mu}.
##'
##'@aliases uh uh.object print.uh
##'@param alpha a nonnegative value, for the constant coefficient.
##'@param tau vector of nonnegative real values, for left knots.
##'@param nu vector of nonnegative values, for masses associated with the left
##'knots.
##'@param eta vector of nonnegative real values, for right knots.
##'@param mu vector of nonnegative real values, for masses associated with the
##'right knots.
##'@param upper a positive value, at which point the hazard starts to become
##'infinite.
##'@param deg nonnegative real number for spline degree (i.e., p in the formula
##'below).
##'@param collapse logical, indicating if identical knots should be collapsed.
##'@param x an object of class \code{uh}.
##'@param ... other auguments for printing.
##'@return
##'
##'\code{uh} returns an object of class \code{uh}. It is a list with components
##'\code{alpha}, \code{tau}, \code{nu}, \code{eta}, \code{mu}, \code{upper} and
##'\code{deg}, which store their corresponding values as described above.
##'@author Yong Wang <yongwang@@auckland.ac.nz>
##'@seealso \code{\link{Uhaz}}, \code{\link{icendata}}, \code{\link{plot.uh}}
##'@references
##'
##'Bray, T. A., Crawford, G. B., and Proschan, F. (1967). \emph{Maximum
##'Likelihood Estimation of a U-shaped Failure Rate Function}. Defense
##'Technical Information Center.
##'
##'Jankowski, H. K. and Wellner, J. A. (2009a). Computation of nonparametric
##'convex hazard estimators via profile methods. \emph{Journal of Nonparametric
##'Statistics}, \bold{21}, 505-518.
##'
##'Jankowski, H. K. and Wellner, J. A. (2009b). Nonparametric estimation of a
##'convex bathtub-shaped hazard function. \emph{Bernoulli}, \bold{15},
##'1010-1035.
##'
##'Wang, Y. and Fani, S. (2018). Nonparametric maximum likelihood
##' computation of a U-shaped hazard function. \emph{Statistics and
##' Computing}, \bold{28}, 187-200.
##'
##'@keywords function
##'@examples
##'
##'(h0 = uh(3, 2, 3, 4, 5, 7, deg=0)) # deg = 0
##'plot(h0, ylim=c(0,20))
##'(h1 = uh(4, 2, 3, 5, 6, 7, deg=1)) # deg = 1
##'plot(h1, add=TRUE, col="green3")
##'(h2 = uh(1, 1:2, 3:4, 5:6, 7:8, 9, deg=2)) # deg = 2
##'plot(h2, add=TRUE, col="red3")
##'
##'@usage
##'uh(alpha, tau, nu, eta, mu, upper=Inf, deg=1, collapse=TRUE)
##'\method{print}{uh}(x, ...)
##'
##'@export uh
##'@export print.uh
uh = function(alpha, tau, nu, eta, mu, upper=Inf, deg=1, collapse=TRUE) {
if(length(tau) == 0) {tau=0; nu=0}
if(length(eta) == 0) {eta=upper; mu=0}
i1 = order(tau)
tau = tau[i1]
nu = nu[i1]
i2 = order(eta)
eta = eta[i2]
mu = mu[i2]
h = list(alpha=alpha, tau=tau, nu=nu, eta=eta, mu=mu, upper=upper, deg=deg)
if(collapse) h = collapse(h)
class(h) = "uh"
h
}
print.uh = function(x, ...) {
cat("$alpha\n")
print(x$alpha, ...)
print(cbind(tau=x$tau, nu=x$nu), ...)
print(cbind(eta=x$eta, mu=x$mu), ...)
cat("$upper\n")
print(x$upper, ...)
cat("$deg\n")
print(x$deg, ...)
}
# Hazard function
##'Distributional Functions given a U-shaped Hazard Function
##'
##'
##'Given an object of class \code{uh}:
##'
##'\code{hazuh} computes the hazard values;
##'
##'\code{chazuh} computes the cumulative hazard values;
##'
##'\code{survuh} computes the survival function values;
##'
##'\code{denuh} computes the density function values.
##'
##'
##'@aliases hazuh chazuh survuh denuh
##'@param t time points at which the function is to be evaluated.
##'@param h an object of class \code{uh}.
##'@return
##'
##'A numeric vector of the function values.
##'@author Yong Wang <yongwang@@auckland.ac.nz>
##'@seealso \code{\link{Uhaz}}, \code{\link{icendata}}, \code{\link{plot.uh}}
##'@references
##'
##'Wang, Y. and Fani, S. (2018). Nonparametric maximum likelihood
##' computation of a U-shaped hazard function. \emph{Statistics and
##' Computing}, \bold{28}, 187-200.
##'
##'@keywords function
##'@examples
##'
##'data(ap)
##'h = Uhaz(icendata(ap), deg=2)$h
##'hazuh(0:15, h) # hazard
##'chazuh(0:15, h) # cumulative hazard
##'survuh(0:15, h) # survival probability
##'denuh(0:15, h) # density
##'
##'@usage
##'hazuh(t, h)
##'chazuh(t, h)
##'survuh(t, h)
##'denuh(t, h)
##'
##'@export hazuh
##'@export chazuh
##'@export survuh
##'@export denuh
hazuh = function(t, h) {
p = h$deg
b = c = 0
if(length(h$tau) > 0) {
if(p == 0) d = outer(h$tau, t, ">=")
else {
d = pmax(outer(h$tau, t, "-"), 0)
if(p != 1) d = d^p
}
b = (h$nu %*% d)[1,]
}
if(length(h$eta) > 0) {
if(p == 0) d = outer(t, h$eta, ">=")
else {
d = pmax(outer(t, h$eta, "-"), 0)
if(p != 1) d = d^p
}
c = (d %*% h$mu)[,1]
}
h$alpha + pmax(b, c)
}
chazuh = function(t, h) {
deg = pmax(h$deg, 1)
p1 = h$deg + 1
a = b = c = 0
if(h$alpha > 0) a = h$alpha * t
if(length(h$tau) > 0) {
tau.t = pmax(outer(h$tau, t, "-"), 0)
b = (h$nu %*% (h$tau^p1 - tau.t^p1))[1,] / p1
}
if(length(h$eta) > 0) {
t.eta = pmax(outer(t, h$eta, "-"), 0)
c = (t.eta^p1 %*% h$mu)[,1] / p1
}
H = a + b + c
H[t > h$upper] = 1e100
H
}
# survival function
survuh = function(t, h) exp(-chazuh(t, h))
# density function
denuh = function(t, h) hazuh(t, h) * survuh(t, h)
## plotting functions
##'Plot Functions for U-shaped Hazard Estimation
##'
##'
##' Functions for plotting various functions in U-shaped hazard estimation
##'
##'\code{plot.Uhaz} and \code{plot.uh} are wrapper functions that can be used
##'to invoke \code{plot.hazuh}, \code{plot.chazuh}, \code{plot.survuh},
##'\code{plot.denuh} or \code{plot.graduh}.
##'
##'\code{plothazuh} plots a U-shaped hazard function.
##'
##'\code{plotchazuh} plots a cumulative hazard function that has a U-shaped
##'hazard function.
##'
##'\code{plotsurvuh} plots the survival function that has a U-shaped hazard
##'function.
##'
##'\code{plotdenuh} plots the density function that has a U-shaped hazard
##'function.
##'
##'\code{plotgraduh} plots the gradient function that has a U-shaped hazard
##'function.
##'
##'
##'A U-shaped hazard function is given by
##'
##'\deqn{h(t) = \alpha + \sum_{j = 1}^k \nu_j(\tau_j - t)_+^p + \sum_{j = 1}^{m} \mu_j (t-\eta_j)_+^p,}{ h(t) = alpha + sum_{j=1}^k nu_j (tau_j - t)_+^p
##' + sum_{j=1}^m mu_j (t - eta_j)_+^p,}
##'
##'where \eqn{\alpha,\nu_j,\mu_j \ge 0}{alpha, nu_j, mu_j \ge 0},
##'\eqn{\tau_1 < \cdots < \tau_k \le \eta_1 < \cdots < \eta_m,}{tau_1 < ... < tau_k <= eta_1 < ... < eta_m,} and \eqn{p \ge 0}{p >= 0}.
##'
##'@aliases plot.Uhaz plot.uh plothazuh plotchazuh plotsurvuh plotdenuh
##'plotgraduh
##'@param x an object of class \code{Uhaz}, i.e., an output of function
##'\code{Uhaz}, or an object of class \code{uh}..
##'@param h an object of class \code{uh}.
##'@param data vector or matrix that stores observations, or an object of class
##'\code{icendata}.
##'@param w additional weights/multiplicities for the observations stored in
##'\code{data}.
##'@param fn function to be plotted. It can be
##'
##'= \code{haz}, for hazard function;
##'
##'= \code{chaz}, for cumulative hazard function;
##'
##'= \code{den}, for density function;
##'
##'= \code{surv}, for survival function;
##'
##'= \code{gradient}, for gradient functions.
##'
##'@param xlim,ylim numeric vectors of length 2, giving the x and y coordinates
##'ranges.
##'@param xlab,ylab x- or y-axis labels.
##'@param add = \code{TRUE}, adds the curve to the existing plot;
##'
##'= \code{FALSE}, plots the curve in a new one.
##'@param col color used for plotting the curve.
##'@param lty line type for plotting the curve.
##'@param lwd line width for plotting the curve.
##'@param len number of points used to plot a curve.
##'@param add.knots logical, indicating if knots are also plotted.
##'@param pch point character/type for plotting knots.
##'@param vert logical, indicating if grey vertical lines are plotted to show
##'the interval that separates the two discrete measures.
##'@param col0 color for gradient function 0, i.e., for the hazard-constant
##'part, or alpha.
##'@param col1 color for gradient function 1, i.e., for the hazard-decreasing
##'part.
##'@param col2 color for gradient function 1, i.e., for the hazard-increasing
##'part.
##'@param order = 0, the gradient functions are plotted;
##'
##'= 1, their first derivatives are plotted;
##'
##'= 2, their second derivatives are plotted.
##'@param ... arguments for other graphical parameters (see \code{par}).
##'@author Yong Wang <yongwang@@auckland.ac.nz>
##'@seealso \code{\link{icendata}}, \code{\link{uh}}, \code{\link{npsurv}}.
##'@references
##'
##'Wang, Y. and Fani, S. (2018). Nonparametric maximum likelihood
##' computation of a U-shaped hazard function. \emph{Statistics and
##' Computing}, \bold{28}, 187-200.
##'
##'@keywords function
##'@examples
##'
##'## Angina Pectoris Survival Data
##'data(ap)
##'plot(r<-Uhaz(ap)) # hazard function for a convex hazard
##'plot(r, fn="c") # cumulative hazard function
##'plot(r, fn="s") # survival function
##'plot(r, fn="d") # density function
##'plot(r, ap, fn="g") # gradient functions
##'plot(r, ap, fn="g", order=1) # first derivatives of gradient functions
##'plot(r, ap, fn="g", order=2) # second derivatives of gradient functions
##'
##'## New Zealand Mortality in 2000
##'data(nzmort)
##'i = nzmort$ethnic == "maori"
##'x = nzmort[i,1:2] # Maori mortality
##'h = Uhaz(x[,1]+0.5, x[,2], deg=2)$h # smooth U-shaped hazard
##'plot(h) # hazard function
##'plot(h, fn="d") # density function
##'plot(h, fn="s") # survival function
##'
##'x2 = nzmort[!i,1:2] # Non-Maori mortality
##'h2 = Uhaz(x2[,1]+0.5, x2[,2], deg=2)$h
##'plot(h2, fn="s", add=TRUE, col="green3")
##'
##'@usage
##'\method{plot}{Uhaz}(x, ...)
##'\method{plot}{uh}(x, data, fn=c("haz","grad","surv","den","chaz"), ...)
##'plothazuh(h, add=FALSE, col="darkblue", lty=1, xlim, ylim,
##' lwd=2, pch=19, len=500, vert=FALSE, add.knots=TRUE,
##' xlab="Time", ylab="Hazard", ...)
##'plotchazuh(h, add=FALSE, lwd=2, len=500, col="darkblue",
##' pch=19, add.knots=TRUE, vert=FALSE, xlim, ylim, ...)
##'plotdenuh(h, add=FALSE, lty=1, lwd=2, col="darkblue",
##' add.knots=TRUE, pch=19, ylim, len=500, vert=FALSE, ...)
##'plotsurvuh(h, add=FALSE, lty=1, lwd=2, len=500, vert=FALSE,
##' col="darkblue", pch=19, add.knots=TRUE, xlim, ylim, ...)
##'plotgraduh(h, data, w=1, len=500, xlim, ylim, vert=TRUE,
##' add=FALSE, xlab="Time", ylab="Gradient",
##' col0="red3", col1="blue3", col2="green3", order=0, ...)
##'
##'@export plot.Uhaz
##'@export plot.uh
##'@export plothazuh
##'@export plotchazuh
##'@export plotsurvuh
##'@export plotdenuh
##'@export plotgraduh
plot.Uhaz = function(x, ...) plot(x$h, ...)
plot.uh = function(x, data, fn=c("haz","grad","surv","den","chaz"), ...) {
fn = match.arg(fn)
fnR = getFunction(paste("plot",fn,"uh",sep=""))
switch(fn,
"haz" =,
"surv" =,
"den" =,
"chaz" = fnR(x, ...),
"grad" = fnR(x, data, ...) )
}
plothazuh = function(h, add=FALSE, col="darkblue", lty=1,
xlim, ylim, lwd=2, pch=19, len=500, vert=FALSE, add.knots=TRUE,
xlab="Time", ylab="Hazard", ...) {
p = h$deg
pt = switch(as.character(p),
"0" = unique(sort(c(0, h$tau, h$eta, h$upper))),
"1" = unique(sort(c(0, h$tau, h$eta, h$upper))),
unique(sort(c(h$tau, h$eta, seq(0, h$upper, len=len)))))
m = length(pt)
knots = unique(c(h$tau, h$eta))
haz = hazuh(pt, h)
max.haz = max(haz)
if(missing(xlim)) xlim = range(pt)
if(missing(ylim)) ylim = c(0, max.haz)
if(!add) plot(pt, haz, type="n", xlim=xlim, ylim=ylim,
xlab=xlab, ylab=ylab, ...)
if(vert) {
lines(rep(max(h$tau),2), ylim, col="grey", lty=2)
lines(rep(min(h$eta),2), ylim, col="grey", lty=2)
}
abline(h=0, col ="black")
if(p == 0) {
lines(c(h$tau[length(h$tau)], h$eta[1]), rep(h$alpha,2),
lwd=lwd, col=col, lty=lty)
lines(c(rep(rev(h$tau),each=2), 0),
c(h$alpha, rep(hazuh(rev(h$tau), h), each=2)),
lwd=lwd, col=col, lty=lty)
lines(c(rep(h$eta,each=2), h$upper),
c(h$alpha, rep(hazuh(h$eta, h), each=2)),
lwd=lwd, col=col, lty=lty)
}
else lines(pt, haz, lwd=lwd, col=col, lty=lty)
if(add.knots && length(knots) > 0)
points(knots, hazuh(knots, h), col=col, pch=pch)
}
plotchazuh = function(h, add=FALSE, lwd=2, len=500, col="darkblue",
pch=19, add.knots=TRUE, vert=FALSE, xlim, ylim, ...) {
pt = unique(sort(c(seq(0, h$upper, len=len), h$tau, h$eta)))
m = length(pt)
H = chazuh(pt, h)
max.H = max(H)
if(missing(xlim)) xlim = range(pt)
if(missing(ylim)) ylim = c(0, max.H)
plot(rep(max(h$tau),2), c(0,max.H), type="n",
xlim=xlim, ylim=ylim, xlab="Time", ylab="Cumulative Hazard", ...)
if(vert) lines(rep(max(h$tau),2), c(0,max.H), type="l", col ="grey", lty=2)
if(vert) lines(rep(min(h$eta),2), c(0,max.H), col ="grey", lty=2)
abline(h=0, col ="black")
lines(pt, H, type="l", lwd=lwd, col=col)
if(add.knots) {
knots = c(h$tau, h$eta)
knots = knots[knots>0 & knots<h$upper]
points(knots, chazuh(knots, h), col=col, pch=pch)
}
}
plotdenuh = function(h, add=FALSE, lty=1, lwd=2,
col="darkblue", add.knots=TRUE, pch=19,
ylim, len=500, vert=FALSE, ...) {
pt = sort(c(seq(0, h$upper, len=len), h$tau, h$eta))
m = length(pt)
d = denuh(pt, h)
max.d = max(d)
type = if(vert) "l" else "n"
if(missing(ylim)) ylim = c(0, max.d)
if(add) lines(rep(max(h$tau),2), ylim, type=type, col ="grey", lty=2)
else plot(rep(max(h$tau),2), ylim, type=type,
col ="grey", lty=2, xlim=range(pt), ylim=ylim,
xlab="Time", ylab="Density", ...)
if(vert) lines(rep(min(h$eta),2), ylim, col="grey", lty=2)
abline(h=0, col ="black")
lines(pt, d, lty=lty, lwd=lwd, col=col)
if(add.knots) {
knots = c(h$tau, h$eta)
knots = knots[knots>0 & knots<h$upper]
points(knots, denuh(knots, h), col=col, pch=pch)
}
}
plotsurvuh = function(h, add=FALSE, lty=1, lwd=2, len=500, vert=FALSE,
col="darkblue", pch=19, add.knots=TRUE, xlim, ylim, ...) {
pt = sort(c(seq(0, h$upper, len=len), h$tau, h$eta))
m = length(pt)
s = survuh(pt, h)
if(missing(xlim)) xlim = range(pt)
if(missing(ylim)) ylim = c(0, 1)
if(add) {
if(vert) lines(rep(max(h$tau),2), c(0,1), col ="grey", lty=2)
}
else {
if(vert)
plot(rep(max(h$tau),2), c(0,1), type="l", col ="grey", lty=2,
xlim=xlim, ylim=ylim, xlab="Time", ylab="Survival Probability",
...)
else plot(0, 0, type="n", xlim=xlim, ylim=ylim,
xlab="Time", ylab="Survival Probability")
}
if(vert) lines(rep(min(h$eta),2), c(0,1), col ="grey", lty=2)
abline(h=0, col ="black")
lines(pt, s, lty=lty, lwd=lwd, col=col)
pt = c(h$tau, h$eta)
if(add.knots) {
knots = c(h$tau, h$eta)
knots = knots[knots>0 & knots<h$upper]
points(knots, survuh(knots, h), col=col, pch=pch)
}
}
plotgraduh = function(h, data, w=1, len=500, xlim, ylim, vert=TRUE,
add=FALSE, xlab="Time", ylab="Gradient",
col0="red3", col1="blue3", col2="green3", order=0, ...) {
if(! order %in% 0:2) stop("'order' must be 0, 1 or 2")
x = icendata(data, w)
maxf = h$upper
if(missing(xlim)) xlim = c(0, maxf)
k = length(h$tau)
tau1 = unique(sort(c(seq(xlim[1], h$eta[1] - if(h$deg!=0) 0 else 1e-4*h$upper,
len=len), h$tau)))
eta1 = unique(sort(c(seq(h$tau[k] + if(h$deg!=0) 0 else 1e-4*h$upper, xlim[2],
len=len), h$eta)))
g11 = grad1(tau1, h, x, order=order)[[1]]
g22 = grad2(eta1, h, x, order=order)[[1]]
g0 = if(order == 0) grad0(h, x) else 0
if(missing (ylim)) ylim= range(g0, g11, g22)
if(add) lines(c(0,maxf), c(0, 0), col="grey")
else plot(c(0,maxf), c(0, 0), type="l", col="grey", xlim=xlim, ylim=ylim,
xlab=xlab, ylab=ylab, ...)
if(vert) {
lines(rep(h$tau[k],2), ylim, col="gray", lty=2)
lines(rep(h$eta[1],2), ylim, col="gray", lty=2)
}
lines(eta1, g22, col=col2, lty=1, lwd=1)
lines(tau1, g11, col=col1, lty=1, lwd=1)
j1 = h$tau != 0 & h$tau != h$upper
j2 = h$eta != 0 & h$eta != h$upper
points(h$tau[j1], grad1(h$tau[j1],h,x,order=order)[[1]], col=col1, pch=2)
points(h$eta[j2], grad2(h$eta[j2],h,x,order=order)[[1]], col=col2, pch=6)
lines(c(h$tau[k], h$eta[1]), rep(g0, 2), col=col0, lty=1, lwd=1)
points((h$tau[k] + h$eta[1])/2, g0, col=col0, pch=20)
if(!add) legend("bottomright",
legend=expression(d[0], d[1], d[2]), col=c(col0, col1, col2),
lty=c(1,1,1), pch=c(20,2,6), bty="n", lwd=c(1,1,1))
}
uh.initial = function(x, deg=1) {
if(!is.icendata(x)) x = icendata(x)
n1 = length(x$t)
n2 = nrow(x$o)
y1 = if(n1 > 0) x$t else numeric()
y2 = if(n2 > 0) {
x$o[x$o[,2] == Inf,2] = 1.6 * x$upper
rowMeans(x$o)
}
else numeric()
beta = weighted.mean(c(y1,y2), c(x$wt, x$wo))
uh(alpha=1/beta, tau=NULL, nu=NULL, eta=NULL, mu=NULL, upper=x$upper,
deg=deg)
}
Try the npsurv package in your browser
Any scripts or data that you put into this service are public.
npsurv documentation built on Oct. 23, 2020, 5:43 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions uh.initial plotgraduh plotsurvuh plotdenuh plotchazuh plothazuh plot.uh plot.Uhaz denuh survuh chazuh hazuh print.uh uh collapse simplify maxgrad dlogLik logLikuh grad2 grad1 grad0 backtrack updatemass Uhaz
Documented in chazuh denuh hazuh logLikuh plotchazuh plotdenuh plotgraduh plothazuh plotsurvuh plot.uh plot.Uhaz print.uh survuh uh Uhaz
############################################
# Estimation of a U-shaped Hazard Function #
############################################
##'U-shaped Hazard Function Estimation
##'
##'
##'\code{Uhaz} computes the nonparametric maximum likelihood esimate (NPMLE) of
##'a U-shaped hazard function from exact or interval-censored data, or a mix of
##'the two types of data.
##'
##'
##'If \code{data} is a vector, it contains only exact observations, with
##'weights given in \code{w}.
##'
##'If \code{data} is a matrix with two columns, it contains interval-censored
##'observations, with the two columns storing their left and right end-points,
##'respectively. If the left and right end-points are equal, then the
##'observation is exact. Weights are provided by \code{w}.
##'
##'If \code{data} is a matrix with three columns, it contains interval-censored
##'observations, with the first two columns storing their left and right
##'end-points, respectively. The weight of each observation is the third-column
##'value multiplied by the corresponding weight value in \code{w}.
##'
##'The algorithm used for the computing the NPMLE of a hazard function under
##'the U-shape restriction is is proposed by Wang and Fani (2015). Such a
##'hazard function is given by
##'
##'A U-shaped hazard function is given by
##'
##'\deqn{h(t) = \alpha + \sum_{j = 1}^k \nu_j(\tau_j - t)_+^p + \sum_{j = 1}^{m} \mu_j (t-\eta_j)_+^p,}{ h(t) = alpha + sum_{j=1}^k nu_j (tau_j - t)_+^p
##' + sum_{j=1}^m mu_j (t - eta_j)_+^p,}
##'
##'where \eqn{\alpha,\nu_j,\mu_j \ge 0}{alpha, nu_j, mu_j \ge 0},
##'\eqn{\tau_1 < \cdots < \tau_k \le \eta_1 < \cdots < \eta_m,}{tau_1 < ... < tau_k <= eta_1 < ... < eta_m,} and \eqn{p \ge 0}{p >= 0} is the the spline degree which
##'determines the smoothness of the U-shaped hazard. As p increases, the family
##'of hazard functions becomes increasingly smoother, but at the time, smaller.
##'When p = 0, the hazard function is U-shaped, as studied by Bray et al.
##'(1967). When p = 1, the hazard function is convex, as studied by Jankowski
##'and Wellner (2009a,b).
##'
##'Note that \code{deg} (i.e., p in the above mathematical display) can take on
##'any nonnegative real value.
##'
##'@aliases Uhaz Uhaz.object
##'@param data vector or matrix, or an object of class \code{icendata}.
##'@param w weights or multiplicities of the observations.
##'@param deg nonnegative real number for spline degree (i.e., p in the formula
##'below).
##'@param maxit maximum number of iterations.
##'@param tol tolerance level for stopping the algorithm. It is used as the
##'threshold on the increase of the log-likelihood after each iteration.
##'@param verb verbosity level for printing intermediate results in each
##'iteration.
##'@return
##'
##'An object of class \code{Uhaz}, which is a list with components:
##'
##'\item{convergence}{= \code{TRUE}, converged successfully;
##'
##'= \code{FALSE}, maximum number of iterations reached.}
##'
##'\item{grad}{gradient values at the knots.}
##'
##'\item{numiter}{number of iterations used.}
##'
##'\item{ll}{log-likelihood value of the NPMLE \code{h}.}
##'
##'\item{h}{NPMLE of the U-shaped hazard function, an object of class
##'\code{uh}.}
##'@author Yong Wang <yongwang@@auckland.ac.nz>
##'@seealso \code{\link{icendata}}, \code{\link{nzmort}}.
##'@references
##'
##'Bray, T. A., Crawford, G. B., and Proschan, F. (1967). \emph{Maximum
##'Likelihood Estimation of a U-shaped Failure Rate Function}. Defense
##'Technical Information Center.
##'
##'Jankowski, H. K. and Wellner, J. A. (2009a). Computation of nonparametric
##'convex hazard estimators via profile methods. \emph{Journal of Nonparametric
##'Statistics}, \bold{21}, 505-518.
##'
##'Jankowski, H. K. and Wellner, J. A. (2009b). Nonparametric estimation of a
##'convex bathtub-shaped hazard function. \emph{Bernoulli}, \bold{15},
##'1010-1035.
##'
##'Wang, Y. and Fani, S. (2018). Nonparametric maximum likelihood
##' computation of a U-shaped hazard function. \emph{Statistics and
##' Computing}, \bold{28}, 187-200.
##'
##'@keywords function
##'@examples
##'
##'## Interval-censored observations
##'data(ap)
##'(r = Uhaz(ap, deg=0))
##'plot(r, ylim=c(0,.3), col=1)
##'for(i in 1:6) plot(Uhaz(ap, deg=i/2), add=TRUE, col=i+1)
##'legend(15, 0.01, paste0("deg = ", 0:6/2), lwd=2, col=1:7, xjust=1, yjust=0)
##'
##'## Exact observations
##'data(nzmort)
##'x = with(nzmort, nzmort[ethnic=="maori",])[,1:2] # Maori mortality
##'(h0 = Uhaz(x[,1]+0.5, x[,2], deg=0)$h) # U-shaped hazard
##'(h1 = Uhaz(x[,1]+0.5, x[,2], deg=1)$h) # convex hazard
##'(h2 <- Uhaz(x[,1]+0.5, x[,2], deg=2)$h) # smooth U-shaped hazard
##'
##'plot(h0, pch=2) # plot hazard functions
##'plot(h1, add=TRUE, col="green3", pch=1)
##'plot(h2, add=TRUE, col="red3", pch=19)
##'
##'age = 0:max(x[,1]) # plot densities
##'count = integer(length(age))
##'count[x[,"age"]+1] = x[,"deaths"]
##'barplot(count/sum(count), space=0, col="lightgrey")
##'axis(1, pos=NA, at=0:10*10)
##'plot(h0, fn="d", add=TRUE, pch=2)
##'plot(h1, fn="d", add=TRUE, col="green3", pch=1)
##'plot(h2, fn="d", add=TRUE, col="red3", pch=19)
##'
##'plot(h0, fn="s", pch=2) # plot survival functions
##'plot(h1, fn="s", add=TRUE, col="green3", pch=1)
##'plot(h2, fn="s", add=TRUE, col="red3", pch=19)
##'
##'## Exact and right-censored observations
##'data(gastric)
##'plot(h0<-Uhaz(gastric, deg=0)$h) # plot hazard functions
##'plot(h1<-Uhaz(gastric, deg=1)$h, add=TRUE, col="green3")
##'plot(h2<-Uhaz(gastric, deg=2)$h, add=TRUE, col="red3")
##'
##'plot(npsurv(gastric), fn="s", col="grey") # plot survival functions
##'plot(h0, fn="s", add=TRUE)
##'plot(h1, fn="s", add=TRUE, col="green3")
##'plot(h2, fn="s", add=TRUE, col="red3")
##'
##'@export Uhaz
Uhaz = function(data, w=1, deg=1, maxit=100, tol=1e-6, verb=0) {
x = icendata(data, w)
h = uh.initial(x, deg)
attr(h, "ll") = logLikuh(h, x)
expdH = NULL
bc = TRUE # boundary change
convergence = 1
for(i in 1:maxit){
h.old = h
if(nrow(x$o) > 0) expdH = exp(chazuh(x$o[,1],h) - chazuh(x$o[,2],h))
maxima = maxgrad(h, x, expdH, bc=bc)
np1 = maxima$np1
np2 = maxima$np2
h = uh(h$alpha, c(h$tau, np1), c(h$nu, double(length(np1))),
c(h$eta, np2), c(h$mu, double(length(np2))),
h$upper, h$deg, collapse=TRUE)
r = updatemass(h, x, expdH, tol=tol)
h = r$h
if(h$deg == 0) {h = simplify(h); attr(h, "ll") = logLikuh(h, x)}
if(verb>0) { cat("##### Iteration", i, "#####\n")
cat("Log-likelihood: ", signif(attr(h,"ll"), 6), "\n")
if(verb>1) cat("Gradient values: ", signif(dlogLik(h, x), 6), "\n")
if(verb>2) {cat("hazard function:\n"); print(h)} }
if(r$convergence == 1) bc = FALSE # backtracking failed.
else if(attr(h, "ll") <= attr(h.old, "ll") + tol) {convergence = 0; break}
}
r = list(convergence=convergence, grad=dlogLik(h, x), numiter=i,
ll=attr(h, "ll"), h=h)
class(r) = "Uhaz"
r
}
# Update masses
updatemass = function(h, x, expdH=NULL, tol=1e-10) {
tau = h$tau
k = length(tau)
j2 = h$eta != h$upper
eta = h$eta = h$eta[j2]
h$mu = h$mu[j2]
m = length(eta)
p = h$deg
D1 = D2 = NULL
t1 = x$t[x$i1]
n1 = length(t1)
if(n1 > 0) {
tau.r = rep(tau, rep.int(n1,k))
dim(tau.r) = c(n1, k)
if(p > 0) tau.t = pmax(tau.r - t1, 0)
if(m > 0) {
eta.r = rep(eta, rep.int(n1,m))
dim(eta.r) = c(n1, m)
if(p > 0) t.eta = pmax(t1 - eta.r, 0)
}
D1 = switch(as.character(p),
"0" = cbind(1, tau.r >= t1, if(m>0) t1 >= eta.r else NULL) / hazuh(t1, h),
"1" = cbind(1, tau.t, if(m>0) t.eta else NULL) / hazuh(t1, h),
"2" = cbind(1, tau.t * tau.t,
if(m>0) t.eta * t.eta else NULL) / hazuh(t1, h),
cbind(1, tau.t^p, if(m>0) t.eta^p else NULL) / hazuh(t1, h) )
}
n2 = nrow(x$o)
if(n2 > 0) {
if(is.null(expdH)) expdH = exp(chazuh(x$o[,1],h) - chazuh(x$o[,2],h))
delta = sqrt(expdH) / (1 - expdH)
tau.r1 = rep(tau, rep.int(n2,k))
dim(tau.r1) = c(n2, k)
tau.x1 = pmax(tau.r1 - x$o[,1], 0)
tau.x2 = pmax(tau.r1 - x$o[,2], 0)
xd0 = x$o[,1] - x$o[,2]
xd1 = switch(as.character(p),
"0" = tau.x2 - tau.x1,
"1" = .5 * (tau.x2 * tau.x2 - tau.x1 * tau.x1),
"2" = (tau.x2 * tau.x2 * tau.x2 - tau.x1 * tau.x1 * tau.x1) / 3,
(tau.x2^(p+1) - tau.x1^(p+1)) / (p+1) )
if(m > 0) {
eta.r2 = rep(eta, rep.int(n2,m))
dim(eta.r2) = c(n2, m)
x1.eta = pmax(x$o[,1] - eta.r2, 0)
x2.eta = pmax(x$o[,2] - eta.r2, 0)
xd2 = switch(as.character(p),
"0" = x1.eta - x2.eta,
"1" = .5 * (x1.eta * x1.eta - x2.eta * x2.eta),
"2" = (x1.eta * x1.eta * x1.eta - x2.eta * x2.eta * x2.eta) / 3,
(x1.eta^(p+1) - x2.eta^(p+1)) / (p+1) )
}
else xd2 = NULL
D2 = cbind(xd0, xd1, xd2) * delta
D2[delta == 0] = 0
}
D = rbind(D1, D2) * sqrt(c(x$wt[x$i1], x$wo))
H = crossprod(D) # Choleski decomposition
v = sqrt(diag(H))
jv = v != 0
Hv = H[jv,jv,drop=FALSE] / tcrossprod(v[jv])
diag(Hv) = diag(Hv) + 1e-10
Rv = chol(Hv)
gv = dlogLik(h, x, expdH, interior=TRUE)[jv] / v[jv]
plus = forwardsolve(Rv, gv, upper.tri=TRUE, transpose=TRUE)
par = c(h$alpha, h$nu, h$mu[j2])[jv] * v[jv]
w.new = double(length(v))
w.new[jv] = nnls(Rv, Rv %*% par + plus)$x / v[jv]
alpha = w.new[1]
nu = if(k > 0) w.new[2:(k+1)] else numeric()
mu = if(m > 0) w.new[(k+2):(k+m+1)] else numeric()
newh = uh(alpha=alpha, tau=h$tau, nu=nu, eta=h$eta, mu=mu,
upper=x$upper, h$deg, collapse=FALSE)
if(h$deg == 0) b = backtrack(h, newh, x, expdH, alpha=0)
else b = backtrack(h, newh, x, expdH)
newh = b$h2
j1 = newh$nu != 0
j2 = newh$mu != 0
h2 = uh(newh$alpha, newh$tau[j1], newh$nu[j1], newh$eta[j2], newh$mu[j2],
upper=h$upper, h$deg, collapse=FALSE)
h = collapse(h2, x, tol=pmax(tol,1e-10))
list(h=h, convergence=b$convergence)
}
# Backtracking line search. h and h2 must have the same knots
backtrack = function(h, h2, x, expdH, tol=1e-10, alpha=0.33){
j = h$eta != h$upper
h2$eta = h$eta = h$eta[j]
h$mu = h$mu[j]
h2$mu = h2$mu[j]
ll.h = logLikuh(h, x)
d = c(h2$alpha - h$alpha, h2$nu - h$nu, h2$mu - h$mu)
g = alpha * sum(dlogLik(h, x, expdH) * d)
convergence = 0
r = 1
repeat {
hr = uh((1-r) * h$alpha + r * h2$alpha,
h$tau, (1-r) * h$nu + r * h2$nu,
h$eta, (1-r)*h$mu + r * h2$mu, upper=h$upper, deg=h$deg)
ll.hr = logLikuh(hr, x)
if (ll.hr >= ll.h + r * g) {convergence =0; break}
r = 0.5 * r
if(r < tol) {r = 0; hr = h; ll.h2 = ll.h; convergence = 1; break}
}
attr(h2, "ll") = ll.hr
list(h2=hr, r=r, convergence=convergence)
}
# Zeroth gradient function
grad0 = function(h, x, expdH=NULL) {
if(length(x$t) > 0)
d0 = sum(x$wt[x$i1] / hazuh(x$t[x$i1], h)) - sum(x$wt * x$t)
else d0 = 0
if(nrow(x$o) > 0) {
if(is.null(expdH)) expdH = exp(chazuh(x$o[,1],h) - chazuh(x$o[,2],h))
Delta = expdH / (1 - expdH)
xd = x$o[,1] - x$o[,2]
xd[Delta == 0] = 0
d0 = d0 - sum(x$wo * (x$o[,1] + xd * Delta))
}
d0
}
# First gradient function
grad1 = function(tau, h, x, expdH=NULL, order=0) {
g = vector("list", length(order))
names(g) = paste("d", order, sep="")
g[1:length(g)] = 0
if(length(tau) == 0) return(NULL)
m = length(tau)
n1 = length(x$t)
p = h$deg
if(n1 > 0) { # for exact observations
tau.r1 = rep(tau, rep.int(n1,m))
dim(tau.r1) = c(n1,m)
ind = tau.r1 >= x$t
tau.t = pmax(tau.r1 - x$t, 0)
if(0 %in% order) {
g$d0 = switch(as.character(p),
"0" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
ind[x$i1,,drop=FALSE]) - crossprod(x$wt, tau.r1 - tau.t))[1,],
"1" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
tau.t[x$i1,,drop=FALSE]) -
0.5 * crossprod(x$wt, tau.r1 * tau.r1 - tau.t * tau.t))[1,],
"2" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
tau.t[x$i1,,drop=FALSE]^2) -
crossprod(x$wt, tau.r1*tau.r1*tau.r1 - tau.t*tau.t*tau.t) / 3)[1,],
(crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
tau.t[x$i1,,drop=FALSE]^p) -
crossprod(x$wt, tau.r1^(p+1) - tau.t^(p+1))[1,] / (p+1))[1,] )
}
if(1 %in% order) {
g$d1 = switch(as.character(p),
"0" = double(m),
"1" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
ind[x$i1,,drop=FALSE]) - crossprod(x$wt, tau.r1 - tau.t))[1,],
"2" = (2 * crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
tau.t[x$i1,,drop=FALSE]) -
crossprod(x$wt, tau.r1 * tau.r1 - tau.t * tau.t))[1,],
{ tau.t.pm1 = tau.t[x$i1,,drop=FALSE]^(p-1)
if(p < 1) tau.t.pm1[tau.t.pm1 == Inf] = 0
(p * crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h), tau.t.pm1) -
crossprod(x$wt, tau.r1^p - tau.t^p))[1,]
} )
}
if(2 %in% order) {
g$d2 = switch(as.character(p),
"0" = double(m),
"1" = - sum(x$wt) + crossprod(x$wt, ind)[1,],
"2" = (2 * crossprod(x$wt[x$i1], ind[x$i1,,drop=FALSE]) -
crossprod(x$wt, tau.r1 - tau.t))[1,],
{ tau.t1.pm2 = tau.t[x$i1,,drop=FALSE]^(p-2)
if(p < 2) tau.t1.pm2[tau.t1.pm2 == Inf] = 0
tau.t.pm1 = tau.t^(p-1)
if(p < 1) tau.t.pm1[tau.t.pm1 == Inf] = 0
(p * (p-1) * crossprod(x$wt[x$i1], tau.t1.pm2) -
p * crossprod(x$wt, tau.r1^(p-1) - tau.t.pm1))[1,]
} )
}
}
n2 = nrow(x$o)
if(n2 > 0) { # for interval-censored observations
if(is.null(expdH)) expdH = exp(chazuh(x$o[,1],h) - chazuh(x$o[,2],h))
Delta = expdH / (1 - expdH)
tau.r2 = rep(tau, rep.int(n2,m))
dim(tau.r2) = c(n2,m)
tau.x1 = pmax(tau.r2 - x$o[,1], 0)
tau.x2 = pmax(tau.r2 - x$o[,2], 0)
ind1 = tau.r2 >= x$o[,1]
ind2 = tau.r2 >= x$o[,2]
if(0 %in% order) {
xd0 = switch(as.character(p),
"0" = {tau.r2.p1 = tau.r2; tau.x2 - (tau.x1.p1 <- tau.x1)},
"1" = {tau.r2.p1 = tau.r2 * tau.r2
tau.x2 * tau.x2 - (tau.x1.p1 <- tau.x1 * tau.x1)},
"2" = {tau.r2.p1 = tau.r2 * tau.r2 * tau.r2
tau.x2 * tau.x2 * tau.x2 -
(tau.x1.p1 <- tau.x1 * tau.x1 * tau.x1)},
{ tau.r2.p1 = tau.r2^(p+1)
tau.x2^(p+1) - (tau.x1.p1 <- tau.x1^(p+1))} )
xd0[Delta == 0] = 0
g$d0 = g$d0 -
crossprod(x$wo, (tau.r2.p1 - tau.x1.p1 + xd0 * Delta))[1,] / (p+1)
}
if(1 %in% order) {
xd1 = switch(as.character(p),
"0" = {tau.r2.p = 1; ind2 - (tau.x1.p <- ind1)},
"1" = {tau.r2.p = tau.r2; tau.x2 - (tau.x1.p <- tau.x1)},
"2" = {tau.r2.p = tau.r2 * tau.r2;
tau.x2 * tau.x2 - (tau.x1.p <- tau.x1 * tau.x1)},
{tau.r2.p = tau.r2^p; tau.x2^p - (tau.x1.p <- tau.x1^p)} )
xd1[Delta == 0] = 0
g$d1 = g$d1 - crossprod(x$wo, tau.r2.p - tau.x1.p + xd1 * Delta)[1,]
}
if(2 %in% order) {
g$d2 = switch(as.character(p),
"0" = double(m),
"1" = g$d2 - sum(x$wo) +
crossprod(x$wo, ind1 - (ind2 - ind1) * Delta)[1,],
"2" = { xd2 = tau.x2 - tau.x1
xd2[Delta == 0] = 0
g$d2 -
2 * crossprod(x$wo, tau.r2 - tau.x1 + xd2 * Delta)[1,] },
{ tau.x1.pm1 = tau.x1^(p-1)
if(p < 1) tau.x1.pm1[tau.x1.pm1 == Inf] = 0
xdp = tau.x2^(p-1) - tau.x1^(p-1)
xdp[Delta == 0] = 0
g$d2 - p * crossprod(x$wo,
tau.r2^(p-1) - tau.x1^(p-1) + xdp * Delta)[1,]
} )
}
}
g
}
# Second gradient function
grad2 = function(eta, h, x, expdH=NULL, order=0) {
g = vector("list", length(order))
names(g) = paste("d", order, sep="")
g[1:length(g)] = 0
if(length(eta) == 0) return(NULL)
m = length(eta)
n1 = length(x$t)
p = h$deg
if(n1 > 0) {
eta.r1 = rep(eta, rep.int(n1,m))
dim(eta.r1) = c(n1,m)
t.eta = pmax(x$t - eta.r1, 0)
ind = x$t >= eta.r1
if(0 %in% order) {
g$d0 = switch(as.character(p),
"0" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
ind[x$i1,,drop=FALSE]) - crossprod(x$wt, t.eta))[1,],
"1" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
t.eta[x$i1,,drop=FALSE]) - 0.5 * crossprod(x$wt, t.eta^2))[1,],
"2" = (crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
t.eta[x$i1,,drop=FALSE]^2) -
crossprod(x$wt, t.eta * t.eta * t.eta) / 3)[1,],
(crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
t.eta[x$i1,,drop=FALSE]^p) -
crossprod(x$wt, t.eta^(p+1)) / (p+1))[1,] )
}
if(1 %in% order) {
g$d1 = switch(as.character(p),
"0" = double(m),
"1" = (-crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
ind[x$i1,,drop=FALSE]) +
crossprod(x$wt, t.eta))[1,],
"2" = (-2 * crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h),
t.eta[x$i1,,drop=FALSE]) +
crossprod(x$wt, t.eta^2))[1,],
{ t.eta.pm1 = t.eta[x$i1,,drop=FALSE]^(p-1)
if(p < 1) t.eta.pm1[t.eta.pm1 == Inf] = 0
(-p * crossprod(x$wt[x$i1] / hazuh(x$t[x$i1], h), t.eta.pm1) +
crossprod(x$wt, t.eta^p))[1,]
} )
}
if(2 %in% order) {
g$d2 = switch(as.character(p),
"0" = double(m),
"1" = - crossprod(x$wt, ind)[1,],
"2" = - 2 * crossprod(x$wt, t.eta)[1,],
{ t.eta.pm1 = t.eta^(p-1)
if(p < 1) t.eta.pm1[t.eta.pm1 == Inf] = 0
- p * crossprod(x$wt, t.eta.pm1)[1,]
} )
}
}
n2 = nrow(x$o)
if(n2 > 0) {
if(is.null(expdH)) expdH = exp(chazuh(x$o[,1],h) - chazuh(x$o[,2],h))
Delta = expdH / (1 - expdH)
eta.r2 = rep(eta, rep.int(n2,m))
dim(eta.r2) = c(n2,m)
x1.eta = pmax(x$o[,1] - eta.r2, 0)
x2.eta = pmax(x$o[,2] - eta.r2, 0)
ind1 = x$o[,1] >= eta.r2
ind2 = x$o[,2] >= eta.r2
if(0 %in% order) {
xd0 = switch(as.character(p),
"0" = (x1.eta.p1 <- x1.eta) - x2.eta,
"1" = (x1.eta.p1 <- x1.eta * x1.eta) - x2.eta * x2.eta,
"2" = (x1.eta.p1 <- x1.eta * x1.eta * x1.eta) -
x2.eta * x2.eta * x2.eta,
(x1.eta.p1 <- x1.eta^(p+1)) - x2.eta^(p+1) )
xd0[Delta == 0] = 0
g$d0 = g$d0 - crossprod(x$wo, x1.eta.p1 + xd0 * Delta)[1,] / (p+1)
}
if(1 %in% order) {
xd1 = switch(as.character(p),
"0" = (x1.eta.p <- ind1) - ind2,
"1" = (x1.eta.p <- x1.eta) - x2.eta,
"2" = (x1.eta.p <- x1.eta * x1.eta) - x2.eta * x2.eta,
(x1.eta.p <- x1.eta^p) - x2.eta^p )
xd1[Delta == 0] = 0
g$d1 = g$d1 + crossprod(x$wo, x1.eta.p + xd1 * Delta)[1,]
}
if(2 %in% order) {
g$d2 = switch(as.character(p),
"0" = double(m),
"1" = g$d2 - crossprod(x$wo, ind1 + (ind1 - ind2) * Delta)[1,],
"2" = { xd2 = x1.eta - x2.eta
xd2[Delta == 0] = 0
g$d2 - 2 * crossprod(x$wo, x1.eta + xd2 * Delta)[1,] },
{ x1.eta.pm1 = x1.eta^(p-1)
x2.eta.pm1 = x2.eta^(p-1)
if(p < 1) {
x1.eta.pm1[x1.eta.pm1 == Inf] = 0
x2.eta.pm1[x2.eta.pm1 == Inf] = 0
}
xdp = x1.eta.pm1 - x2.eta.pm1
xdp[Delta == 0] = 0
g$d2 - p * crossprod(x$wo, x1.eta.pm1 + xdp * Delta)[1,]
} )
}
}
g
}
##'Computes the Log-likelihood Value of a U-shaped Hazard Function
##'
##'
##'\code{logLikuh} returns the log-likelihood value of a U-shaped hazard
##'function, given a data set.
##'
##'
##'@param h an object of class \code{uh}.
##'@param data numeric vector or matrix for exact or interval-censored
##'observations, or an object of class \code{icendata}.
##'@return
##'
##'Log-likelihood value evaluated at \code{h}, given \code{data}.
##'@author Yong Wang <yongwang@@auckland.ac.nz>
##'@seealso \code{\link{Uhaz}}, \code{\link{icendata}}, \code{\link{plot.uh}}
##'@references
##'
##'Wang, Y. and Fani, S. (2018). Nonparametric maximum likelihood
##' computation of a U-shaped hazard function. \emph{Statistics and
##' Computing}, \bold{28}, 187-200.
##'
##'@keywords function
##'@examples
##'
##'data(ap)
##'(h0 = uh(.2, NULL, NULL, NULL, NULL, 15, 1)) # Uniform hazard
##'plot(h0, ylim=c(0,.3))
##'logLikuh(h0, ap)
##'
##'r = Uhaz(ap, deg=2)
##'r$ll
##'logLikuh(r$h, ap)
##'plot(r$h, add=TRUE, col="red3")
##'
##'@export logLikuh
logLikuh = function(h, data) {
x = icendata(data)
if(length(x$t) > 0) ll = sum(x$wt[x$i1] * log(hazuh(x$t[x$i1], h))) -
sum(x$wt * chazuh(x$t, h))
else ll = 0
if(nrow(x$o) > 0) {
ll = ll + sum(x$wo *
(log(exp(-chazuh(x$o[,1], h)) - exp(-chazuh(x$o[,2], h)))))
}
ll
}
dlogLik = function(h, x, expdH=NULL, interior=FALSE) {
if(interior) eta = h$eta[h$eta != h$upper]
else eta = h$eta
m = length(eta)
d = c(grad0(h, x, expdH), grad1(h$tau, h, x, expdH)$d0,
if(m>0) grad2(eta, h, x, expdH)$d0 else NULL)
names(d) = c("alpha", paste("nu",1:length(h$tau),sep=""),
if(m>0) paste("mu",1:m,sep="") else NULL)
d
}
### Finds all local maxima of the gradient functions
# bc boundary change
maxgrad = function(h, x, expdH=NULL, maxit=100, grid=100, bc=TRUE) {
if(!is.icendata(x)) x = icendata(x, w=1)
u = sort(unique(c(x$u, h$tau, h$eta)))
p = h$deg
if(p > 0 && p < 0.1)
u = rep(u, rep.int(21, length(u))) + seq(-h$upper*1e-2, h$upper*1e-2, len=21)
if(p == 1) maxit = 1
tau1 = c(0, h$tau[h$tau > 0])
k = length(tau1) - 1
eta1 = c(h$eta[h$eta<h$upper], h$upper)
m = length(eta1) - 1
if(p < 1) grid = grid * log2(2/(p*p))
## grad1
if(p < 0.1) { # optimal points
if(bc) pt1 = u[u<=eta1[1]]
else pt1 = u[u<eta1[1]]
gpt1 = grad1(pt1, h, x, expdH)$d0
}
else {
if(p == 1) u1 = u[u <= eta1[1]]
else u1 = seq(0, eta1[1], len=grid)
if(!bc) u1 = u1[u1 < eta1[1]]
m1 = length(u1)
pt1 = gpt1 = numeric()
if(p == 1) jd = rep(TRUE, m1-1)
else {
g = grad1(u1, h, x, expdH, order=1)
jd = g$d1[-m1] > 0 & g$d1[-1] < 0
}
if(any(jd)) {
left = u1[-m1][jd]
right = u1[-1][jd]
pt1 = (left + right) * .5
for(i in 1:maxit) {
g = grad1(pt1, h, x, expdH, order=1:2)
left[g$d1>0] = pt1[g$d1>0]
right[g$d1<0] = pt1[g$d1<0]
pt1.old = pt1
pt1 = pt1 - g$d1 / g$d2
j = is.na(pt1) | pt1 < left | pt1 > right
pt1[j] = (left[j] + right[j]) * .5
if( max(abs(pt1 - pt1.old)) <= 1e-14 * h$upper) break
}
if(p == 1) pt1 = pt1[!j]
gpt1 = grad1(pt1, h, x, expdH, order=0)$d0
}
}
i = pt1 > 0 & pt1 <= tau1[k+1]
pt1i = pt1[i]
gpt1i = gpt1[i]
if(k > 0 && length(pt1i) > 0) {
grp = apply(outer(pt1i, tau1[-(k+1)], ">") & outer(pt1i, tau1[-1], "<="),
1, which.max)
r1 = aggregate(gpt1i, by=list(group=grp), which.max)
j = integer(k)
j[r1[,1]] = r1[,2]
j = j + c(0,cumsum(tabulate(grp, nbins=k))[-k])
np1 = pt1i[j]
gnp1 = gpt1i[j]
j0 = gnp1 > 0
np1 = np1[j0]
gnp1 = gnp1[j0]
}
else np1 = gnp1 = numeric()
## grad2
if(p < 0.1) {
if(bc) pt2 = u[u>=tau1[k+1] & u<=h$upper]
else pt2 = u[u>tau1[k+1] & u<h$upper]
gpt2 = grad2(pt2, h, x, expdH)$d0
}
else {
if(p == 1) u2 = u[u >= tau1[k+1]]
else u2 = seq(tau1[k+1], h$upper, len=grid)
if(!bc) u2 = u2[u2 > tau1[k+1]]
m2 = length(u2)
pt2 = gpt2 = numeric()
if(p == 1) jd = rep(TRUE, m2-1)
else {
g = grad2(u2, h, x, expdH, order=1)
jd = g$d1[-m2] > 0 & g$d1[-1] < 0
}
if(any(jd)) {
left = u2[-m2][jd]
right = u2[-1][jd]
pt2 = (left + right) * .5
for(i in 1:maxit) {
g = grad2(pt2, h, x, expdH, order=1:2)
left[g$d1>0] = pt2[g$d1>0]
right[g$d1<0] = pt2[g$d1<0]
pt2.old = pt2
pt2 = pt2 - g$d1 / g$d2
j = is.na(pt2) | pt2 < left | pt2 > right
pt2[j] = (left[j] + right[j]) * .5
if( max(abs(pt2 - pt2.old)) <= 1e-14 * h$upper) break
}
if(p == 1) pt2 = pt2[!j]
gpt2 = grad2(pt2, h, x, expdH, order=0)$d0
}
}
i = pt2 >= eta1[1] & pt2 < eta1[m+1]
pt2i = pt2[i]
gpt2i = gpt2[i]
if(m > 0 && length(pt2i) > 0) {
grp = apply(outer(pt2i, eta1[-(m+1)], ">=") & outer(pt2i, eta1[-1], "<"), 1,
which.max)
r2 = aggregate(gpt2i, by=list(group=grp), which.max)
j = integer(m)
j[r2[,1]] = r2[,2]
j = j + c(0, cumsum(tabulate(grp, nbins=m))[-m])
np2 = pt2i[j]
gnp2 = gpt2i[j]
j0 = gnp2 > 0
np2 = np2[j0]
gnp2 = gnp2[j0]
}
else np2 = gnp2 = numeric()
## grad1 and grad2
if(max(h$tau) != h$eta[1]) {
jj1 = pt1 >= tau1[k+1] & pt1 <= eta1[1]
if(p == 0) { uj1 = pt1[jj1]; gj1 = gpt1[jj1] }
else {
uj1 = c(tau1[k+1], if(bc) eta1[1] else NULL, pt1[jj1])
gj1 = c(grad1(c(tau1[k+1], if(bc) eta1[1] else NULL),
h, x, expdH)$d0, gpt1[jj1])
}
jmax = which.max(gj1)
np31 = uj1[jmax]
gnp31 = gj1[jmax]
jj2 = pt2 >= tau1[k+1] & pt2 <= eta1[1]
if(p == 0) { uj2 = pt2[jj2]; gj2 = gpt2[jj2] }
else {
uj2 = c(if(bc) tau1[k+1] else NULL, eta1[1], pt2[jj2])
gj2 = c(grad2(c(if(bc) tau1[k+1] else NULL, eta1[1]),
h, x, expdH)$d0, gpt2[jj2])
}
jmax = which.max(gj2)
np32 = uj2[jmax]
gnp32 = gj2[jmax]
if(gnp31 > gnp32) {np1 = c(np1, np31); gnp1 = c(gnp1, gnp31)}
else {np2 = c(np2, np32); gnp2 = c(gnp2, gnp32)}
}
list(np1=np1, gnp1=gnp1, np2=np2, gnp2=gnp2)
}
simplify = function(h) {
i1 = order(h$tau) # remove identical knots
tau = h$tau[i1]
nu = h$nu[i1]
i2 = order(h$eta)
eta = h$eta[i2]
mu = h$mu[i2]
if(h$deg != 0 || length(tau) == 0 || length(eta) == 0) return(h)
if(tau[length(tau)] != eta[1]) return(h)
if(nu[length(tau)] < mu[1]) {
if(eta[1] == 0) {
h$alpha = h$alpha + mu[1]
eta = eta[-1]
mu = mu[-1]
}
else {
h$alpha = h$alpha + nu[length(nu)]
mu[1] = mu[1] - nu[length(nu)]
tau = tau[-length(tau)]
nu = nu[-length(nu)]
}
}
else {
h$alpha = h$alpha + mu[1]
nu[length(nu)] = nu[length(nu)] - mu[1]
eta = eta[-1]
mu = mu[-1]
}
uh(h$alpha, tau, nu, eta, mu, h$upper, h$deg, collapse=FALSE)
}
# Collapse similar knots
collapse = function(h, x, tol=0){
ll = attr(h, "ll")
i1 = order(h$tau) # remove identical knots
tau = h$tau[i1]
ind1 = cumsum(!duplicated(c(tau)))
tau = unique(tau)
nu = aggregate(h$nu[i1], by=list(group=ind1), sum)[,2]
nu[tau == 0] = 0
i2 = order(h$eta)
eta = h$eta[i2]
ind2 = cumsum(!duplicated(eta))
eta = unique(eta)
mu = aggregate(h$mu[i2], by=list(group=ind2), sum)[,2]
mu[eta == h$upper] = 0
h = uh(h$alpha, tau, nu, eta, mu, h$upper, h$deg, collapse=FALSE)
if(tol > 0) {
if(is.null(ll)) ll = logLikuh(h, x)
# if(h$deg < 1) {attr(h, "ll") = ll; return(h)} # why?
h2 = h
break1 = break2 = FALSE
repeat {
if(!break1 && length(h2$nu) > 1) {
j = which.min(diff(h$tau))
h2$nu[j] = h2$nu[j] + h2$nu[j+1]
h2$nu = h2$nu[-(j+1)]
h2$tau[j] = (h2$tau[j] + h2$tau[j+1]) * .5
h2$tau = h2$tau[-(j+1)]
ll2 = logLikuh(h2, x)
if(ll2 + tol >= ll) {h = h2; ll = ll2; break1 = FALSE}
else {h2 = h; break1 = TRUE}
}
else break1 = TRUE
if(!break2 && length(h2$mu) > 1) {
j = which.min(diff(h2$eta))
h2$mu[j] = h2$mu[j] + h2$mu[j+1]
h2$mu = h2$mu[-(j+1)]
h2$eta[j] = (h2$eta[j] + h2$eta[j+1]) * .5
h2$eta = h2$eta[-(j+1)]
ll2 = logLikuh(h2, x)
if(ll2 + tol >= ll) {h = h2; ll = ll2; break2 = FALSE}
else {h2 = h; break2 = TRUE}
}
else break2 = TRUE
if(break1 && break2) break
}
attr(h, "ll") = ll
}
h
}
# deg - polynomial degree
##'U-shaped Hazard Function
##'
##'
##' Class \code{uh} can be used to store U-shaped hazard functions.
##' There are a couple of functions associated with the class.
##'
##'\code{uh} creates an object of class \code{uh}, which stores a U-shaped
##'hazard function.
##'
##'\code{print.uh} prints an object of class \code{uh}.
##'
##'
##'A U-shape hazard function, as generalized by Wang and Fani (2018), is given
##'by
##'
##'\deqn{h(t) = \alpha + \sum_{j = 1}^k \nu_j(\tau_j - t)_+^p + \sum_{j = 1}^{m} \mu_j (t-\eta_j)_+^p,}{ h(t) = alpha + sum_{j=1}^k nu_j (tau_j - t)_+^p
##' + sum_{j=1}^m mu_j (t - eta_j)_+^p,}
##'
##'where \eqn{\alpha,\nu_j,\mu_j \ge 0}{alpha, nu_j, mu_j \ge 0},
##'\eqn{\tau_1 < \cdots < \tau_k \le \eta_1 < \cdots < \eta_m,}{tau_1 < ... < tau_k <= eta_1 < ... < eta_m,} and \eqn{p \ge 0}{p >= 0} is the the spline degree which
##'determines the smoothness of the U-shaped hazard. As p increases, the family
##'of hazard functions becomes increasingly smoother, but at the same time,
##'smaller. When \eqn{p = 0}{p = 0}, the hazard function is U-shaped, as
##'studied by Bray et al. (1967). When \eqn{p = 1}{p = 1}, the hazard function
##'is convex, as studied by Jankowski and Wellner (2009a,b).
##'
##'\code{print.uh} prints an object of class \code{uh}. While \code{alpha},
##'\code{upper} and \code{deg} are printed as they are, \code{tau} and
##'\code{nu} are printed as a two-column matrix, and so are \code{eta} and
##'\code{mu}.
##'
##'@aliases uh uh.object print.uh
##'@param alpha a nonnegative value, for the constant coefficient.
##'@param tau vector of nonnegative real values, for left knots.
##'@param nu vector of nonnegative values, for masses associated with the left
##'knots.
##'@param eta vector of nonnegative real values, for right knots.
##'@param mu vector of nonnegative real values, for masses associated with the
##'right knots.
##'@param upper a positive value, at which point the hazard starts to become
##'infinite.
##'@param deg nonnegative real number for spline degree (i.e., p in the formula
##'below).
##'@param collapse logical, indicating if identical knots should be collapsed.
##'@param x an object of class \code{uh}.
##'@param ... other auguments for printing.
##'@return
##'
##'\code{uh} returns an object of class \code{uh}. It is a list with components
##'\code{alpha}, \code{tau}, \code{nu}, \code{eta}, \code{mu}, \code{upper} and
##'\code{deg}, which store their corresponding values as described above.
##'@author Yong Wang <yongwang@@auckland.ac.nz>
##'@seealso \code{\link{Uhaz}}, \code{\link{icendata}}, \code{\link{plot.uh}}
##'@references
##'
##'Bray, T. A., Crawford, G. B., and Proschan, F. (1967). \emph{Maximum
##'Likelihood Estimation of a U-shaped Failure Rate Function}. Defense
##'Technical Information Center.
##'
##'Jankowski, H. K. and Wellner, J. A. (2009a). Computation of nonparametric
##'convex hazard estimators via profile methods. \emph{Journal of Nonparametric
##'Statistics}, \bold{21}, 505-518.
##'
##'Jankowski, H. K. and Wellner, J. A. (2009b). Nonparametric estimation of a
##'convex bathtub-shaped hazard function. \emph{Bernoulli}, \bold{15},
##'1010-1035.
##'
##'Wang, Y. and Fani, S. (2018). Nonparametric maximum likelihood
##' computation of a U-shaped hazard function. \emph{Statistics and
##' Computing}, \bold{28}, 187-200.
##'
##'@keywords function
##'@examples
##'
##'(h0 = uh(3, 2, 3, 4, 5, 7, deg=0)) # deg = 0
##'plot(h0, ylim=c(0,20))
##'(h1 = uh(4, 2, 3, 5, 6, 7, deg=1)) # deg = 1
##'plot(h1, add=TRUE, col="green3")
##'(h2 = uh(1, 1:2, 3:4, 5:6, 7:8, 9, deg=2)) # deg = 2
##'plot(h2, add=TRUE, col="red3")
##'
##'@usage
##'uh(alpha, tau, nu, eta, mu, upper=Inf, deg=1, collapse=TRUE)
##'\method{print}{uh}(x, ...)
##'
##'@export uh
##'@export print.uh
uh = function(alpha, tau, nu, eta, mu, upper=Inf, deg=1, collapse=TRUE) {
if(length(tau) == 0) {tau=0; nu=0}
if(length(eta) == 0) {eta=upper; mu=0}
i1 = order(tau)
tau = tau[i1]
nu = nu[i1]
i2 = order(eta)
eta = eta[i2]
mu = mu[i2]
h = list(alpha=alpha, tau=tau, nu=nu, eta=eta, mu=mu, upper=upper, deg=deg)
if(collapse) h = collapse(h)
class(h) = "uh"
h
}
print.uh = function(x, ...) {
cat("$alpha\n")
print(x$alpha, ...)
print(cbind(tau=x$tau, nu=x$nu), ...)
print(cbind(eta=x$eta, mu=x$mu), ...)
cat("$upper\n")
print(x$upper, ...)
cat("$deg\n")
print(x$deg, ...)
}
# Hazard function
##'Distributional Functions given a U-shaped Hazard Function
##'
##'
##'Given an object of class \code{uh}:
##'
##'\code{hazuh} computes the hazard values;
##'
##'\code{chazuh} computes the cumulative hazard values;
##'
##'\code{survuh} computes the survival function values;
##'
##'\code{denuh} computes the density function values.
##'
##'
##'@aliases hazuh chazuh survuh denuh
##'@param t time points at which the function is to be evaluated.
##'@param h an object of class \code{uh}.
##'@return
##'
##'A numeric vector of the function values.
##'@author Yong Wang <yongwang@@auckland.ac.nz>
##'@seealso \code{\link{Uhaz}}, \code{\link{icendata}}, \code{\link{plot.uh}}
##'@references
##'
##'Wang, Y. and Fani, S. (2018). Nonparametric maximum likelihood
##' computation of a U-shaped hazard function. \emph{Statistics and
##' Computing}, \bold{28}, 187-200.
##'
##'@keywords function
##'@examples
##'
##'data(ap)
##'h = Uhaz(icendata(ap), deg=2)$h
##'hazuh(0:15, h) # hazard
##'chazuh(0:15, h) # cumulative hazard
##'survuh(0:15, h) # survival probability
##'denuh(0:15, h) # density
##'
##'@usage
##'hazuh(t, h)
##'chazuh(t, h)
##'survuh(t, h)
##'denuh(t, h)
##'
##'@export hazuh
##'@export chazuh
##'@export survuh
##'@export denuh
hazuh = function(t, h) {
p = h$deg
b = c = 0
if(length(h$tau) > 0) {
if(p == 0) d = outer(h$tau, t, ">=")
else {
d = pmax(outer(h$tau, t, "-"), 0)
if(p != 1) d = d^p
}
b = (h$nu %*% d)[1,]
}
if(length(h$eta) > 0) {
if(p == 0) d = outer(t, h$eta, ">=")
else {
d = pmax(outer(t, h$eta, "-"), 0)
if(p != 1) d = d^p
}
c = (d %*% h$mu)[,1]
}
h$alpha + pmax(b, c)
}
chazuh = function(t, h) {
deg = pmax(h$deg, 1)
p1 = h$deg + 1
a = b = c = 0
if(h$alpha > 0) a = h$alpha * t
if(length(h$tau) > 0) {
tau.t = pmax(outer(h$tau, t, "-"), 0)
b = (h$nu %*% (h$tau^p1 - tau.t^p1))[1,] / p1
}
if(length(h$eta) > 0) {
t.eta = pmax(outer(t, h$eta, "-"), 0)
c = (t.eta^p1 %*% h$mu)[,1] / p1
}
H = a + b + c
H[t > h$upper] = 1e100
H
}
# survival function
survuh = function(t, h) exp(-chazuh(t, h))
# density function
denuh = function(t, h) hazuh(t, h) * survuh(t, h)
## plotting functions
##'Plot Functions for U-shaped Hazard Estimation
##'
##'
##' Functions for plotting various functions in U-shaped hazard estimation
##'
##'\code{plot.Uhaz} and \code{plot.uh} are wrapper functions that can be used
##'to invoke \code{plot.hazuh}, \code{plot.chazuh}, \code{plot.survuh},
##'\code{plot.denuh} or \code{plot.graduh}.
##'
##'\code{plothazuh} plots a U-shaped hazard function.
##'
##'\code{plotchazuh} plots a cumulative hazard function that has a U-shaped
##'hazard function.
##'
##'\code{plotsurvuh} plots the survival function that has a U-shaped hazard
##'function.
##'
##'\code{plotdenuh} plots the density function that has a U-shaped hazard
##'function.
##'
##'\code{plotgraduh} plots the gradient function that has a U-shaped hazard
##'function.
##'
##'
##'A U-shaped hazard function is given by
##'
##'\deqn{h(t) = \alpha + \sum_{j = 1}^k \nu_j(\tau_j - t)_+^p + \sum_{j = 1}^{m} \mu_j (t-\eta_j)_+^p,}{ h(t) = alpha + sum_{j=1}^k nu_j (tau_j - t)_+^p
##' + sum_{j=1}^m mu_j (t - eta_j)_+^p,}
##'
##'where \eqn{\alpha,\nu_j,\mu_j \ge 0}{alpha, nu_j, mu_j \ge 0},
##'\eqn{\tau_1 < \cdots < \tau_k \le \eta_1 < \cdots < \eta_m,}{tau_1 < ... < tau_k <= eta_1 < ... < eta_m,} and \eqn{p \ge 0}{p >= 0}.
##'
##'@aliases plot.Uhaz plot.uh plothazuh plotchazuh plotsurvuh plotdenuh
##'plotgraduh
##'@param x an object of class \code{Uhaz}, i.e., an output of function
##'\code{Uhaz}, or an object of class \code{uh}..
##'@param h an object of class \code{uh}.
##'@param data vector or matrix that stores observations, or an object of class
##'\code{icendata}.
##'@param w additional weights/multiplicities for the observations stored in
##'\code{data}.
##'@param fn function to be plotted. It can be
##'
##'= \code{haz}, for hazard function;
##'
##'= \code{chaz}, for cumulative hazard function;
##'
##'= \code{den}, for density function;
##'
##'= \code{surv}, for survival function;
##'
##'= \code{gradient}, for gradient functions.
##'
##'@param xlim,ylim numeric vectors of length 2, giving the x and y coordinates
##'ranges.
##'@param xlab,ylab x- or y-axis labels.
##'@param add = \code{TRUE}, adds the curve to the existing plot;
##'
##'= \code{FALSE}, plots the curve in a new one.
##'@param col color used for plotting the curve.
##'@param lty line type for plotting the curve.
##'@param lwd line width for plotting the curve.
##'@param len number of points used to plot a curve.
##'@param add.knots logical, indicating if knots are also plotted.
##'@param pch point character/type for plotting knots.
##'@param vert logical, indicating if grey vertical lines are plotted to show
##'the interval that separates the two discrete measures.
##'@param col0 color for gradient function 0, i.e., for the hazard-constant
##'part, or alpha.
##'@param col1 color for gradient function 1, i.e., for the hazard-decreasing
##'part.
##'@param col2 color for gradient function 1, i.e., for the hazard-increasing
##'part.
##'@param order = 0, the gradient functions are plotted;
##'
##'= 1, their first derivatives are plotted;
##'
##'= 2, their second derivatives are plotted.
##'@param ... arguments for other graphical parameters (see \code{par}).
##'@author Yong Wang <yongwang@@auckland.ac.nz>
##'@seealso \code{\link{icendata}}, \code{\link{uh}}, \code{\link{npsurv}}.
##'@references
##'
##'Wang, Y. and Fani, S. (2018). Nonparametric maximum likelihood
##' computation of a U-shaped hazard function. \emph{Statistics and
##' Computing}, \bold{28}, 187-200.
##'
##'@keywords function
##'@examples
##'
##'## Angina Pectoris Survival Data
##'data(ap)
##'plot(r<-Uhaz(ap)) # hazard function for a convex hazard
##'plot(r, fn="c") # cumulative hazard function
##'plot(r, fn="s") # survival function
##'plot(r, fn="d") # density function
##'plot(r, ap, fn="g") # gradient functions
##'plot(r, ap, fn="g", order=1) # first derivatives of gradient functions
##'plot(r, ap, fn="g", order=2) # second derivatives of gradient functions
##'
##'## New Zealand Mortality in 2000
##'data(nzmort)
##'i = nzmort$ethnic == "maori"
##'x = nzmort[i,1:2] # Maori mortality
##'h = Uhaz(x[,1]+0.5, x[,2], deg=2)$h # smooth U-shaped hazard
##'plot(h) # hazard function
##'plot(h, fn="d") # density function
##'plot(h, fn="s") # survival function
##'
##'x2 = nzmort[!i,1:2] # Non-Maori mortality
##'h2 = Uhaz(x2[,1]+0.5, x2[,2], deg=2)$h
##'plot(h2, fn="s", add=TRUE, col="green3")
##'
##'@usage
##'\method{plot}{Uhaz}(x, ...)
##'\method{plot}{uh}(x, data, fn=c("haz","grad","surv","den","chaz"), ...)
##'plothazuh(h, add=FALSE, col="darkblue", lty=1, xlim, ylim,
##' lwd=2, pch=19, len=500, vert=FALSE, add.knots=TRUE,
##' xlab="Time", ylab="Hazard", ...)
##'plotchazuh(h, add=FALSE, lwd=2, len=500, col="darkblue",
##' pch=19, add.knots=TRUE, vert=FALSE, xlim, ylim, ...)
##'plotdenuh(h, add=FALSE, lty=1, lwd=2, col="darkblue",
##' add.knots=TRUE, pch=19, ylim, len=500, vert=FALSE, ...)
##'plotsurvuh(h, add=FALSE, lty=1, lwd=2, len=500, vert=FALSE,
##' col="darkblue", pch=19, add.knots=TRUE, xlim, ylim, ...)
##'plotgraduh(h, data, w=1, len=500, xlim, ylim, vert=TRUE,
##' add=FALSE, xlab="Time", ylab="Gradient",
##' col0="red3", col1="blue3", col2="green3", order=0, ...)
##'
##'@export plot.Uhaz
##'@export plot.uh
##'@export plothazuh
##'@export plotchazuh
##'@export plotsurvuh
##'@export plotdenuh
##'@export plotgraduh
plot.Uhaz = function(x, ...) plot(x$h, ...)
plot.uh = function(x, data, fn=c("haz","grad","surv","den","chaz"), ...) {
fn = match.arg(fn)
fnR = getFunction(paste("plot",fn,"uh",sep=""))
switch(fn,
"haz" =,
"surv" =,
"den" =,
"chaz" = fnR(x, ...),
"grad" = fnR(x, data, ...) )
}
plothazuh = function(h, add=FALSE, col="darkblue", lty=1,
xlim, ylim, lwd=2, pch=19, len=500, vert=FALSE, add.knots=TRUE,
xlab="Time", ylab="Hazard", ...) {
p = h$deg
pt = switch(as.character(p),
"0" = unique(sort(c(0, h$tau, h$eta, h$upper))),
"1" = unique(sort(c(0, h$tau, h$eta, h$upper))),
unique(sort(c(h$tau, h$eta, seq(0, h$upper, len=len)))))
m = length(pt)
knots = unique(c(h$tau, h$eta))
haz = hazuh(pt, h)
max.haz = max(haz)
if(missing(xlim)) xlim = range(pt)
if(missing(ylim)) ylim = c(0, max.haz)
if(!add) plot(pt, haz, type="n", xlim=xlim, ylim=ylim,
xlab=xlab, ylab=ylab, ...)
if(vert) {
lines(rep(max(h$tau),2), ylim, col="grey", lty=2)
lines(rep(min(h$eta),2), ylim, col="grey", lty=2)
}
abline(h=0, col ="black")
if(p == 0) {
lines(c(h$tau[length(h$tau)], h$eta[1]), rep(h$alpha,2),
lwd=lwd, col=col, lty=lty)
lines(c(rep(rev(h$tau),each=2), 0),
c(h$alpha, rep(hazuh(rev(h$tau), h), each=2)),
lwd=lwd, col=col, lty=lty)
lines(c(rep(h$eta,each=2), h$upper),
c(h$alpha, rep(hazuh(h$eta, h), each=2)),
lwd=lwd, col=col, lty=lty)
}
else lines(pt, haz, lwd=lwd, col=col, lty=lty)
if(add.knots && length(knots) > 0)
points(knots, hazuh(knots, h), col=col, pch=pch)
}
plotchazuh = function(h, add=FALSE, lwd=2, len=500, col="darkblue",
pch=19, add.knots=TRUE, vert=FALSE, xlim, ylim, ...) {
pt = unique(sort(c(seq(0, h$upper, len=len), h$tau, h$eta)))
m = length(pt)
H = chazuh(pt, h)
max.H = max(H)
if(missing(xlim)) xlim = range(pt)
if(missing(ylim)) ylim = c(0, max.H)
plot(rep(max(h$tau),2), c(0,max.H), type="n",
xlim=xlim, ylim=ylim, xlab="Time", ylab="Cumulative Hazard", ...)
if(vert) lines(rep(max(h$tau),2), c(0,max.H), type="l", col ="grey", lty=2)
if(vert) lines(rep(min(h$eta),2), c(0,max.H), col ="grey", lty=2)
abline(h=0, col ="black")
lines(pt, H, type="l", lwd=lwd, col=col)
if(add.knots) {
knots = c(h$tau, h$eta)
knots = knots[knots>0 & knots<h$upper]
points(knots, chazuh(knots, h), col=col, pch=pch)
}
}
plotdenuh = function(h, add=FALSE, lty=1, lwd=2,
col="darkblue", add.knots=TRUE, pch=19,
ylim, len=500, vert=FALSE, ...) {
pt = sort(c(seq(0, h$upper, len=len), h$tau, h$eta))
m = length(pt)
d = denuh(pt, h)
max.d = max(d)
type = if(vert) "l" else "n"
if(missing(ylim)) ylim = c(0, max.d)
if(add) lines(rep(max(h$tau),2), ylim, type=type, col ="grey", lty=2)
else plot(rep(max(h$tau),2), ylim, type=type,
col ="grey", lty=2, xlim=range(pt), ylim=ylim,
xlab="Time", ylab="Density", ...)
if(vert) lines(rep(min(h$eta),2), ylim, col="grey", lty=2)
abline(h=0, col ="black")
lines(pt, d, lty=lty, lwd=lwd, col=col)
if(add.knots) {
knots = c(h$tau, h$eta)
knots = knots[knots>0 & knots<h$upper]
points(knots, denuh(knots, h), col=col, pch=pch)
}
}
plotsurvuh = function(h, add=FALSE, lty=1, lwd=2, len=500, vert=FALSE,
col="darkblue", pch=19, add.knots=TRUE, xlim, ylim, ...) {
pt = sort(c(seq(0, h$upper, len=len), h$tau, h$eta))
m = length(pt)
s = survuh(pt, h)
if(missing(xlim)) xlim = range(pt)
if(missing(ylim)) ylim = c(0, 1)
if(add) {
if(vert) lines(rep(max(h$tau),2), c(0,1), col ="grey", lty=2)
}
else {
if(vert)
plot(rep(max(h$tau),2), c(0,1), type="l", col ="grey", lty=2,
xlim=xlim, ylim=ylim, xlab="Time", ylab="Survival Probability",
...)
else plot(0, 0, type="n", xlim=xlim, ylim=ylim,
xlab="Time", ylab="Survival Probability")
}
if(vert) lines(rep(min(h$eta),2), c(0,1), col ="grey", lty=2)
abline(h=0, col ="black")
lines(pt, s, lty=lty, lwd=lwd, col=col)
pt = c(h$tau, h$eta)
if(add.knots) {
knots = c(h$tau, h$eta)
knots = knots[knots>0 & knots<h$upper]
points(knots, survuh(knots, h), col=col, pch=pch)
}
}
plotgraduh = function(h, data, w=1, len=500, xlim, ylim, vert=TRUE,
add=FALSE, xlab="Time", ylab="Gradient",
col0="red3", col1="blue3", col2="green3", order=0, ...) {
if(! order %in% 0:2) stop("'order' must be 0, 1 or 2")
x = icendata(data, w)
maxf = h$upper
if(missing(xlim)) xlim = c(0, maxf)
k = length(h$tau)
tau1 = unique(sort(c(seq(xlim[1], h$eta[1] - if(h$deg!=0) 0 else 1e-4*h$upper,
len=len), h$tau)))
eta1 = unique(sort(c(seq(h$tau[k] + if(h$deg!=0) 0 else 1e-4*h$upper, xlim[2],
len=len), h$eta)))
g11 = grad1(tau1, h, x, order=order)[[1]]
g22 = grad2(eta1, h, x, order=order)[[1]]
g0 = if(order == 0) grad0(h, x) else 0
if(missing (ylim)) ylim= range(g0, g11, g22)
if(add) lines(c(0,maxf), c(0, 0), col="grey")
else plot(c(0,maxf), c(0, 0), type="l", col="grey", xlim=xlim, ylim=ylim,
xlab=xlab, ylab=ylab, ...)
if(vert) {
lines(rep(h$tau[k],2), ylim, col="gray", lty=2)
lines(rep(h$eta[1],2), ylim, col="gray", lty=2)
}
lines(eta1, g22, col=col2, lty=1, lwd=1)
lines(tau1, g11, col=col1, lty=1, lwd=1)
j1 = h$tau != 0 & h$tau != h$upper
j2 = h$eta != 0 & h$eta != h$upper
points(h$tau[j1], grad1(h$tau[j1],h,x,order=order)[[1]], col=col1, pch=2)
points(h$eta[j2], grad2(h$eta[j2],h,x,order=order)[[1]], col=col2, pch=6)
lines(c(h$tau[k], h$eta[1]), rep(g0, 2), col=col0, lty=1, lwd=1)
points((h$tau[k] + h$eta[1])/2, g0, col=col0, pch=20)
if(!add) legend("bottomright",
legend=expression(d[0], d[1], d[2]), col=c(col0, col1, col2),
lty=c(1,1,1), pch=c(20,2,6), bty="n", lwd=c(1,1,1))
}
uh.initial = function(x, deg=1) {
if(!is.icendata(x)) x = icendata(x)
n1 = length(x$t)
n2 = nrow(x$o)
y1 = if(n1 > 0) x$t else numeric()
y2 = if(n2 > 0) {
x$o[x$o[,2] == Inf,2] = 1.6 * x$upper
rowMeans(x$o)
}
else numeric()
beta = weighted.mean(c(y1,y2), c(x$wt, x$wo))
uh(alpha=1/beta, tau=NULL, nu=NULL, eta=NULL, mu=NULL, upper=x$upper,
deg=deg)
}
Try the npsurv package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.