The R Project SVN
R-packages
(root)/trunk/Old/survival/Examples/testreg/testall.save – Rev 7500
Rev 104 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
S-PLUS : Copyright (c) 1988, 1996 MathSoft, Inc.S : Copyright AT&T.Version 3.4 Release 1 for Sun SPARC, SunOS 5.3 : 1996Mayo local startup...loading survivalloading dateloading slocal.miscset missing value action to 'na.omit'set contrasts to ('contr.treatment', 'contr.poly')stop automatic character -> factor conversion in data framesWorking 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-27voltage -53.9 5.545 -9.72 2.56e-22Log(scale) 4.8 0.105 45.56 0.00e+00Scale= 121Gaussian distributionLoglik(model)= -361.9 Loglik(intercept only)= -420.1Chisq= 116.33 on 1 degrees of freedom, p= 0Number of Newton-Raphson Iterations: 5n= 125Correlation of Coefficients:(Intercept) voltagevoltage -0.996Log(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-30voltage -62.21 5.967 -10.4 1.88e-25Log(scale) 4.53 0.108 41.9 0.00e+00Scale= 92.9Extreme value distributionLoglik(model)= -360 Loglik(intercept only)= -427.1Chisq= 134.25 on 1 degrees of freedom, p= 0Number of Newton-Raphson Iterations: 6n= 125Correlation of Coefficients:(Intercept) voltagevoltage -0.998Log(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-34voltage -55.48 4.986 -11.1 9.39e-29Log(scale) 4.19 0.117 35.8 2.03e-280Scale= 66.3Logistic distributionLoglik(model)= -360.4 Loglik(intercept only)= -423.7Chisq= 126.5 on 1 degrees of freedom, p= 0Number of Newton-Raphson Iterations: 5n= 125Correlation of Coefficients:(Intercept) voltagevoltage -0.996Log(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.00295Log(scale) -0.237 0.346 -0.687 0.49222Scale= 0.789Weibull distributionLoglik(model)= -51 Loglik(intercept only)= -51Number of Newton-Raphson Iterations: 9n= 210Correlation 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.684751Points out of bounds X= 2.814988 Y= -0.01247539Points 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.684751Lines out of bounds X= 2.814988 Y= -0.01247539Lines out of bounds X= 2.51267 Y= -0.182087Lines 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 foundDumped> > #> # 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)> > lfit0Call:survreg(formula = Surv(time, status) ~ 1, data = lung)Coefficients:(Intercept)6.034903Scale= 0.7593932Loglik(model)= -1153.9 Loglik(intercept only)= -1153.9n= 228> lfit1Call: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.00000age -0.00783 0.00687 0.00687 1.3 1 0.25000ridge(ph.ecog) -0.32032 0.08484 0.08405 14.2 1 0.00016Scale= 0.738Iterations: 1 outer, 4 Newton-RaphsonDegrees of freedom for terms= 1 1 1 1Likelihood ratio test=18.6 on 2 df, p=8.73e-05n=227 (1 observations deleted due to missing)> lfit2Call: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+00sex 0.40096 0.12371 0.12371 10.50 1 1.2e-03ridge(age) -0.00746 0.00675 0.00674 1.22 1 2.7e-01ridge(ph.ecog) -0.33848 0.08329 0.08314 16.51 1 4.8e-05Scale= 0.731Iterations: 1 outer, 5 Newton-RaphsonDegrees of freedom for terms= 1 1 2 1Likelihood ratio test=30 on 3 df, p=1.37e-06n=227 (1 observations deleted due to missing)> lfit3Call:survreg(formula = Surv(time, status) ~ sex + age + ph.ecog, data = lung)Coefficients:(Intercept) sex age ph.ecog6.273435 0.4010877 -0.007475331 -0.3396365Scale= 0.7311049Loglik(model)= -1132.4 Loglik(intercept only)= -1147.4Chisq= 29.98 on 3 degrees of freedom, p= 1.4e-06n=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)> > lfit4Call:survreg(formula = Surv(time, status) ~ xx, data = lung)Coefficients:(Intercept) xx1 xx2 xx3 xx4 xx513.5507 -7.615118 -7.423983 -7.532781 -7.570687 -14.52685Scale= 0.7557376Loglik(model)= -1150.1 Loglik(intercept only)= -1153.9Chisq= 7.52 on 5 degrees of freedom, p= 0.19n= 228> lfit5Call:survreg(formula = Surv(time, status) ~ age, data = lung)Coefficients:(Intercept) age6.887117 -0.01360819Scale= 0.7587492Loglik(model)= -1151.9 Loglik(intercept only)= -1153.9Chisq= 3.91 on 1 degrees of freedom, p= 0.048n= 228> lfit6Call: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.000pspline(age, df = 2), lin -0.0136 0.00687 0.00687 3.94 1.00 0.047pspline(age, df = 2), non 0.78 1.06 0.400Scale= 0.756Iterations: 4 outer, 9 Newton-RaphsonTheta= 0.926Degrees of freedom for terms= 0.4 2.1 1.0Likelihood 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-88factor(grp)2 0.786 0.177 4.44 8.79e-06factor(grp)3 0.728 0.183 3.97 7.09e-05factor(grp)4 -1.598 0.218 -7.32 2.48e-13factor(grp)5 -0.500 0.218 -2.29 2.21e-02factor(grp)6 0.475 0.170 2.79 5.23e-03Log(scale) -1.684 0.257 -6.54 6.09e-11Scale= 0.186Weibull distributionLoglik(model)= -26.7 Loglik(intercept only)= -40.7Chisq= 28.18 on 5 degrees of freedom, p= 3.4e-05Number of Newton-Raphson Iterations: 8n= 19Correlation of Coefficients:(Intercept) factor(grp)2 factor(grp)3 factor(grp)4 factor(grp)5factor(grp)2 -0.668factor(grp)3 -0.683 0.463factor(grp)4 -0.527 0.352 0.360factor(grp)5 -0.527 0.352 0.360 0.278factor(grp)6 -0.617 0.406 0.400 0.325 0.325Log(scale) -0.364 0.285 0.380 0.192 0.192factor(grp)6factor(grp)2factor(grp)3factor(grp)4factor(grp)5factor(grp)6Log(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)> ffitCall: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.0000frailty(grp, theta = 0.1) 10.4 2.15 0.0067Scale= 0.301Iterations: 1 outer, 6 Newton-RaphsonVariance of random effect= 0.1 EM likelihood = -11.8Degrees of freedom for terms= 0.3 2.2 0.7Likelihood 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.00000frailty(grp, df = 3, dist 16.4 3 0.00096Scale= 0.251Iterations: 7 outer, 24 Newton-RaphsonVariance of random effect= 0.197Degrees of freedom for terms= 0.1 3.0 0.6Likelihood 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 6N 19.0 3.000 6.00 6.00 1.00 1.00 2.00factor fit 2.3 0.018 0.80 0.75 -1.60 -0.48 0.49frailty 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.17654982 -1.56198696 -0.17654983 -0.87785012 -0.13012454 0.23871941 0.13369575 0.49650010 -0.17654986 -0.02255551 -0.13012457 -0.06575616 0.18012108 -0.31442628 0.18012109 -0.68264179 -0.176549810 0.08414645 0.180121011 -0.05034745 -0.130124512 0.51611042 0.133695713 -0.29527617 0.180121014 -0.07556559 0.133695715 1.43981531 -0.176549816 1.11151929 -0.130124517 1.47203044 -0.130124518 0.98566722 -0.176549819 0.19249326 0.133695720 1.01928857 0.133695721 -0.29177342 0.180121022 -1.56291591 -0.130124523 1.11035170 -0.176549824 -0.60115796 0.180121025 -0.70389695 0.133695726 -0.18273628 0.1336957attr($fit, "constant"):[1] 6.890663$se.fit:age ridge(ecog.ps, rx)1 0.356016316 0.17386872 0.403622262 0.17386873 0.226839188 0.19429114 0.061685835 0.18038855 0.128297162 0.17386876 0.005828413 0.19429117 0.016991596 0.18722828 0.081248722 0.18722829 0.176396749 0.173868710 0.021743704 0.187228211 0.013009936 0.194291112 0.133364529 0.180388513 0.076300276 0.187228214 0.019526381 0.180388515 0.372052731 0.173868716 0.287220023 0.194291117 0.380377220 0.194291118 0.254699460 0.173868719 0.049740854 0.180388520 0.263387319 0.180388521 0.075395153 0.187228222 0.403862307 0.194291123 0.286918315 0.173868724 0.155341078 0.187228225 0.181889150 0.180388526 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")> fit1Call:survreg(formula = Surv(time, status) ~ cage + cage^2, data = stanford2, dist ="lognormal")Coefficients:(Intercept) cage I(cage^2)6.717596 -0.06190903 -0.003504326Scale= 2.362866Loglik(model)= -863.6 Loglik(intercept only)= -868.8Chisq= 10.5 on 2 degrees of freedom, p= 0.0053n= 184> ldcase <- resid(fit1, type="ldcase")> ldresp <- resid(fit1, type="ldresp")> print(ldresp)139 159 181 119 74 120 990.1379203 0.145245 0.02628074 0.07320179 0.07624326 0.0399479 0.06328466108 179 43 134 160 177 1530.0612898 0.009685606 0.04767553 0.02980549 0.1036051 0.008990546 0.02114946136 133 176 66 157 114 460.0255769 0.1591464 0.008618358 0.03389346 0.01141316 0.01990885 0.0204497865 184 88 182 180 1630.02480539 1.085676e-05 0.0547439 0.001786473 0.002574794 0.00765404784 90 68 48 174 151 1250.02024456 0.08561197 0.03894985 0.07007566 0.0037674 0.008314653 0.0124855273 105 117 96 39 38 1060.01954895 0.01831982 0.01739301 0.01789441 0.02406183 0.02364314 0.0471718514 123 135 111 83 143 690.02051897 0.04763901 0.01663805 0.01367015 0.03204509 0.01857902 0.0205886827 113 167 156 141 300.03896725 0.03775025 0.005091493 0.01528402 0.008682116 0.01746136144 158 79 102 77 360.02593291 0.006620361 0.01375918 0.01547852 0.01786267 0.0233067183 122 162 121 87 23.720795e-05 0.01696469 0.005954799 0.01233286 0.01655939 0.108948964 150 85 71 19 21 1750.06015393 0.007469416 0.016665 0.01893414 0.02645489 0.18433 0.01789942169 148 138 98 104 103 120.004379942 0.007619682 0.009332594 0.014288 0.01445961 0.01449499 0.0340429989 3 100 55 142 63 1680.03358406 0.03113308 0.01412657 0.01179741 0.00864158 0.01426955 0.0045540372 137 10 124 17 94 820.01094162 0.009645953 0.01226565 0.01222511 0.01088511 0.01493685 0.0184422170 149 42 128 67 109 750.03988065 0.03038322 0.02127744 0.01439502 0.01285836 0.00894498 0.019977926 97 58 178 140 32 1260.02757124 0.02549339 0.02356049 0.002057497 0.01269584 0.01103393 0.012530351 101 29 33 164 600.0143023 0.01637414 0.02201027 0.01118993 0.006417554 0.008492277152 145 112 76 47 1180.008651511 0.009608663 0.01609214 0.02168279 0.02622512 0.022742765 129 31 35 40 1300.01184996 0.009391146 0.008772106 0.008526052 0.009451662 0.0129599728 56 91 44 23 37 700.01285987 0.01536642 0.02031498 0.02807957 0.01965943 0.01733256 0.009129009132 9 81 59 127 1310.009121716 0.009083024 0.01025238 0.01032187 0.01183694 0.0140329880 20 25 165 24 172 1460.02363944 0.02181251 0.02723395 0.02043511 0.02019542 0.01152649 0.0126590686 107 95 116 41 61 1550.01538527 0.02107502 0.0229847 0.02128395 0.01791007 0.01763098 0.01345059166 154 4 92 93 62 340.01285115 0.0121809 0.01470506 0.02599207 0.03098464 0.03037746 0.0216652215 173 171 52 110 500.01478523 0.007517959 0.008681577 0.01679632 0.02540017 0.0347067145 53 54 147 115 16 10.03229507 0.03017737 0.02416304 0.01870027 0.02172489 0.1164272 0.042577996 7 57 78 161 11 80.02459122 0.03585529 0.03587691 0.02865161 0.02603297 0.05640971 0.0433825149 13 22 180.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))> > fit1Call:survreg(formula = Surv(time, status) ~ strata(sex), data = lung)Coefficients:(Intercept)6.06217Scale:sex=1 sex=20.8167547 0.6533025Loglik(model)= -1152.5 Loglik(intercept only)= -1152.5n= 228> fit2Call:survreg(formula = Surv(time, status) ~ strata(sex) + sex, data = lung)Coefficients:(Intercept) sex5.49441 0.3801714Scale:sex=1 sex=20.8084286 0.6355802Loglik(model)= -1147.1 Loglik(intercept only)= -1152.5Chisq= 10.9 on 1 degrees of freedom, p= 0.00096n= 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)> fit1Call: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-16pspline(age, theta = 0.92 -0.0124 0.0067 0.0067 3.45 1.00 6.3e-02pspline(age, theta = 0.92 2.53 2.65 4.0e-01Scale:sex=1 sex=20.807 0.654Iterations: 1 outer, 3 Newton-RaphsonTheta= 0.92Degrees of freedom for terms= 0.5 3.6 2.0Likelihood ratio test=6.54 on 3.1 df, p=0.0937 n= 228> fit2Call: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-14pspline(age, theta = 0.92 -0.0111 0.00666 0.00666 2.77 1.00 9.6e-02pspline(age, theta = 0.92 2.46 2.68 4.2e-01sex 0.3686 0.11711 0.11685 9.91 1.00 1.6e-03Scale:sex=1 sex=20.8 0.636Iterations: 1 outer, 4 Newton-RaphsonTheta= 0.92Degrees of freedom for terms= 0.5 3.7 1.0 2.0Likelihood 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) rx4.983121 -0.2385013Scale= 0.2637831Loglik(model)= -242.3 Loglik(intercept only)= -246.3Chisq= 8 on 1 degrees of freedom, p= 0.0047n= 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"))> tempIntercept rx scalerfit0 4.983121 -0.2385013 0.2637831rfit1 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)^2surv surv.sparse cox cox.sparse0.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".