knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
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.
cequre is available on CRAN:
install.packages("cequre")
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.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.