The R Project SVN
R-packages
(root)/trunk/Old/survival/Examples/testreg/Rtest.out.save – Rev 7500
Rev 4821 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R version 2.7.0 Under development (unstable) (2007-10-09 r43132)Copyright (C) 2007 The R Foundation for Statistical ComputingISBN 3-900051-07-0R is free software and comes with ABSOLUTELY NO WARRANTY.You are welcome to redistribute it under certain conditions.Type 'license()' or 'licence()' for distribution details.Natural language support but running in an English localeR is a collaborative project with many contributors.Type 'contributors()' for more information and'citation()' on how to cite R or R packages in publications.Type 'demo()' for some demos, 'help()' for on-line help, or'help.start()' for an HTML browser interface to help.Type 'q()' to quit R.> #attach("../.Data")> #dyn.load('../loadmod.o')> postscript(file='Rtest.ps')> library(survival)> if (R.version$minor>2) options(na.action="na.exclude")> #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: 6n= 125>> 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.349 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: 7n= 125>> 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: 6n= 125>> 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.972 0.00295Log(scale) -0.237 0.346 -0.687 0.49232Scale= 0.789Weibull distributionLoglik(model)= -51 Loglik(intercept only)= -51Number of Newton-Raphson Iterations: 10n= 210>> #> # 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))+ }> points(donpath[,1], donpath[,2])> lines(donpath[,1], donpath[,2], col=4)>> rm(beta0, logsig, niter, fit, npt, donlog, clev)> #lfit1 <- censorReg(censor(time, status) ~ age + ph.ecog + strata(sex),lung)> data(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] TRUE> #aeq(c(lfit4$scale, lfit5$scale), sapply(lfit1, function(x) x$scale))>> #> # 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.034904Scale= 0.7593936Loglik(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, 5 Newton-RaphsonDegrees of freedom for terms= 1 1 1 1Likelihood ratio test=18.6 on 2 df, p=8.73e-05n=227 (1 observation deleted due to missingness)> 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, 6 Newton-RaphsonDegrees of freedom for terms= 1 1 2 1Likelihood ratio test=30 on 3 df, p=1.37e-06n=227 (1 observation deleted due to missingness)> lfit3Call:survreg(formula = Surv(time, status) ~ sex + age + ph.ecog, data = lung)Coefficients:(Intercept) sex age ph.ecog6.27343525 0.40109054 -0.00747544 -0.33963810Scale= 0.731109Loglik(model)= -1132.4 Loglik(intercept only)= -1147.4Chisq= 29.98 on 3 degrees of freedom, p= 1.4e-06n=227 (1 observation deleted due to missingness)>>> 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.551290 -7.615741 -7.424565 -7.533378 -7.571272 -14.527489Scale= 0.755741Loglik(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.88712062 -0.01360829Scale= 0.7587515Loglik(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, 12 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)1.478892e-09>> 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: 9n= 19>> # 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, 7 Newton-RaphsonVariance of random effect= 0.1 I-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.066 119.7 1 0.00000frailty(grp, df = 3, dist 16.4 3 0.00096Scale= 0.251Iterations: 8 outer, 33 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> #> data(ovarian)> ofit1 <- survreg(Surv(futime, fustat) ~ age + ridge(ecog.ps, rx), ovarian)>> predict(ofit1)[1] 207.7548 172.7986 358.7735 1426.6498 1353.7357 843.8627 1102.1652[8] 859.5084 416.3290 1280.4094 820.7329 1882.7269 876.1267 1041.8963[15] 3477.0615 2622.9894 3761.5364 2207.8901 1362.2021 3113.9802 879.2010[22] 180.8417 2501.0794 645.2434 555.8305 936.0103> predict(ofit1, type='response')[1] 207.7548 172.7986 358.7735 1426.6498 1353.7357 843.8627 1102.1652[8] 859.5084 416.3290 1280.4094 820.7329 1882.7269 876.1267 1041.8963[15] 3477.0615 2622.9894 3761.5364 2207.8901 1362.2021 3113.9802 879.2010[22] 180.8417 2501.0794 645.2434 555.8305 936.0103> predict(ofit1, type='terms', se=T)$fitage ridge(ecog.ps, rx)1 -1.37776207 -0.17654822 -1.56199426 -0.17654823 -0.87785423 -0.13012354 0.23872053 0.13369475 0.49650242 -0.17654826 -0.02255561 -0.13012357 -0.06575647 0.18011938 -0.31442775 0.18011939 -0.68264498 -0.176548210 0.08414685 0.180119311 -0.05034768 -0.130123512 0.51611283 0.133694713 -0.29527755 0.180119314 -0.07556594 0.133694715 1.43982204 -0.176548216 1.11152449 -0.130123517 1.47203732 -0.130123518 0.98567183 -0.176548219 0.19249416 0.133694720 1.01929333 0.133694721 -0.29177478 0.180119322 -1.56292322 -0.130123523 1.11035690 -0.176548224 -0.60116078 0.180119325 -0.70390024 0.133694726 -0.18273714 0.1336947attr(,"constant")(Intercept)11.14419$se.fitage ridge(ecog.ps, rx)1 0.0303419941 0.041195402 0.0343992782 0.041195403 0.0193326906 0.041760164 0.0052572625 0.041597205 0.0109343070 0.041195406 0.0004967347 0.041760167 0.0014481328 0.041292878 0.0069245373 0.041292879 0.0150336625 0.0411954010 0.0018531379 0.0412928711 0.0011087902 0.0417601612 0.0113661806 0.0415972013 0.0065027990 0.0412928714 0.0016641635 0.0415972015 0.0317087202 0.0411954016 0.0244787327 0.0417601617 0.0324181865 0.0417601618 0.0217071217 0.0411954019 0.0042392346 0.0415972020 0.0224475568 0.0415972021 0.0064256586 0.0412928722 0.0344197363 0.0417601623 0.0244530192 0.0411954024 0.0132391631 0.0412928725 0.0155017601 0.0415972026 0.0040243590 0.04159720>> temp1 <- predict(ofit1,type="link", se=T)> temp2 <- predict(ofit1, type= 'response', se=T)> all.equal(temp2$se.fit, temp1$se.fit* exp(temp1$fit))[1] TRUE> #> # 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 <- survreg(Surv(time, status) ~ cage + I(cage^2), stanford2,+ dist='lognormal')> fit1Call:survreg(formula = Surv(time, status) ~ cage + I(cage^2), data = stanford2,dist = "lognormal")Coefficients:(Intercept) cage I(cage^2)6.717591081 -0.061908619 -0.003504315Scale= 2.362872Loglik(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 1201.379202e-01 1.452449e-01 2.628088e-02 7.320182e-02 7.624325e-02 3.994793e-0299 108 179 43 134 1606.328460e-02 6.128977e-02 9.685668e-03 4.767550e-02 2.980553e-02 1.036051e-01177 153 136 133 176 668.990601e-03 2.114950e-02 2.557694e-02 1.591463e-01 8.618405e-03 3.389342e-02157 114 46 65 184 881.141319e-02 1.990887e-02 2.044977e-02 2.480540e-02 1.085737e-05 5.474389e-02182 180 163 84 90 681.786494e-03 2.574818e-03 7.654075e-03 2.024457e-02 8.561193e-02 3.894985e-0248 174 151 125 73 1057.007563e-02 3.767424e-03 8.314677e-03 1.248554e-02 1.954896e-02 1.831982e-02117 96 39 38 106 141.739299e-02 1.789440e-02 2.406180e-02 2.364315e-02 4.717184e-02 2.051895e-02123 135 111 83 143 694.763900e-02 1.663805e-02 1.367017e-02 3.204509e-02 1.857902e-02 2.058865e-0227 113 167 156 141 303.896724e-02 3.775024e-02 5.091513e-03 1.528402e-02 8.682136e-03 1.746134e-02144 158 79 102 77 362.593291e-02 6.620379e-03 1.375918e-02 1.547851e-02 1.786268e-02 2.330671e-02183 122 162 121 87 23.720939e-05 1.696469e-02 5.954816e-03 1.233287e-02 1.655940e-02 1.089489e-0164 150 85 71 19 216.015390e-02 7.469432e-03 1.666501e-02 1.893415e-02 2.645489e-02 1.843298e-01175 169 148 138 98 1041.789942e-02 4.379957e-03 7.619698e-03 9.332609e-03 1.428798e-02 1.445962e-02103 12 89 3 100 551.449500e-02 3.404298e-02 3.358405e-02 3.113308e-02 1.412657e-02 1.179741e-02142 63 168 72 137 108.641594e-03 1.426955e-02 4.554043e-03 1.094162e-02 9.645966e-03 1.226564e-02124 17 94 82 170 1491.222512e-02 1.088510e-02 1.493684e-02 1.844221e-02 3.988063e-02 3.038321e-0242 128 67 109 75 262.127744e-02 1.439502e-02 1.285836e-02 8.944979e-03 1.997791e-02 2.757124e-0297 58 178 140 32 1262.549339e-02 2.356049e-02 2.057505e-03 1.269584e-02 1.103392e-02 1.253031e-0251 101 29 33 164 601.430228e-02 1.637415e-02 2.201027e-02 1.118993e-02 6.417566e-03 8.492275e-03152 145 112 76 47 1188.651522e-03 9.608673e-03 1.609215e-02 2.168279e-02 2.622512e-02 2.274275e-025 129 31 35 40 1301.184996e-02 9.391147e-03 8.772106e-03 8.526052e-03 9.451658e-03 1.295998e-0228 56 91 44 23 371.285986e-02 1.536640e-02 2.031499e-02 2.807956e-02 1.965943e-02 1.733256e-0270 132 9 81 59 1279.129009e-03 9.121716e-03 9.083024e-03 1.025238e-02 1.032186e-02 1.183693e-02131 80 20 25 165 241.403298e-02 2.363944e-02 2.181249e-02 2.723391e-02 2.043510e-02 2.019541e-02172 146 86 107 95 1161.152649e-02 1.265907e-02 1.538526e-02 2.107503e-02 2.298470e-02 2.128394e-0241 61 155 166 154 41.791007e-02 1.763098e-02 1.345059e-02 1.285115e-02 1.218090e-02 1.470506e-0292 93 62 34 15 1732.599207e-02 3.098462e-02 3.037745e-02 2.166521e-02 1.478524e-02 7.517967e-03171 52 110 50 45 538.681586e-03 1.679632e-02 2.540017e-02 3.470668e-02 3.229506e-02 3.017736e-0254 147 115 16 1 62.416303e-02 1.870027e-02 2.172488e-02 1.164271e-01 4.257797e-02 2.459122e-027 57 78 161 11 83.585527e-02 3.587688e-02 2.865161e-02 2.603297e-02 5.640968e-02 4.338250e-0249 13 22 183.425475e-02 6.262791e-02 1.029315e-01 1.442429e-01> # 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)++ plot(ldcase, xlab="Case Number", ylab="(1/2) A")+ title (main="Case weight pertubations")+ plot(ldresp, xlab="Case Number", ylab="(1/2) A")+ title(main="Response pertubations")+ }>> plot1()> #> # 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] TRUE>> 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] TRUE>> aeq(p4$fit, p5$fit)[1] TRUE> #!aeq(p3$se.fit, p4$se.fit) #this one should be false> aeq(p4$se.fit, p5$se.fit) #this one true[1] TRUE>> #> # 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)> #all.equal(fit1$loglik, fit2$loglik, tolerance= tol)> all.equal(fit1$coef, fit2$coef)[1] TRUE> all.equal(fit1$loglik, fit2$loglik)[1] TRUE>> # 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)> ##all.equal(fit1$loglik[2], fit2$loglik[2], tolerance=tol)> all.equal(fit1$coef, fit2$coef)[1] TRUE> all.equal(fit1$loglik[2], fit2$loglik[2])[1] TRUE>> # 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)> #all.equal(fit1$loglik[2], fit2$loglik[2], tolerance=tol)> all.equal(fit1$coef, fit2$coef)[1] "Mean relative difference: 0.0002457368"> all.equal(fit1$loglik[2], fit2$loglik[2])[1] "Mean relative difference: 4.971155e-07">> 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.062171Scale:sex=1 sex=20.8167551 0.6533036Loglik(model)= -1152.5 Loglik(intercept only)= -1152.5n= 228> fit2Call:survreg(formula = Surv(time, status) ~ strata(sex) + sex, data = lung)Coefficients:(Intercept) sex5.494409 0.380171Scale:sex=1 sex=20.8084294 0.6355816Loglik(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] TRUE> aeq(fit2$loglik[2], (fit3a$loglik + fit3b$loglik)[2], tolerance=tol)[1] TRUE> aeq(fit2$coef[1] + 1:2*fit2$coef[2], c(fit3a$coef, fit3b$coef), tolerance=tol)[1] TRUE>> #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, 4 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.800 0.636Iterations: 1 outer, 5 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.001270879"> aeq(fit3$loglik[2], (fit3a$loglik + fit3b$loglik)[2], tolerance=tol*10)[1] "Mean relative difference: 8.673582e-05"> pred <- predict(fit3)> aeq(pred[lung$sex==1] , predict(fit3a), tolerance=tol*10)[1] TRUE> aeq(pred[lung$sex==2], predict(fit3b), tolerance=tol*10)###????FIXME[1] "Mean relative difference: 0.01158253">>>>> #> # 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.9831358 -0.2385112Scale= 0.2637875Loglik(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.983136 -0.2385112 0.2637875rfit1 4.902438 -0.2189410 0.2025434>>> 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.1346218 0.1346218 0.1224049 0.1207863> apply(temp,2,mean)surv surv.sparse cox cox.sparse6.546887e-19 1.242665e-18 -1.301043e-17 -2.279128e-18>> # 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)> options(na.action="na.exclude")> 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)> index<-index[!is.na(yy[index,1])]> matplot(mpip$ef[index], yy[index,], type='l', lty=c(1,2,2), col=1,+ 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)> 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,+ 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,+ xlab="Ejection Fraction", ylab="Log-normal model predictor",+ main="Post-Infarction Survival")> q()