1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
|
#SCCS @(#)plot.cox.zph.s 4.6 08/13/96
plot.cox.zph <- function(x, resid=TRUE, se=TRUE, df=4, nsmo=40, var, ...) {
##require(splines)
xx <- x$x
yy <- x$y
d <- nrow(yy)
df <- max(df) #error proofing
nvar <- ncol(yy)
pred.x <- seq(from=min(xx), to=max(xx), length=nsmo)
temp <- c(pred.x, xx)
lmat <- ns(temp, df=df, intercept=TRUE)
pmat <- lmat[1:nsmo,] # for prediction
xmat <- lmat[-(1:nsmo),]
qmat <- qr(xmat)
if (se) {
bk <- backsolve(qmat$qr[1:df, 1:df], diag(df))
xtx <- bk %*% t(bk)
seval <- d*((pmat%*% xtx) *pmat) %*% rep(1, df)
}
ylab <- paste("Beta(t) for", dimnames(yy)[[2]])
if (missing(var)) var <- 1:nvar
else {
if (is.character(var)) var <- match(var, dimnames(yy)[[2]])
if (any(is.na(var)) || max(var)>nvar || min(var) <1)
stop("Invalid variable requested")
}
#
# Figure out a 'good' set of x-axis labels. Find 8 equally spaced
# values on the 'transformed' axis. Then adjust until they correspond
# to rounded 'true time' values. Avoid the edges of the x axis, or
# approx() may give a missing value
if (x$transform == 'log') {
xx <- exp(xx)
pred.x <- exp(pred.x)
}
else if (x$transform != 'identity') {
xtime <- as.numeric(dimnames(yy)[[1]])
apr1 <- approx(xx, xtime, seq(min(xx), max(xx), length=17)[2*(1:8)])
temp <- signif(apr1$y,2)
apr2 <- approx(xtime, xx, temp)
xaxisval <- apr2$y
xaxislab <- rep("",8)
for (i in 1:8) xaxislab[i] <- format(temp[i])
}
for (i in var) {
y <- yy[,i]
yhat <- pmat %*% qr.coef(qmat, y)
if (resid) yr <-range(yhat, y)
else yr <-range(yhat)
if (se) {
temp <- 2* sqrt(x$var[i,i]*seval)
yup <- yhat + temp
ylow<- yhat - temp
yr <- range(yr, yup, ylow)
}
if (x$transform=='identity')
plot(range(xx), yr, type='n', xlab="Time", ylab=ylab[i], ...)
else if (x$transform=='log')
plot(range(xx), yr, type='n', xlab="Time", ylab=ylab[i], log='x',
...)
else {
plot(range(xx), yr, type='n', xlab="Time", ylab=ylab[i],
axes=FALSE,...)
axis(1, xaxisval, xaxislab)
axis(2)
box()
}
if (resid) points(xx, y)
lines(pred.x, yhat)
if (se) {
lines(pred.x, yup,lty=2)
lines(pred.x, ylow, lty=2)
}
}
}
|
