Nothing
options(na.action=na.exclude) # preserve missings
options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
library(survival)
#
# Reproduce example 1 in the SAS lifereg documentation
#
# this fit doesn't give the same log-lik that they claim
fit1 <- survreg(Surv(time, status) ~ I(1000/(273.2+temp)), imotor,
subset=(temp>150), dist='lognormal')
summary(fit1)
# This one, with the loglik on the transformed scale (the inappropriate
# scale, Ripley & Venables would argue) does agree.
# All coefs are of course identical.
fit2 <- survreg(Surv(log(time), status) ~ I(1000/(273.2+temp)), imotor,
subset=(temp>150), dist='gaussian')
# Give the quantile estimates, which is the lower half of "output 48.1.5"
# in the SAS 9.2 manual
pp1 <- predict(fit1, newdata=list(temp=c(130,150)), p=c(.1, .5, .9),
type='quantile', se=T)
pp2 <- predict(fit1, newdata=list(temp=c(130,150)), p=c(.1, .5, .9),
type='uquantile', se=T)
pp1
temp130 <- matrix(0, nrow=3, ncol=6)
temp130[,1] <- pp1$fit[1,]
temp130[,2] <- pp1$se.fit[1,]
temp130[,3] <- pp2$fit[1,]
temp130[,4] <- pp2$se.fit[1,]
temp130[,5] <- exp(pp2$fit[1,] - 1.64*pp2$se.fit[1,])
temp130[,6] <- exp(pp2$fit[1,] + 1.64*pp2$se.fit[1,])
dimnames(temp130) <- list(c("p=.1", "p=.2", "p=.3"),
c("Time", "se(time)", "log(time)", "se[log(time)]",
"lower 90", "upper 90"))
print(temp130)
# A set of examples, copied from the manual pages of SAS procedure
# "reliability", which is part of their QC product.
#
color <- c("black", "red", "green", "blue", "magenta", "red4",
"orange", "DarkGreen", "cyan2", "DarkViolet")
palette(color)
pdf(file='reliability.pdf')
#
# Insulating fluids example
#
# Adding a -1 to the fit just causes the each group to have it's own
# intercept, rather than a global intercept + constrasts. The strata
# statement allows each to have a separate scale
ffit <- survreg(Surv(time) ~ factor(voltage) + strata(voltage) -1, ifluid)
# Get predicted quantiles at each of the voltages
# By default predict() would give a line of results for each observation,
# I only want the unique set of x's, i.e., only 4 cases
uvolt <- sort(unique(ifluid$voltage)) #the unique levels
plist <- c(1, 2, 5, 1:9 *10, 95, 99)/100
pred <- predict(ffit, type='quantile', p=plist,
newdata=data.frame(voltage=uvolt))
tfun <- function(x) log(-log(1-x))
matplot(t(pred), tfun(plist), type='l', log='x', lty=1,
col=1:4, yaxt='n',
xlab="Predicted", ylab="Quantile")
axis(2, tfun(plist), format(100*plist), adj=1)
kfit <- survfit(Surv(time) ~ voltage, ifluid, type='fleming') #KM fit
for (i in 1:4) {
temp <- kfit[i]
points(temp$time, tfun(1-temp$surv), col=i, pch=i)
}
# Now a table
temp <- array(0, dim=c(4,4,4)) #4 groups by 4 parameters by 4 stats
temp[,1,1] <- ffit$coef # "EV Location" in SAS manual
temp[,2,1] <- ffit$scale # "EV scale"
temp[,3,1] <- exp(ffit$coef) # "Weibull Scale"
temp[,4,1] <- 1/ffit$scale # "Weibull Shape"
temp[,1,2] <- sqrt(diag(ffit$var))[1:4] #standard error
temp[,2,2] <- sqrt(diag(ffit$var))[5:8] * ffit$scale
temp[,3,2] <- temp[,1,2] * temp[,3,1]
temp[,4,2] <- temp[,2,2] / (temp[,2,1])^2
temp[,1,3] <- temp[,1,1] - 1.96*temp[,1,2] #lower conf limits
temp[,1,4] <- temp[,1,1] + 1.96*temp[,1,2] # upper
# log(scale) is the natural parameter, in which the routine did its fitting
# and on which the std errors were computed
temp[,2, 3] <- exp(log(ffit$scale) - 1.96*sqrt(diag(ffit$var))[5:8])
temp[,2, 4] <- exp(log(ffit$scale) + 1.96*sqrt(diag(ffit$var))[5:8])
temp[,3, 3:4] <- exp(temp[,1,3:4])
temp[,4, 3:4] <- 1/temp[,2,4:3]
dimnames(temp) <- list(uvolt, c("EV Location", "EV Scale", "Weibull scale",
"Weibull shape"),
c("Estimate", "SE", "lower 95% CI", "uppper 95% CI"))
print(aperm(temp, c(2,3,1)), digits=5)
rm(temp, uvolt, plist, pred, ffit, kfit)
#####################################################################
# Turbine cracks data
crack2 <- with(cracks, data.frame(day1=c(NA, days), day2=c(days, NA),
n=c(fail, 167-sum(fail))))
cfit <- survreg(Surv(day1, day2, type='interval2') ~1,
dist='weibull', data=crack2, weights=n)
summary(cfit)
#Their output also has Wiebull scale = exp(cfit$coef), shape = 1/(cfit$scale)
# Draw the SAS plot
# The "type=fleming" argument reflects that they estimate hazards rather than
# survival, and forces a Nelson-Aalen hazard estimate
#
plist <- c(1, 2, 5, 1:8 *10)/100
plot(qsurvreg(plist, cfit$coef, cfit$scale), tfun(plist), log='x',
yaxt='n', type='l',
xlab="Weibull Plot for Time", ylab="Percent")
axis(2, tfun(plist), format(100*plist), adj=1)
kfit <- survfit(Surv(day1, day2, type='interval2') ~1, data=crack2,
weight=n, type='fleming')
# Only plot point where n.event > 0
# Why? I'm trying to match them. Personally, all should be plotted.
who <- (kfit$n.event > 0)
points(kfit$time[who], tfun(1-kfit$surv[who]), pch='+')
points(kfit$time[who], tfun(1-kfit$upper[who]), pch='-')
points(kfit$time[who], tfun(1-kfit$lower[who]), pch='-')
text(rep(3,6), seq(.5, -1.0, length=6),
c("Scale", "Shape", "Right Censored", "Left Censored",
"Interval Censored", "Fit"), adj=0)
text(rep(9,6), seq(.5, -1.0, length=6),
c(format(round(exp(cfit$coef), 2)),
format(round(1/cfit$scale, 2)),
format(tapply(crack2$n, cfit$y[,3], sum)), "ML"), adj=1)
# Now a portion of his percentiles table
# I don't get the same SE as SAS, I haven't checked out why. The
# estimates and se for the underlying Weibull model are the same.
temp <- predict(cfit, type='quantile', p=plist, se=T)
tempse <- sqrt(temp$se[1,])
mat <- cbind(temp$fit[1,], tempse,
temp$fit[1,] -1.96*tempse, temp$fit[1,] + 1.96*tempse)
dimnames(mat) <- list(plist*100, c("Estimate", "SE", "Lower .95", "Upper .95"))
print(mat)
#
# The cracks data has a particularly easy estimate, so use
# it to double check code
time <- c(crack2$day2[1], (crack2$day1 + crack2$day2)[2:8]/2,
crack2$day1[9])
cdf <- cumsum(crack2$n)/sum(crack2$n)
all.equal(kfit$time, time)
all.equal(kfit$surv, 1-cdf[c(1:8,8)])
rm(time, cdf, kfit)
#######################################################
#
# Replacement of valve seats in diesel engines
# The input data has id, time, and an indicator of whether there was an
# event at that time: 0=no, 1=yes. No one has an event at their last time.
# The input data has two engines with dual failures: 328 loses 2 valves at
# time 653, and number 402 loses 2 at time 139. For each, fudge the first
# time to be .1 days earlier.
#
ties <- which(diff(valveSeat$time)==0 & diff(valveSeat$id)==0)
temp <- valveSeat$time
temp[ties] <- temp[ties] - .1
n <- length(temp)
first <- !duplicated(valveSeat$id)
vtemp <- with(valveSeat, data.frame(id =id,
time1= ifelse(first, 0, c(0, temp[-n])),
time2= temp, status=status))
kfit <- survfit(Surv(time1, time2, status) ~1, vtemp, id=id)
plot(kfit, fun='cumhaz', ylab="Sample Mean Cumulative Failures", xlab='Time')
title("Valve replacement data")
# The summary.survfit function doesn't have an option for printing out
# cumulative hazards instead of survival --- need to add that
# so I just reprise the central code of print.summary.survfit
xx <- summary(kfit)
temp <- cbind(xx$time, xx$n.risk, xx$n.event, xx$cumhaz,
xx$std.chaz, -log(xx$upper), -log(xx$lower))
dimnames(temp) <- list(rep("", nrow(temp)),
c("time", "n.risk", "n.event", "Cum haz", "std.err",
"lower 95%", "upper 95%"))
print(temp, digits=2)
# Note that I have the same estimates but different SE's than SAS. We are using
# a different estimator. It's a statistical argument as to which is
# better (one could defend both sides): SAS the more standard estimate found
# in the reliability literature, mine the estimate from statistics literature
rm(temp, kfit, xx)
######################################################
# Turbine data, lognormal fit
turbine <- read.table('data.turbine',
col.names=c("time1", "time2", "n"))
tfit <- survreg(Surv(time1, time2, type='interval2') ~1, turbine,
dist='lognormal', weights=n, subset=(n>0))
summary(tfit)
# Now, do his plot, but put bootstrap confidence bands on it!
# First, make a simple data set without weights
tdata <- turbine[rep(1:nrow(turbine), turbine$n),]
qstat <- function(data) {
temp <- survreg(Surv(time1, time2, type='interval2') ~1, data=data,
dist='lognormal')
qsurvreg(plist, temp$coef, temp$scale, dist='lognormal')
}
{if (exists('bootstrap')) {
set.seed(1953) # a good year :-)
bfit <- bootstrap(tdata, qstat, B=1000)
bci <- limits.bca(bfit, probs=c(.025, .975))
}
else {
values <- matrix(0, nrow=1000, ncol=length(plist))
n <- nrow(tdata)
for (i in 1:1000) {
subset <- sample(1:n, n, replace=T)
values[i,] <- qstat(tdata[subset,])
}
bci <- t(apply(values,2, quantile, c(.05, .95)))
}
}
xmat <- cbind(qsurvreg(plist, tfit$coef, tfit$scale, dist='lognormal'),
bci)
matplot(xmat, qnorm(plist),
type='l', lty=c(1,2,2), col=c(1,1,1),
log='x', yaxt='n', ylab='Percent',
xlab='Time of Cracking (Hours x 100)')
axis(2, qnorm(plist), format(100*plist), adj=1)
title("Turbine Data")
kfit <- survfit(Surv(time1, time2, type='interval2') ~1, data=tdata)
points(kfit$time, qnorm(1-kfit$surv), pch='+')
dev.off() #close the plot file
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.