###################################### #### parametric regression models #### ###################################### ######### ### R ### ######### ############################## ## example 12.2 on page 399 ## ############################## library(survival) library(KMsurv) data(larynx) ## fit a weibull model to the whole data ## f1 <- survreg(Surv(time, delta)~ 1, data=larynx, dist="weibull") summary(f1) # gives MLE of mu and log(sigma) as well as Wald test # ## LR test on sigma=1 ## f1b <- survreg(Surv(time, delta)~ 1, data=larynx, dist="exponential") anova(f1, f1b) ## MLE of lambda and alpha ## mu <- f1\$coef sigma <- f1\$scale lambda <- exp(-mu/sigma) alpha <- 1/sigma ## covariance matrix of mu-sigma ## msvar <- f1\$var * matrix(c(1, sigma, sigma, sigma^2), nrow=2, byrow=TRUE) dimnames(msvar) <- list(c("mu", "sigma"), c("mu", "sigma")) ## variance and covariance of lambda and alpha using 12.2.7-12.2.9## var.alpha <- msvar[2,2]/sigma^4 var.lambda <- exp(-2*mu/sigma)*(msvar[1,1]/sigma^2+(mu^2)*msvar[2,2]/sigma^4-2*mu*msvar[1,2]/(sigma^3)) lacov <- exp(-mu/sigma)*(msvar[1,2]/sigma^3-mu*msvar[2,2]/sigma^4) ## median time ## mdn <- (-log(0.5)/lamda)^sigma der1 <- ((log(2))^sigma)*(-sigma)*(lambda^(-sigma-1)) #derivative wrt lambda der2 <- -log(log(2)/lambda)*((log(2)/lambda)^sigma)*(sigma^2) #derivative wrt sigma mdn.var <- var.lambda*der1^2 + var.alpha*der2^2 + 2*lacov*der1*der2 ## or... ## mdn2 <- predict(f1, type="quantile", p=0.5, se.fit=TRUE) mdn2\$fit[1] mdn2\$se[1] ## incorporate covariates in the model ## larynx <- within(larynx, { z1=(stage==2)*1 z2=(stage==3)*1 z3=(stage==4)*1 z4=age }) f2 <- survreg(Surv(time, delta)~ z1 + z2 + z3 + z4, data=larynx, dist="weibull") ## LR test on stage ## f2b <- survreg(Surv(time, delta)~ z4, data=larynx, dist="weibull") anova(f2, f2b) ## reparameterization (table 12.2) ## mu <- f2\$coef[1] sigma <- f2\$scale msvar <- f2\$var[c(1,6), c(1,6)] * matrix(c(1, sigma, sigma, sigma^2), nrow=2, byrow=TRUE) lambda <- exp(-mu/sigma) var.lambda <- exp(-2*mu/sigma)*(msvar[1,1]/sigma^2+(mu^2)*msvar[2,2]/sigma^4-2*mu*msvar[1,2]/(sigma^3)) alpha <- 1/sigma var.alpha <- msvar[2,2]/sigma^4 beta <- -f2\$coef[2:5]/sigma bscov <- f2\$var[2:5,6]*sigma var.beta <- diag(f2\$var)[2:5]/sigma^2 + (f2\$coef[2:5]^2)*msvar[2,2]/sigma^4 - 2*f2\$coef[2:5]*bscov/sigma^3 print(rbind(lambda=c(lambda, sqrt(var.lambda)), alpha=c(alpha, sqrt(var.alpha)), cbind(beta, sqrt(var.beta)))) ## model diagnosis ## ## CoxSnell residual ## cs <- lambda*exp(with(larynx, cbind(z1, z2, z3, z4)) %*% beta)*(larynx\$time^alpha) s.cs <- survfit(Surv(cs, larynx\$delta) ~ 1, type="fleming-harrington") plot(s.cs, mark.time=FALSE, fun="cumhaz", xlab="Cox-Snell Residual", ylab="Estimated Cumulative Hazard Rates", xlim=c(0,3), conf.int=FALSE) abline(a=0, b=1) ## deviance residual ## dres <- -resid(f2, type="deviance") #make the sign of deviance the same as martingale residual# riskscore <- predict(f2, type="linear") - f2\$coef[1] plot(riskscore, dres, pch=20) ************* *** STATA *** ************* insheet using larynx.csv stset time, failure(delta) ** AFT format ** streg z1 z2 z3 z4, dist(weibull) time ** PH format ** streg z1 z2 z3 z4, dist(weibull) nohr ** median time ** predict mdn, time median ** cox snell residual ** predict csres, csnell stset csres, fail(delta) sts gen H=na line H csres, sort ** deviance residual ** predict dres, deviance scatter dres time