S-PLUS : Copyright (c) 1988, 1996 MathSoft, Inc. S : Copyright AT&T. Version 3.4 Release 1 for Sun SPARC, SunOS 5.3 : 1996 Mayo local startup... loading survival loading date loading slocal.misc set missing value action to 'na.omit' set contrasts to ('contr.treatment', 'contr.poly') stop automatic character -> factor conversion in data frames Working data will be in .Data > attach("../.Data") > dyn.load("../loadmod.o") > postscript(file="testall.ps") > options(na.action="na.omit", contrasts="contr.treatment") > # > # This data set caused problems for Splus 3.4 due to a mistake in > # my initial value code. Data courtesy Bob Treder at Statsci > # > capacitor <- read.table("data.capacitor", row.names=1, + col.names=c("", "days", "event", "voltage")) > > fitig <- survreg(Surv(days, event)~voltage, + dist = "gaussian", data = capacitor) > summary(fitig) Call: survreg(formula = Surv(days, event) ~ voltage, data = capacitor, dist = "gaussian") Value Std. Error z p (Intercept) 1764.9 163.387 10.80 3.36e-27 voltage -53.9 5.545 -9.72 2.56e-22 Log(scale) 4.8 0.105 45.56 0.00e+00 Scale= 121 Gaussian distribution Loglik(model)= -361.9 Loglik(intercept only)= -420.1 Chisq= 116.33 on 1 degrees of freedom, p= 0 Number of Newton-Raphson Iterations: 5 n= 125 Correlation of Coefficients: (Intercept) voltage voltage -0.996 Log(scale) 0.412 -0.384 > > fitix <- survreg(Surv(days, event)~voltage, + dist = "extreme", data = capacitor) > summary(fitix) Call: survreg(formula = Surv(days, event) ~ voltage, data = capacitor, dist = "extreme") Value Std. Error z p (Intercept) 2055.59 180.348 11.4 4.28e-30 voltage -62.21 5.967 -10.4 1.88e-25 Log(scale) 4.53 0.108 41.9 0.00e+00 Scale= 92.9 Extreme value distribution Loglik(model)= -360 Loglik(intercept only)= -427.1 Chisq= 134.25 on 1 degrees of freedom, p= 0 Number of Newton-Raphson Iterations: 6 n= 125 Correlation of Coefficients: (Intercept) voltage voltage -0.998 Log(scale) 0.425 -0.420 > > fitil <- survreg(Surv(days, event)~voltage, + dist = "logistic", data = capacitor) > summary(fitil) Call: survreg(formula = Surv(days, event) ~ voltage, data = capacitor, dist = "logistic") Value Std. Error z p (Intercept) 1811.56 148.853 12.2 4.48e-34 voltage -55.48 4.986 -11.1 9.39e-29 Log(scale) 4.19 0.117 35.8 2.03e-280 Scale= 66.3 Logistic distribution Loglik(model)= -360.4 Loglik(intercept only)= -423.7 Chisq= 126.5 on 1 degrees of freedom, p= 0 Number of Newton-Raphson Iterations: 5 n= 125 Correlation of Coefficients: (Intercept) voltage voltage -0.996 Log(scale) 0.343 -0.321 > > rm(fitil, fitig, fitix) > # > # Good initial values are key to this data set > # It killed v4 of survreg; > # data courtesy of Deborah Donnell, Fred Hutchinson Cancer Center > # > > donnell <- scan("data.donnell", what=list(time1=0, time2=0, status=0)) > donnell <- data.frame(donnell) > > dfit <- survreg(Surv(time1, time2, status, type="interval") ~1, donnell) > summary(dfit) Call: survreg(formula = Surv(time1, time2, status, type = "interval") ~ 1, data = donnell) Value Std. Error z p (Intercept) 2.390 0.804 2.973 0.00295 Log(scale) -0.237 0.346 -0.687 0.49222 Scale= 0.789 Weibull distribution Loglik(model)= -51 Loglik(intercept only)= -51 Number of Newton-Raphson Iterations: 9 n= 210 Correlation of Coefficients: (Intercept) Log(scale) 0.955 > > # > # Do a contour plot of the donnell data > # > npt <- 25 > beta0 <- seq(.4, 2.4, length=npt) > logsig <- seq(-1.4, 0.41, length=npt) > donlog <- matrix(0,npt, npt) > > for (i in 1:npt) { + for (j in 1:npt) { + fit <- survreg(Surv(time1, time2, status, type="interval") ~1, + donnell, init=c(beta0[i],logsig[j]), + control=list(maxiter=0));+ donlog[i,j] <- fit$log[1];+ };+ } > > clev <- -c(51, 51.5, 52:60, 65, 75, 85, 100, 150) > contour(beta0, logsig, pmax(donlog, -200), levels=clev, xlab="Intercept", + ylab="Log(sigma)") > points(2.39, log(.7885), pch=1, col=2) > title("Donnell data") > > # > # Compute the path of the iteration > # Step 2 isn't so good, and is followed by 3 iters of step-halving > # > niter <- 14 > donpath <- matrix(0,niter+1,2) > for (i in 0:niter){ + fit <- survreg(Surv(time1, time2, status, type="interval") ~1, + donnell, maxiter=i);+ donpath[i+1,] <- c(fit$coef, log(fit$scale));+ } Warning messages: 1: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights, offset, init = init, controlvals = controlvals, dist .... 2: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights, offset, init = init, controlvals = controlvals, dist .... 3: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights, offset, init = init, controlvals = controlvals, dist .... 4: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights, offset, init = init, controlvals = controlvals, dist .... 5: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights, offset, init = init, controlvals = controlvals, dist .... 6: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights, offset, init = init, controlvals = controlvals, dist .... 7: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights, offset, init = init, controlvals = controlvals, dist .... > points(donpath[,1], donpath[,2]) Points out of bounds X= 3.890545 Y= 0.684751 Points out of bounds X= 2.814988 Y= -0.01247539 Points out of bounds X= 2.51267 Y= -0.182087 > lines(donpath[,1], donpath[,2], col=4) Lines out of bounds X= 3.890545 Y= 0.684751 Lines out of bounds X= 2.814988 Y= -0.01247539 Lines out of bounds X= 2.51267 Y= -0.182087 Lines out of bounds X= 2.385003 Y= -0.2381598 > > rm(beta0, logsig, niter, fit, npt, donlog, clev) > #lfit1 <- censorReg(censor(time, status) ~ age + ph.ecog + strata(sex),lung) > lfit2 <- survreg(Surv(time, status) ~ age + ph.ecog + strata(sex), lung) > lfit3 <- survreg(Surv(time, status) ~ sex + (age+ph.ecog)*strata(sex), lung) > > lfit4 <- survreg(Surv(time, status) ~ age + ph.ecog , lung, + subset=(sex==1)) > lfit5 <- survreg(Surv(time, status) ~ age + ph.ecog , lung, + subset=(sex==2)) > > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y)) > #aeq(lfit4$coef, lfit1[[1]]$coef) > #aeq(lfit4$scale, lfit1[[1]]$scale) > aeq(c(lfit4$scale, lfit5$scale), lfit3$scale ) [1] "Mean relative difference: 1.364018e-07" > aeq(c(lfit4$scale, lfit5$scale), sapply(lfit1, function(x) x$scale)) Error: Object "lfit1" not found Dumped > > # > # Test out ridge regression and splines > # > lfit0 <- survreg(Surv(time, status) ~1, lung) > lfit1 <- survreg(Surv(time, status) ~ age + ridge(ph.ecog, theta=5), lung) > lfit2 <- survreg(Surv(time, status) ~ sex + ridge(age, ph.ecog, theta=1), lung) > lfit3 <- survreg(Surv(time, status) ~ sex + age + ph.ecog, lung) > > lfit0 Call: survreg(formula = Surv(time, status) ~ 1, data = lung) Coefficients: (Intercept) 6.034903 Scale= 0.7593932 Loglik(model)= -1153.9 Loglik(intercept only)= -1153.9 n= 228 > lfit1 Call: survreg(formula = Surv(time, status) ~ age + ridge(ph.ecog, theta = 5), data = lung) coef se(coef) se2 Chisq DF p (Intercept) 6.83082 0.42860 0.42860 254.0 1 0.00000 age -0.00783 0.00687 0.00687 1.3 1 0.25000 ridge(ph.ecog) -0.32032 0.08484 0.08405 14.2 1 0.00016 Scale= 0.738 Iterations: 1 outer, 4 Newton-Raphson Degrees of freedom for terms= 1 1 1 1 Likelihood ratio test=18.6 on 2 df, p=8.73e-05 n=227 (1 observations deleted due to missing) > lfit2 Call: survreg(formula = Surv(time, status) ~ sex + ridge(age, ph.ecog, theta = 1), data = lung) coef se(coef) se2 Chisq DF p (Intercept) 6.27163 0.45280 0.45210 191.84 1 0.0e+00 sex 0.40096 0.12371 0.12371 10.50 1 1.2e-03 ridge(age) -0.00746 0.00675 0.00674 1.22 1 2.7e-01 ridge(ph.ecog) -0.33848 0.08329 0.08314 16.51 1 4.8e-05 Scale= 0.731 Iterations: 1 outer, 5 Newton-Raphson Degrees of freedom for terms= 1 1 2 1 Likelihood ratio test=30 on 3 df, p=1.37e-06 n=227 (1 observations deleted due to missing) > lfit3 Call: survreg(formula = Surv(time, status) ~ sex + age + ph.ecog, data = lung) Coefficients: (Intercept) sex age ph.ecog 6.273435 0.4010877 -0.007475331 -0.3396365 Scale= 0.7311049 Loglik(model)= -1132.4 Loglik(intercept only)= -1147.4 Chisq= 29.98 on 3 degrees of freedom, p= 1.4e-06 n=227 (1 observations deleted due to missing) > > > xx <- pspline(lung$age, nterm=3, theta=.3) > xx <- matrix(unclass(xx), ncol=ncol(xx)) # the raw matrix > lfit4 <- survreg(Surv(time, status) ~xx, lung) > lfit5 <- survreg(Surv(time, status) ~age, lung) > > lfit6 <- survreg(Surv(time, status)~pspline(age, df=2), lung) > plot(lung$age, predict(lfit6), xlab="Age", ylab="Spline prediction") > title("Lung Data") > > lfit7 <- survreg(Surv(time, status) ~ offset(lfit6$lin), lung) > > lfit4 Call: survreg(formula = Surv(time, status) ~ xx, data = lung) Coefficients: (Intercept) xx1 xx2 xx3 xx4 xx5 13.5507 -7.615118 -7.423983 -7.532781 -7.570687 -14.52685 Scale= 0.7557376 Loglik(model)= -1150.1 Loglik(intercept only)= -1153.9 Chisq= 7.52 on 5 degrees of freedom, p= 0.19 n= 228 > lfit5 Call: survreg(formula = Surv(time, status) ~ age, data = lung) Coefficients: (Intercept) age 6.887117 -0.01360819 Scale= 0.7587492 Loglik(model)= -1151.9 Loglik(intercept only)= -1153.9 Chisq= 3.91 on 1 degrees of freedom, p= 0.048 n= 228 > lfit6 Call: survreg(formula = Surv(time, status) ~ pspline(age, df = 2), data = lung) coef se(coef) se2 Chisq DF p (Intercept) 6.5918 0.63681 0.41853 107.15 1.00 0.000 pspline(age, df = 2), lin -0.0136 0.00687 0.00687 3.94 1.00 0.047 pspline(age, df = 2), non 0.78 1.06 0.400 Scale= 0.756 Iterations: 4 outer, 9 Newton-Raphson Theta= 0.926 Degrees of freedom for terms= 0.4 2.1 1.0 Likelihood ratio test=5.2 on 1.5 df, p=0.0441 n= 228 > lfit7$coef (Intercept) -6.008005e-07 > > rm(lfit1, lfit2, lfit3, lfit4, lfit5, lfit6, lfit7) > rm(xx, lfit0) > # > # Data courtesy of Bercedis Peterson, Duke University. > # v4 of survreg fails due to 2 groups that have only 1 subject; the coef > # for them easily gets out of hand. In fact, this data set is my toughest > # test of the minimizer. > # > # A shrinkage model for this coefficient is therefore interesting > > > peterson <- data.frame( + scan("data.peterson", what=list(grp=0, time=0, status=0))) > > fitp <- survreg(Surv(time, status) ~ factor(grp), peterson) > summary(fitp) Call: survreg(formula = Surv(time, status) ~ factor(grp), data = peterson) Value Std. Error z p (Intercept) 2.291 0.115 19.92 2.93e-88 factor(grp)2 0.786 0.177 4.44 8.79e-06 factor(grp)3 0.728 0.183 3.97 7.09e-05 factor(grp)4 -1.598 0.218 -7.32 2.48e-13 factor(grp)5 -0.500 0.218 -2.29 2.21e-02 factor(grp)6 0.475 0.170 2.79 5.23e-03 Log(scale) -1.684 0.257 -6.54 6.09e-11 Scale= 0.186 Weibull distribution Loglik(model)= -26.7 Loglik(intercept only)= -40.7 Chisq= 28.18 on 5 degrees of freedom, p= 3.4e-05 Number of Newton-Raphson Iterations: 8 n= 19 Correlation of Coefficients: (Intercept) factor(grp)2 factor(grp)3 factor(grp)4 factor(grp)5 factor(grp)2 -0.668 factor(grp)3 -0.683 0.463 factor(grp)4 -0.527 0.352 0.360 factor(grp)5 -0.527 0.352 0.360 0.278 factor(grp)6 -0.617 0.406 0.400 0.325 0.325 Log(scale) -0.364 0.285 0.380 0.192 0.192 factor(grp)6 factor(grp)2 factor(grp)3 factor(grp)4 factor(grp)5 factor(grp)6 Log(scale) 0.083 > > # Now a shrinkage model. Give the group coefficients > # about 1/2 the scale parameter of the original model, i.e., .18. > # > ffit <- survreg(Surv(time, status) ~ frailty(grp, theta=.1), peterson) > ffit Call: survreg(formula = Surv(time, status) ~ frailty(grp, theta = 0.1), data = peterson) coef se(coef) se2 Chisq DF p (Intercept) 2.62 0.172 0.0874 232.0 1.00 0.0000 frailty(grp, theta = 0.1) 10.4 2.15 0.0067 Scale= 0.301 Iterations: 1 outer, 6 Newton-Raphson Variance of random effect= 0.1 EM likelihood = -11.8 Degrees of freedom for terms= 0.3 2.2 0.7 Likelihood ratio test=13.8 on 1.1 df, p=0.00027 n= 19 > > # > # Try 3 degrees of freedom Gaussian fit, since there are 6 groups. > # Compare them to the unconstrained ones. The frailty coefs are > # on a "sum to 0" constraint rather than "first coef=0", so > # some conversion is neccessary > # > ffit3 <- survreg(Surv(time, status) ~ frailty(grp, df=3, dist="gauss"), + peterson) > print(ffit3) Call: survreg(formula = Surv(time, status) ~ frailty(grp, df = 3, dist = "gauss"), data = peterson) coef se(coef) se2 Chisq DF p (Intercept) 2.44 0.223 0.0661 119.8 1 0.00000 frailty(grp, df = 3, dist 16.4 3 0.00096 Scale= 0.251 Iterations: 7 outer, 24 Newton-Raphson Variance of random effect= 0.197 Degrees of freedom for terms= 0.1 3.0 0.6 Likelihood ratio test=20.1 on 1.7 df, p=2.79e-05 n= 19 > > temp <- mean(c(0, fitp$coef[-1])) > temp2 <- c(fitp$coef[1] + temp, c(0,fitp$coef[-1]) - temp) > xx <- rbind(c(nrow(peterson), table(peterson$grp)), + temp2, + c(ffit3$coef, ffit3$frail)) > dimnames(xx) <- list(c("N", "factor fit", "frailty fit"), + c("Intercept", paste("grp", 1:6))) > signif(xx,2) Intercept grp 1 grp 2 grp 3 grp 4 grp 5 grp 6 N 19.0 3.000 6.00 6.00 1.00 1.00 2.00 factor fit 2.3 0.018 0.80 0.75 -1.60 -0.48 0.49 frailty fit 2.4 -0.180 0.58 0.55 -0.77 -0.44 0.26 > # > # All but the first coef are shrunk towards zero. > # > rm(ffit, ffit3, temp, temp2, xx, fitp) > > # > # Look at predicted values > # > ofit1 <- survreg(Surv(futime, fustat) ~ age + ridge(ecog.ps, rx), ovarian) > > predict(ofit1) [1] 207.7546 172.7985 358.7725 1426.6414 1353.7225 843.8571 1102.1610 [8] 859.5061 416.3272 1280.4037 820.7276 1882.7133 876.1244 1041.8917 [15] 3477.0123 2622.9581 3761.4852 2207.8635 1362.1943 3113.9504 879.1986 [22] 180.8418 2501.0478 645.2425 555.8297 936.0066 > predict(ofit1, type="response") [1] 207.7546 172.7985 358.7725 1426.6414 1353.7225 843.8571 1102.1610 [8] 859.5061 416.3272 1280.4037 820.7276 1882.7133 876.1244 1041.8917 [15] 3477.0123 2622.9581 3761.4852 2207.8635 1362.1943 3113.9504 879.1986 [22] 180.8418 2501.0478 645.2425 555.8297 936.0066 > predict(ofit1, type="terms", se=T) $fit: age ridge(ecog.ps, rx) 1 -1.37775562 -0.1765498 2 -1.56198696 -0.1765498 3 -0.87785012 -0.1301245 4 0.23871941 0.1336957 5 0.49650010 -0.1765498 6 -0.02255551 -0.1301245 7 -0.06575616 0.1801210 8 -0.31442628 0.1801210 9 -0.68264179 -0.1765498 10 0.08414645 0.1801210 11 -0.05034745 -0.1301245 12 0.51611042 0.1336957 13 -0.29527617 0.1801210 14 -0.07556559 0.1336957 15 1.43981531 -0.1765498 16 1.11151929 -0.1301245 17 1.47203044 -0.1301245 18 0.98566722 -0.1765498 19 0.19249326 0.1336957 20 1.01928857 0.1336957 21 -0.29177342 0.1801210 22 -1.56291591 -0.1301245 23 1.11035170 -0.1765498 24 -0.60115796 0.1801210 25 -0.70389695 0.1336957 26 -0.18273628 0.1336957 attr($fit, "constant"): [1] 6.890663 $se.fit: age ridge(ecog.ps, rx) 1 0.356016316 0.1738687 2 0.403622262 0.1738687 3 0.226839188 0.1942911 4 0.061685835 0.1803885 5 0.128297162 0.1738687 6 0.005828413 0.1942911 7 0.016991596 0.1872282 8 0.081248722 0.1872282 9 0.176396749 0.1738687 10 0.021743704 0.1872282 11 0.013009936 0.1942911 12 0.133364529 0.1803885 13 0.076300276 0.1872282 14 0.019526381 0.1803885 15 0.372052731 0.1738687 16 0.287220023 0.1942911 17 0.380377220 0.1942911 18 0.254699460 0.1738687 19 0.049740854 0.1803885 20 0.263387319 0.1803885 21 0.075395153 0.1872282 22 0.403862307 0.1942911 23 0.286918315 0.1738687 24 0.155341078 0.1872282 25 0.181889150 0.1803885 26 0.047219621 0.1803885 > > temp1 <- predict(ofit1, se=T) > temp2 <- predict(ofit1, type= "response", se=T) > all.equal(temp2$se.fit, temp1$se.fit* sqrt(exp(temp1$fit))) [1] "Mean relative difference: Inf" > # > # The Stanford data from 1980 is used in Escobar and Meeker > # t5 = T5 mismatch score > # Their case numbers correspond to a data set sorted by age > # > stanford2 <- read.table("data.stanford", + col.names=c("id", "time", "status", "age", "t5")) > > stanford2$t5 <- ifelse(stanford2$t5 <0, NA, stanford2$t5) > stanford2 <- stanford2[order(stanford2$age, stanford2$time),] > stanford2$time <- ifelse(stanford2$time==0, .5, stanford2$time) > > cage <- stanford2$age - mean(stanford2$age) > fit1 <- survreg(Surv(time, status) ~ cage + cage^2, stanford2, + dist="lognormal") > fit1 Call: survreg(formula = Surv(time, status) ~ cage + cage^2, data = stanford2, dist = "lognormal") Coefficients: (Intercept) cage I(cage^2) 6.717596 -0.06190903 -0.003504326 Scale= 2.362866 Loglik(model)= -863.6 Loglik(intercept only)= -868.8 Chisq= 10.5 on 2 degrees of freedom, p= 0.0053 n= 184 > ldcase <- resid(fit1, type="ldcase") > ldresp <- resid(fit1, type="ldresp") > print(ldresp) 139 159 181 119 74 120 99 0.1379203 0.145245 0.02628074 0.07320179 0.07624326 0.0399479 0.06328466 108 179 43 134 160 177 153 0.0612898 0.009685606 0.04767553 0.02980549 0.1036051 0.008990546 0.02114946 136 133 176 66 157 114 46 0.0255769 0.1591464 0.008618358 0.03389346 0.01141316 0.01990885 0.02044978 65 184 88 182 180 163 0.02480539 1.085676e-05 0.0547439 0.001786473 0.002574794 0.007654047 84 90 68 48 174 151 125 0.02024456 0.08561197 0.03894985 0.07007566 0.0037674 0.008314653 0.01248552 73 105 117 96 39 38 106 0.01954895 0.01831982 0.01739301 0.01789441 0.02406183 0.02364314 0.04717185 14 123 135 111 83 143 69 0.02051897 0.04763901 0.01663805 0.01367015 0.03204509 0.01857902 0.02058868 27 113 167 156 141 30 0.03896725 0.03775025 0.005091493 0.01528402 0.008682116 0.01746136 144 158 79 102 77 36 0.02593291 0.006620361 0.01375918 0.01547852 0.01786267 0.0233067 183 122 162 121 87 2 3.720795e-05 0.01696469 0.005954799 0.01233286 0.01655939 0.1089489 64 150 85 71 19 21 175 0.06015393 0.007469416 0.016665 0.01893414 0.02645489 0.18433 0.01789942 169 148 138 98 104 103 12 0.004379942 0.007619682 0.009332594 0.014288 0.01445961 0.01449499 0.03404299 89 3 100 55 142 63 168 0.03358406 0.03113308 0.01412657 0.01179741 0.00864158 0.01426955 0.00455403 72 137 10 124 17 94 82 0.01094162 0.009645953 0.01226565 0.01222511 0.01088511 0.01493685 0.0184422 170 149 42 128 67 109 75 0.03988065 0.03038322 0.02127744 0.01439502 0.01285836 0.00894498 0.0199779 26 97 58 178 140 32 126 0.02757124 0.02549339 0.02356049 0.002057497 0.01269584 0.01103393 0.0125303 51 101 29 33 164 60 0.0143023 0.01637414 0.02201027 0.01118993 0.006417554 0.008492277 152 145 112 76 47 118 0.008651511 0.009608663 0.01609214 0.02168279 0.02622512 0.02274276 5 129 31 35 40 130 0.01184996 0.009391146 0.008772106 0.008526052 0.009451662 0.01295997 28 56 91 44 23 37 70 0.01285987 0.01536642 0.02031498 0.02807957 0.01965943 0.01733256 0.009129009 132 9 81 59 127 131 0.009121716 0.009083024 0.01025238 0.01032187 0.01183694 0.01403298 80 20 25 165 24 172 146 0.02363944 0.02181251 0.02723395 0.02043511 0.02019542 0.01152649 0.01265906 86 107 95 116 41 61 155 0.01538527 0.02107502 0.0229847 0.02128395 0.01791007 0.01763098 0.01345059 166 154 4 92 93 62 34 0.01285115 0.0121809 0.01470506 0.02599207 0.03098464 0.03037746 0.02166522 15 173 171 52 110 50 0.01478523 0.007517959 0.008681577 0.01679632 0.02540017 0.03470671 45 53 54 147 115 16 1 0.03229507 0.03017737 0.02416304 0.01870027 0.02172489 0.1164272 0.04257799 6 7 57 78 161 11 8 0.02459122 0.03585529 0.03587691 0.02865161 0.02603297 0.05640971 0.04338251 49 13 22 18 0.03425475 0.06262793 0.1029315 0.1442429 > # The ldcase and ldresp should be compared to table 1 in Escobar and > # Meeker, Biometrics 1992, p519; the colum they label as (1/2) A_{ii} > > plot1 <- function() { + # make their figure 1, 2, and 6 + plot(stanford2$age, stanford2$time, log="y", xlab="Age", ylab="Days", + ylim=c(.01, 10^6));+ temp <- predict(fit1, type="response", se.fit=T) ;+ matlines(stanford2$age, cbind(temp$fit, temp$fit-1.96*temp$se.fit, + temp$fit+1.96*temp$se.fit), + lty=c(1,2,2));+ # these are the wrong CI lines, he plotted std dev, I plotted std err + # here are the right ones + # Using uncentered age gives different coefs, but makes prediction over an + # extended range somewhat simpler + refit <- survreg(Surv(time,status)~ age + age^2, stanford2, + dist="lognormal");+ plot(stanford2$age, stanford2$time, log="y", xlab="Age", ylab="Days", + ylim=c(.01, 10^6), xlim=c(0,75));+ temp2 <- predict(refit, list(age=1:75), type="quantile", p=c(.05, .5, .95));+ matlines(1:75, temp2, lty=c(1,2,2), col=2);+ + tsplot(ldcase, xlab="Case Number", ylab="(1/2) A");+ title (main="Case weight pertubations");+ tsplot(ldresp, xlab="Case Number", ylab="(1/2) A");+ title(main="Response pertubations");+ } > > plot1() Warning: Data values <=0 omitted from logarithmic plot > # > # Stanford predictions in other ways > # > fit2 <- survreg(Surv(time, status) ~ poly(age,2), stanford2, + dist="lognormal") > > p1 <- predict(fit1, type="response") > p2 <- predict(fit2, type="response") > aeq(p1, p2) [1] T > > p3 <- predict(fit2, type="terms", se=T) > p4 <- predict(fit2, type="lp", se=T) > p5 <- predict(fit1, type="lp", se=T) > aeq(p3$fit + attr(p3$fit, "constant"), p4$fit) [1] T > aeq(p4$fit, p5$fit) [1] T > aeq(p3$se.fit, p4$se.fit) #this one should be false [1] "Mean relative difference: 0.358807" > aeq(p4$se.fit, p5$se.fit) #this one true [1] T > > # > # Verify that scale can be fixed at a value > # coefs will differ slightly due to different iteration paths > tol <- survreg.control()$rel.tolerance > > # Intercept only models > fit1 <- survreg(Surv(time,status) ~ 1, lung) > fit2 <- survreg(Surv(time,status) ~ 1, lung, scale=fit1$scale) > all.equal(fit1$coef, fit2$coef, tolerance= tol) [1] T > all.equal(fit1$loglik, fit2$loglik, tolerance= tol) [1] T > > # multiple covariates > fit1 <- survreg(Surv(time,status) ~ age + ph.karno, lung) > fit2 <- survreg(Surv(time,status) ~ age + ph.karno, lung, + scale=fit1$scale) > all.equal(fit1$coef, fit2$coef, tolerance=tol) [1] T > all.equal(fit1$loglik[2], fit2$loglik[2], tolerance=tol) [1] T > > # penalized models > fit1 <- survreg(Surv(time, status) ~ pspline(age), lung) > fit2 <- survreg(Surv(time, status) ~ pspline(age), lung, scale=fit1$scale) > all.equal(fit1$coef, fit2$coef, tolerance=tol) [1] "Mean relative difference: 0.0002487206" > all.equal(fit1$loglik[2], fit2$loglik[2], tolerance=tol) [1] T > > rm(fit1, fit2, tol) > > # > # Test out the strata capabilities > # > tol <- survreg.control()$rel.tolerance > aeq <- function(x,y,...) all.equal(as.vector(x), as.vector(y), ...) > > # intercept only models > fit1 <- survreg(Surv(time, status) ~ strata(sex), lung) > fit2 <- survreg(Surv(time, status) ~ strata(sex) + sex, lung) > fit3a<- survreg(Surv(time,status) ~1, lung, subset=(sex==1)) > fit3b<- survreg(Surv(time,status) ~1, lung, subset=(sex==2)) > > fit1 Call: survreg(formula = Surv(time, status) ~ strata(sex), data = lung) Coefficients: (Intercept) 6.06217 Scale: sex=1 sex=2 0.8167547 0.6533025 Loglik(model)= -1152.5 Loglik(intercept only)= -1152.5 n= 228 > fit2 Call: survreg(formula = Surv(time, status) ~ strata(sex) + sex, data = lung) Coefficients: (Intercept) sex 5.49441 0.3801714 Scale: sex=1 sex=2 0.8084286 0.6355802 Loglik(model)= -1147.1 Loglik(intercept only)= -1152.5 Chisq= 10.9 on 1 degrees of freedom, p= 0.00096 n= 228 > aeq(fit2$scale, c(fit3a$scale, fit3b$scale), tolerance=tol) [1] T > aeq(fit2$loglik[2], (fit3a$loglik + fit3b$loglik)[2], tolerance=tol) [1] T > aeq(fit2$coef[1] + 1:2*fit2$coef[2], c(fit3a$coef, fit3b$coef), tolerance=tol) [1] T > > #penalized models > fit1 <- survreg(Surv(time, status) ~ pspline(age, theta=.92)+strata(sex), lung) > fit2 <- survreg(Surv(time, status) ~ pspline(age, theta=.92)+ + strata(sex) + sex, lung) > fit1 Call: survreg(formula = Surv(time, status) ~ pspline(age, theta = 0.92) + strata(sex), data = lung) coef se(coef) se2 Chisq DF p (Intercept) 6.9036 0.8469 0.5688 66.45 1.00 3.3e-16 pspline(age, theta = 0.92 -0.0124 0.0067 0.0067 3.45 1.00 6.3e-02 pspline(age, theta = 0.92 2.53 2.65 4.0e-01 Scale: sex=1 sex=2 0.807 0.654 Iterations: 1 outer, 3 Newton-Raphson Theta= 0.92 Degrees of freedom for terms= 0.5 3.6 2.0 Likelihood ratio test=6.54 on 3.1 df, p=0.0937 n= 228 > fit2 Call: survreg(formula = Surv(time, status) ~ pspline(age, theta = 0.92) + strata(sex) + sex, data = lung) coef se(coef) se2 Chisq DF p (Intercept) 6.3729 0.84471 0.59118 56.92 1.00 4.5e-14 pspline(age, theta = 0.92 -0.0111 0.00666 0.00666 2.77 1.00 9.6e-02 pspline(age, theta = 0.92 2.46 2.68 4.2e-01 sex 0.3686 0.11711 0.11685 9.91 1.00 1.6e-03 Scale: sex=1 sex=2 0.8 0.636 Iterations: 1 outer, 4 Newton-Raphson Theta= 0.92 Degrees of freedom for terms= 0.5 3.7 1.0 2.0 Likelihood ratio test=16.8 on 4.2 df, p=0.00245 n= 228 > > age1 <- ifelse(lung$sex==1, lung$age, mean(lung$age)) > age2 <- ifelse(lung$sex==2, lung$age, mean(lung$age)) > fit3 <- survreg(Surv(time,status) ~ pspline(age1, theta=.92) + + pspline(age2, theta=.95) + sex + strata(sex), lung, + rel.tol=1e-6) > fit3a<- survreg(Surv(time,status) ~pspline(age, theta=.92), lung, + subset=(sex==1)) > fit3b<- survreg(Surv(time,status) ~pspline(age, theta=.95), lung, + subset=(sex==2)) > > # relax the tolerance a little, since the above has lots of parameters > # I still don't exactly match the second group, but very close > aeq(fit3$scale, c(fit3a$scale, fit3b$scale), tolerance=tol*10) [1] "Mean relative difference: 0.001270825" > aeq(fit3$loglik[2], (fit3a$loglik + fit3b$loglik)[2], tolerance=tol*10) [1] T > pred <- predict(fit3) > aeq(pred[lung$sex==1] , predict(fit3a), tolerance=tol*10) [1] T > aeq(pred[lung$sex==2], predict(fit3b), tolerance=tol*10) [1] "Mean relative difference: 0.01158256" > > > > > # > # Some tests using the rat data > # > rats <- read.table("../testfrail/data.rats", + col.names=c("litter", "rx", "time", "status")) > > rfitnull <- survreg(Surv(time, status) ~1, rats) > temp <- rfitnull$scale^2 * pi^2/6 > cat("Effective n =", round(temp*(solve(rfitnull$var))[1,1],1), "\n") Effective n = 65.8 > > rfit0 <- survreg(Surv(time, status) ~ rx , rats) > print(rfit0) Call: survreg(formula = Surv(time, status) ~ rx, data = rats) Coefficients: (Intercept) rx 4.983121 -0.2385013 Scale= 0.2637831 Loglik(model)= -242.3 Loglik(intercept only)= -246.3 Chisq= 8 on 1 degrees of freedom, p= 0.0047 n= 150 > > rfit1 <- survreg(Surv(time, status) ~ rx + factor(litter), rats) > temp <- rbind(c(rfit0$coef, rfit0$scale), c(rfit1$coef[1:2], rfit1$scale)) > dimnames(temp) <- list(c("rfit0", "rfit1"), c("Intercept", "rx", "scale")) > temp Intercept rx scale rfit0 4.983121 -0.2385013 0.2637831 rfit1 4.902437 -0.2189405 0.2025429 > > > rfit2a <- survreg(Surv(time, status) ~ rx + + frailty.gaussian(litter, df=13, sparse=F), rats ) > rfit2b <- survreg(Surv(time, status) ~ rx + + frailty.gaussian(litter, df=13, sparse=T), rats ) > > rfit3a <- coxph(Surv(time,status) ~ rx + + frailty.gaussian(litter, df=13, sparse=F), rats ) > rfit3b <- coxph(Surv(time,status) ~ rx + + frailty(litter, df=13, dist="gauss"), rats) > > temp <- cbind(rfit2a$coef[3:52], rfit2b$frail, rfit3a$coef[2:51], rfit3b$frail) > dimnames(temp) <- list(NULL, c("surv","surv.sparse","cox","cox.sparse")) > pairs(temp) > apply(temp,2,var)/c(rfit2a$scale, rfit2b$scale, 1,1)^2 surv surv.sparse cox cox.sparse 0.1346009 0.1346009 0.1224049 0.1207863 > apply(temp,2,mean) surv surv.sparse cox cox.sparse -7.979728e-19 6.938894e-20 -1.096345e-17 1.054712e-17 > > # The parametric model gives the coefficients less variance for the > # two fits, for the same df, but the scaled results are similar. > # 13 df is near to the rmle for the rats > > rm(temp, rfit2a, rfit2b, rfit3a, rfit3b, rfitnull, rfit0, rfit1) > > temp <- matrix(scan("data.mpip", skip=23), ncol=13, byrow=T) > dimnames(temp) <- list(NULL, c("ved", "angina", "education", "prior.mi", + "nyha", "rales", "ef", "ecg", "angina2", "futime", + "status", "admit", "betab")) > > mpip <- data.frame(temp) > lved <- log(mpip$ved + .02) > > fit1 <- coxph(Surv(futime, status) ~ pspline(lved) + factor(nyha) + + rales + pspline(ef), mpip) > > temp <- predict(fit1, type="terms", se.fit=T) > yy <- cbind(temp$fit[,4], temp$fit[,4] + 1.96*temp$se[,4], + temp$fit[,4] - 1.96*temp$se[,4]) > index <- order(mpip$ef) > matplot(mpip$ef[index], yy[index,], type="l", lty=c(1,2,2), col=1) > title(xlab="Ejection Fraction", ylab="Cox model risk", + main="Post-Infarction Survival") > > fit2 <- coxph(Surv(futime, status) ~ lved + factor(nyha) + rales + + pspline(ef, df=0), mpip) Warning messages: 1: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <- (df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 .... 2: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <- (df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 .... 3: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <- (df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 .... 4: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <- (df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 .... 5: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <- (df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 .... 6: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <- (df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 .... 7: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <- (df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 .... 8: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <- (df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 .... 9: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <- (df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 .... 10: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <- (df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 .... > temp <- predict(fit2, type="terms", se.fit=T) > yy <- cbind(temp$fit[,4], temp$fit[,4] + 1.96*temp$se[,4], + temp$fit[,4] - 1.96*temp$se[,4]) > matplot(mpip$ef[index], yy[index,], type="l", lty=c(1,2,2), col=1) > title(xlab="Ejection Fraction", ylab="Cox model risk", + main="Post-Infarction Survival, AIC") > > > fit3 <- survreg(Surv(futime, status) ~ lved + factor(nyha) + rales + + pspline(ef, df=2), mpip, dist="lognormal") > temp <- predict(fit3, type="terms", se.fit=T) > yy <- cbind(temp$fit[,4], temp$fit[,4] + 1.96*temp$se[,4], + temp$fit[,4] - 1.96*temp$se[,4]) > matplot(mpip$ef[index], yy[index,], type="l", lty=c(1,2,2), col=1) > title(xlab="Ejection Fraction", ylab="Log-normal model predictor", + main="Post-Infarction Survival") > q() Generated postscript file "testall.ps".