Package cequre
performs censored quantile regression of
Huang (2010), and restores monotonicity respecting via adaptive
interpolation for dynamic regression of Huang (2017). The
monotonicity-respecting restoration applies to general dynamic
regression models including (uncensored or censored) quantile regression
model, additive hazards model, and dynamic survival models of Peng and
Huang (2007), among others.
This procedure is illustrated with the Mayo primary biliary
cholangitis dataset as given in package survival
.
## Mayo PBC data
library(survival)
pbc_analy <- as.matrix(na.omit(pbc[,c("time","status","age","edema","bili","albumin","protime")]))
# log transformation for time, bili, albumin, and protime
pbc_analy[,c(1,5:7)] <- log(pbc_analy[,c(1,5:7)])
colnames(pbc_analy)[c(1,5:7)] <- paste("log",colnames(pbc_analy)[c(1,5:7)])
# convert status to censoring indicator
pbc_analy[,2] <- pbc_analy[,2]>1
## Censored quantile regression
library(cequre)
taus <- .1*(1:9)
fit <- cequre(pbc_analy[,1],pbc_analy[,2],pbc_analy[,-c(1,2)],
taus=taus,res=200)
Plot the estimated regression coefficient processes and associated pointwise 95% CI’s:
var.names <- c("baseline",colnames(pbc_analy)[-(1:2)])
par(mfrow=c(3,2),mai=c(.5,.5,.1,.1))
for(k in 1:6) {
plot(fit$curve[7,],fit$curve[k,],type="s",lwd=2,xlab="",ylab="",
ylim=c(-1,1)*max(abs(fit$curve[k,])))
lines(c(0,1),c(0,0),lty=2)
text(1,max(abs(fit$curve[k,])),var.names[k],adj=c(1,1))
# pointwise 95% CI
for(i in 1:length(taus)) lines(c(1,1)*taus[i],fit$bt[k,i]+
c(-1,1)*qnorm(.975)*sqrt(fit$va[k,k,i]))
}
mtext("Estimated regression coefficient",side=2,line=-1,outer=T,cex=.8)
mtext("Probability",side=1,line=-1,outer=T,cex=.8)
The monotonicity-respecting restoration is illustrated with the preceding censored quantile regression using the Mayo primary biliary cholangitis dataset:
# zch needs to include constant 1 as being the covariate for the intercept
mfit <- monodr(fit$curve,
zch = cbind(rep(1,dim(pbc_analy)[1]),pbc_analy[,-(1:2)]),
initau=fit$tau.bnd/2)
# plot the original regression coefficient (in black) and the one after the
# monotonicity-respecting restoration (in orange)
par(mfrow=c(3,2),mai=c(.5,.5,.1,.1))
for(k in 1:6) {
plot(fit$curve[7,],fit$curve[k,],type="s",xlab="",ylab="",
ylim=c(-1,1)*max(abs(fit$curve[k,])))
lines(c(0,1),c(0,0),lty=2)
text(1,max(abs(fit$curve[k,])),var.names[k],adj=c(1,1))
lines(mfit$airc[7,],mfit$airc[k,],lwd=2,col="dark orange")
}
mtext("Estimated regression coefficient",side=2,line=-1,outer=T,cex=.8)
mtext("Probability",side=1,line=-1,outer=T,cex=.8)
Huang, Y. (2010) Quantile calculus and censored regression, The Annals of Statistics 38, 1607–1637.
Huang, Y. (2017) Restoration of monotonicity respecting in dynamic regression. Journal of the American Statistical Association 112, 613–622.