The R Project SVN
R-packages
(root)/trunk/Old/survival/Examples/testfrail/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 .DataWarning: Cannot open audit file> #> # Set up for the test> #> dyn.load("../loadmod.o")> attach("../.Data")> options(na.action="na.omit", contrasts="contr.treatment")> #> # Test the logic of the new program, by fitting some no-frailty models> # (theta=0). It should give exactly the same answers as 'ordinary' coxph.> # By default frailty models run with eps=1e-7, ordinary with 1e-4. I match> # these to get the same number of iterations.> #> test1 <- data.frame(time= c(4, 3,1,1,2,2,3),+ status=c(1,NA,1,0,1,1,0),+ x= c(0, 2,1,1,1,0,0))> > test2 <- data.frame(start=c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8),+ stop =c(2, 3, 6, 7, 8, 9, 9, 9,14,17),+ event=c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),+ x =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0) )> > zz <- rep(0, nrow(test1))> tfit1 <- coxph(Surv(time,status) ~x, test1, eps=1e-7)> tfit2 <- coxph(Surv(time,status) ~x + frailty(zz, theta=0, sparse=T), test1)> tfit3 <- coxph(Surv(zz,time,status) ~x + frailty(zz, theta=0,sparse=T), test1)> > temp <- c("coefficients", "var", "loglik", "linear.predictors",+ "means", "n")> > all.equal(tfit1[temp], tfit2[temp])[1] T> all.equal(tfit1[temp], tfit3[temp])[1] T> > zz <- rep(0, nrow(test2))> tfit1 <- coxph(Surv(start, stop, event) ~x, test2, eps=1e-7)> tfit2 <- coxph(Surv(start, stop, event) ~ x + frailty(zz, theta=0, sparse=T),+ test2)> all.equal(tfit1[temp], tfit2[temp])[1] T> > > > # Tests using the rats data> #> # (Female rats, from Mantel et al, Cancer Research 37,> # 3863-3868, November 77)>> rats <- read.table("data.rats", col.names=c("litter", "rx", "time",+ "status"))> > rfit <- coxph(Surv(time,status) ~ rx + frailty(litter), rats,+ method="breslow")> names(rfit)[1] "coefficients" "var" "var2"[4] "loglik" "iter" "linear.predictors"[7] "residuals" "means" "method"[10] "frail" "fvar" "df"[13] "df2" "penalty" "pterms"[16] "assign2" "history" "coxlist1"[19] "printfun" "n" "terms"[22] "assign" "wald.test" "y"[25] "formula" "call"> rfitCall:coxph(formula = Surv(time, status) ~ rx + frailty(litter), data = rats, method= "breslow")coef se(coef) se2 Chisq DF prx 0.906 0.323 0.319 7.88 1.0 0.005frailty(litter) 16.89 13.8 0.250Iterations: 6 outer, 20 Newton-RaphsonVariance of random effect= 0.474 EM likelihood = -181.1Degrees of freedom for terms= 1.0 13.9Likelihood ratio test=36.3 on 14.83 df, p=0.00145 n= 150> > rfit$iter[1] 6 20> rfit$df[1] 0.9759431 13.8548423> rfit$history[[1]]$theta:[1] 0.4742848$done:c.loglikT$history:theta loglik c.loglik[1,] 0.0000000 -181.8451 -181.8451[2,] 1.0000000 -168.3683 -181.5458[3,] 0.5000000 -173.3117 -181.0788[4,] 0.3090061 -175.9446 -181.1490[5,] 0.4645720 -173.7590 -181.0775[6,] 0.4736209 -173.6431 -181.0773$c.loglik:[1] -181.0773> > rfit1 <- coxph(Surv(time,status) ~ rx + frailty(litter, theta=1), rats,+ method="breslow")> rfit1Call:coxph(formula = Surv(time, status) ~ rx + frailty(litter, theta = 1), data =rats, method = "breslow")coef se(coef) se2 Chisq DF prx 0.918 0.327 0.321 7.85 1.0 0.0051frailty(litter, theta = 1 27.25 22.7 0.2300Iterations: 1 outer, 5 Newton-RaphsonVariance of random effect= 1 EM likelihood = -181.5Degrees of freedom for terms= 1.0 22.7Likelihood ratio test=50.7 on 23.67 df, p=0.001 n= 150> > rfit2 <- coxph(Surv(time,status) ~ frailty(litter), rats)> rfit2Call:coxph(formula = Surv(time, status) ~ frailty(litter), data = rats)coef se(coef) se2 Chisq DF pfrailty(litter) 18 14.6 0.24Iterations: 6 outer, 17 Newton-RaphsonVariance of random effect= 0.504 EM likelihood = -184.8Degrees of freedom for terms= 14.6Likelihood ratio test=30 on 14.61 df, p=0.0101 n= 150> #> # Here is a test case with multiple smoothing terms> #>> fit0 <- coxph(Surv(time, status) ~ ph.ecog + age, lung)> fit1 <- coxph(Surv(time, status) ~ ph.ecog + pspline(age,3), lung)> fit2 <- coxph(Surv(time, status) ~ ph.ecog + pspline(age,4), lung)> fit3 <- coxph(Surv(time, status) ~ ph.ecog + pspline(age,8), lung)> > > > fit4 <- coxph(Surv(time, status) ~ ph.ecog + pspline(wt.loss,3), lung)> > fit5 <-coxph(Surv(time, status) ~ ph.ecog + pspline(age,3) ++ pspline(wt.loss,3), lung)> > fit1Call:coxph(formula = Surv(time, status) ~ ph.ecog + pspline(age, 3), data = lung)coef se(coef) se2 Chisq DF pph.ecog 0.4480 0.11707 0.11678 14.64 1.00 0.00013pspline(age, 3), linear 0.0113 0.00928 0.00928 1.47 1.00 0.22000pspline(age, 3), nonlin 2.08 2.08 0.37000Iterations: 4 outer, 10 Newton-RaphsonTheta= 0.861Degrees of freedom for terms= 1.0 3.1Likelihood ratio test=21.9 on 4.08 df, p=0.000227n=227 (1 observations deleted due to missing)> fit2Call:coxph(formula = Surv(time, status) ~ ph.ecog + pspline(age, 4), data = lung)coef se(coef) se2 Chisq DF pph.ecog 0.4505 0.11766 0.11723 14.66 1.00 0.00013pspline(age, 4), linear 0.0112 0.00927 0.00927 1.45 1.00 0.23000pspline(age, 4), nonlin 2.96 3.08 0.41000Iterations: 4 outer, 10 Newton-RaphsonTheta= 0.797Degrees of freedom for terms= 1.0 4.1Likelihood ratio test=22.7 on 5.07 df, p=0.000412n=227 (1 observations deleted due to missing)> fit3Call:coxph(formula = Surv(time, status) ~ ph.ecog + pspline(age, 8), data = lung)coef se(coef) se2 Chisq DF pph.ecog 0.4764 0.12024 0.11925 15.70 1.00 7.4e-05pspline(age, 8), linear 0.0117 0.00923 0.00923 1.61 1.00 2.0e-01pspline(age, 8), nonlin 6.93 6.99 4.3e-01Iterations: 5 outer, 13 Newton-RaphsonTheta= 0.691Degrees of freedom for terms= 1 8Likelihood ratio test=27.6 on 8.97 df, p=0.00108n=227 (1 observations deleted due to missing)> fit4Call:coxph(formula = Surv(time, status) ~ ph.ecog + pspline(wt.loss, 3), data = lung)coef se(coef) se2 Chisq DF pph.ecog 0.51545 0.12960 0.12737 15.82 1.00 0.00007pspline(wt.loss, 3), line -0.00702 0.00655 0.00655 1.15 1.00 0.28000pspline(wt.loss, 3), nonl 2.45 2.09 0.31000Iterations: 3 outer, 8 Newton-RaphsonTheta= 0.776Degrees of freedom for terms= 1.0 3.1Likelihood ratio test=21.1 on 4.06 df, p=0.000326n=213 (15 observations deleted due to missing)> fit5Call:coxph(formula = Surv(time, status) ~ ph.ecog + pspline(age, 3) + pspline(wt.loss, 3), data = lung)coef se(coef) se2 Chisq DF pph.ecog 0.47422 0.13495 0.13206 12.35 1.00 0.00044pspline(age, 3), linear 0.01368 0.00976 0.00974 1.96 1.00 0.16000pspline(age, 3), nonlin 1.90 2.07 0.40000pspline(wt.loss, 3), line -0.00717 0.00661 0.00660 1.18 1.00 0.28000pspline(wt.loss, 3), nonl 2.08 2.03 0.36000Iterations: 4 outer, 10 Newton-RaphsonTheta= 0.85Theta= 0.779Degrees of freedom for terms= 1.0 3.1 3.0Likelihood ratio test=25.2 on 7.06 df, p=0.000726n=213 (15 observations deleted due to missing)> > rm(fit1, fit2, fit3, fit4, fit5)> #> # Test on the ovarian data>> fit1 <- coxph(Surv(futime, fustat) ~ rx + age, ovarian)> fit2 <- coxph(Surv(futime, fustat) ~ rx + pspline(age, df=2),+ data=ovarian)> fit2$iter[1] 2 7> > fit2$df[1] 0.9426611 1.9293052> > fit2$history$"pspline(age, df = 2)":$"pspline(age, df = 2)"$theta:[1] 0.4468868$"pspline(age, df = 2)"$done:[1] T$"pspline(age, df = 2)"$history:thetas dfs[1,] 1.0000000 1.000000[2,] 0.0000000 5.000000[3,] 0.6000000 1.734267[4,] 0.4845205 1.929305$"pspline(age, df = 2)"$half:[1] 0> > fit4 <- coxph(Surv(futime, fustat) ~ rx + pspline(age, df=4),+ data=ovarian)> fit4Call:coxph(formula = Surv(futime, fustat) ~ rx + pspline(age, df = 4), data =ovarian)coef se(coef) se2 Chisq DF prx -0.373 0.761 0.749 0.24 1.00 0.6200pspline(age, df = 4), lin 0.139 0.044 0.044 9.98 1.00 0.0016pspline(age, df = 4), non 2.59 2.93 0.4500Iterations: 3 outer, 13 Newton-RaphsonTheta= 0.242Degrees of freedom for terms= 1.0 3.9Likelihood ratio test=19.4 on 4.9 df, p=0.00149 n= 26> > > # Simulation for the ovarian data set> #> fit1 <- coxph(Surv(futime, fustat) ~ rx + ridge(age, ecog.ps, theta=1),+ ovarian)> > dfs <- eigen(solve(fit1$var, fit1$var2))$values> > temp <- matrix(rnorm(30000), ncol=3)> temp2 <- apply((temp^2) %*% dfs, 1, sum)> > round(rbind(quantile(temp2, c(.8, .9, .95, .99)),+ qchisq( c(.8, .9, .95, .99), sum(fit1$df))), 3)80% 90% 95% 99%[1,] 4.287 5.779 7.156 10.458[2,] 4.313 5.874 7.399 10.861> # From: McGilchrist and Aisbett, Biometrics 47, 461-66, 1991> # Data on the recurrence times to infection, at the point of insertion of> # the catheter, for kidney patients using portable dialysis equipment.> # Catheters may be removed for reasons other than infection, in which case> # the observation is censored. Each patient has exactly 2 observations.>> # Variables: patient, time, status, age,> # sex (1=male, 2=female),> # disease type (0=GN, 1=AN, 2=PKD, 3=Other)> # author's estimate of the frailty>> # I don't match their answers, and I think that I'm right>> kidney <- read.table("data.kidney", col.names=c("id", "time", "status",+ "age", "sex", "disease", "frail"))> kidney$disease <- factor(kidney$disease, levels=c(3, 0:2),+ labels=c("Other", "GN", "AN", "PKD"))> > kfit <- coxph(Surv(time, status)~ age + sex + disease + frailty(id), kidney)> kfit1<- coxph(Surv(time, status) ~age + sex + disease ++ frailty(id, theta=1), kidney, iter=20)> kfit0 <- coxph(Surv(time, status)~ age + sex + disease, kidney)> temp <- coxph(Surv(time, status) ~age + sex + disease ++ frailty(id, theta=1, sparse=F), kidney)> > > # Check out the EM based score equations> # temp1 and kfit1 should have essentially the same coefficients> # temp2 should equal kfit1$frail> # equality won't be exact because of the different iteration paths> temp1 <- coxph(Surv(time, status) ~ age + sex + disease ++ offset(kfit1$frail[id]), kidney)> rr <- tapply(resid(temp1), kidney$id, sum)> temp2 <- log(rr/1 +1)> all.equal(temp1$coef, kfit1$coef)[1] "Mean relative difference: 9.243907e-06"> all.equal(temp2, kfit1$frail)[1] "names for target but not for current"[2] "Mean relative difference: 0.002409231"> > > > kfitCall:coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id), data =kidney)coef se(coef) se2 Chisq DF page 0.00318 0.0111 0.0111 0.08 1 7.8e-01sex -1.48314 0.3582 0.3582 17.14 1 3.5e-05diseaseGN 0.08796 0.4064 0.4064 0.05 1 8.3e-01diseaseAN 0.35079 0.3997 0.3997 0.77 1 3.8e-01diseasePKD -1.43111 0.6311 0.6311 5.14 1 2.3e-02frailty(id) 0.00 0 9.3e-01Iterations: 6 outer, 28 Newton-RaphsonVariance of random effect= 5e-07 EM likelihood = -179.1Degrees of freedom for terms= 1 1 3 0Likelihood ratio test=17.6 on 5 df, p=0.00342 n= 76> kfit1Call:coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id, theta =1), data = kidney, iter.max = 20)coef se(coef) se2 Chisq DF page 0.00389 0.0196 0.00943 0.04 1.0 0.84000sex -2.00788 0.5910 0.41061 11.54 1.0 0.00068diseaseGN 0.35334 0.7165 0.38015 0.24 1.0 0.62000diseaseAN 0.52363 0.7229 0.40462 0.52 1.0 0.47000diseasePKD -0.45980 1.0898 0.66091 0.18 1.0 0.67000frailty(id, theta = 1) 28.48 18.8 0.06900Iterations: 1 outer, 10 Newton-RaphsonVariance of random effect= 1 EM likelihood = -182.5Degrees of freedom for terms= 0.2 0.5 1.1 18.8Likelihood ratio test=63.8 on 20.55 df, p=2.55e-06 n= 76> kfit0Call:coxph(formula = Surv(time, status) ~ age + sex + disease, data = kidney)coef exp(coef) se(coef) z page 0.00318 1.003 0.0111 0.285 7.8e-01sex -1.48319 0.227 0.3582 -4.140 3.5e-05diseaseGN 0.08796 1.092 0.4064 0.216 8.3e-01diseaseAN 0.35078 1.420 0.3997 0.878 3.8e-01diseasePKD -1.43111 0.239 0.6311 -2.268 2.3e-02Likelihood ratio test=17.6 on 5 df, p=0.00342 n= 76> tempCall:coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id, theta =1, sparse = F), data = kidney)coef se(coef) se2 Chisq DF page 0.00389 0.0186 0.0112 0.04 1.0 0.83000sex -2.00763 0.5762 0.4080 12.14 1.0 0.00049diseaseGN 0.35335 0.6786 0.4315 0.27 1.0 0.60000diseaseAN 0.52340 0.6891 0.4404 0.58 1.0 0.45000diseasePKD -0.45934 1.0139 0.7130 0.21 1.0 0.65000frailty(id, theta = 1, sp 26.23 18.7 0.12000Iterations: 1 outer, 5 Newton-RaphsonVariance of random effect= 1 EM likelihood = -182.5Degrees of freedom for terms= 0.4 0.5 1.4 18.7Likelihood ratio test=63.8 on 20.97 df, p=3.27e-06 n= 76> > #> # Now fit the data using REML> #> kfitm1 <- coxph(Surv(time,status) ~ age + sex + disease ++ frailty(id, dist="gauss"), kidney)> kfitm2 <- coxph(Surv(time,status) ~ age + sex + disease ++ frailty(id, dist="gauss", sparse=F), kidney)> kfitm1Call:coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id, dist ="gauss"), data = kidney)coef se(coef) se2 Chisq DF page 0.00489 0.015 0.0106 0.11 1.0 0.74000sex -1.69703 0.461 0.3617 13.56 1.0 0.00023diseaseGN 0.17980 0.545 0.3927 0.11 1.0 0.74000diseaseAN 0.39283 0.545 0.3982 0.52 1.0 0.47000diseasePKD -1.13630 0.825 0.6173 1.90 1.0 0.17000frailty(id, dist = "gauss 17.89 12.1 0.12000Iterations: 6 outer, 30 Newton-RaphsonVariance of random effect= 0.493Degrees of freedom for terms= 0.5 0.6 1.7 12.1Likelihood ratio test=47.5 on 14.89 df, p=2.82e-05 n= 76> summary(kfitm2)Call:coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id, dist ="gauss", sparse = F), data = kidney)n= 76coef se(coef) se2 Chisq DF page 0.00492 0.0149 0.0108 0.11 1.0 0.74000sex -1.70204 0.4631 0.3613 13.51 1.0 0.00024diseaseGN 0.18173 0.5413 0.4017 0.11 1.0 0.74000diseaseAN 0.39442 0.5428 0.4052 0.53 1.0 0.47000diseasePKD -1.13160 0.8175 0.6298 1.92 1.0 0.17000frailty(id, dist = "gauss 18.13 12.3 0.12000exp(coef) exp(-coef) lower .95 upper .95age 1.005 0.995 0.9760 1.035sex 0.182 5.485 0.0736 0.452diseaseGN 1.199 0.834 0.4151 3.465diseaseAN 1.484 0.674 0.5120 4.299diseasePKD 0.323 3.101 0.0650 1.601gamma:1 1.701 0.588 0.5181 5.586gamma:2 1.424 0.702 0.3851 5.266gamma:3 1.159 0.863 0.3828 3.511gamma:4 0.623 1.606 0.2340 1.657gamma:5 1.254 0.797 0.3981 3.953gamma:6 1.135 0.881 0.3834 3.360gamma:7 1.973 0.507 0.5694 6.834gamma:8 0.620 1.614 0.2166 1.772gamma:9 0.823 1.215 0.2888 2.346gamma:10 0.503 1.988 0.1747 1.448gamma:11 0.757 1.322 0.2708 2.113gamma:12 1.105 0.905 0.3343 3.651gamma:13 1.302 0.768 0.4275 3.967gamma:14 0.591 1.691 0.1854 1.885gamma:15 0.545 1.835 0.1858 1.598gamma:16 1.044 0.958 0.3142 3.470gamma:17 0.914 1.095 0.3000 2.782gamma:18 0.918 1.089 0.3248 2.597gamma:19 0.643 1.556 0.1951 2.117gamma:20 1.170 0.855 0.3453 3.963gamma:21 0.334 2.997 0.1020 1.091gamma:22 0.687 1.455 0.2353 2.006gamma:23 1.478 0.677 0.4756 4.592gamma:24 1.017 0.983 0.3156 3.278gamma:25 0.810 1.235 0.2749 2.384gamma:26 0.614 1.627 0.2149 1.757gamma:27 1.088 0.919 0.3282 3.610gamma:28 1.542 0.649 0.4923 4.829gamma:29 1.379 0.725 0.4377 4.342gamma:30 1.375 0.727 0.4444 4.253gamma:31 1.445 0.692 0.4703 4.438gamma:32 1.199 0.834 0.3521 4.085gamma:33 1.945 0.514 0.5523 6.849gamma:34 0.862 1.161 0.2769 2.682gamma:35 1.703 0.587 0.5266 5.508gamma:36 0.827 1.209 0.2281 3.002gamma:37 1.471 0.680 0.3894 5.555gamma:38 1.048 0.954 0.3068 3.579Iterations: 6 outer, 17 Newton-RaphsonVariance of random effect= 0.509Degrees of freedom for terms= 0.5 0.6 1.7 12.3Rsquare= 0.788 (max possible= 0.997 )Likelihood ratio test= 118 on 15.14 df, p=0Wald test = 37.4 on 15.14 df, p=0.00119> #> # Fit the kidney data using AIC> #>> # gamma, corrected aic> coxph(Surv(time, status) ~ age + sex + frailty(id, method="aic", caic=T),+ kidney)Call:coxph(formula = Surv(time, status) ~ age + sex + frailty(id, method = "aic",caic = T), data = kidney)coef se(coef) se2 Chisq DF page 0.00364 0.0105 0.00891 0.12 1.0 0.73000sex -1.31907 0.3955 0.32493 11.13 1.0 0.00085frailty(id, method = "aic 13.54 7.8 0.08700Iterations: 9 outer, 47 Newton-RaphsonVariance of random effect= 0.202 EM likelihood = -182.1Degrees of freedom for terms= 0.7 0.7 7.8Likelihood ratio test=33.3 on 9.2 df, p=0.000137 n= 76> > coxph(Surv(time, status) ~ age + sex + frailty(id, dist="t"), kidney)Call:coxph(formula = Surv(time, status) ~ age + sex + frailty(id, dist = "t"), data= kidney)coef se(coef) se2 Chisq DF page 0.00558 0.012 0.00873 0.22 1.0 0.6400sex -1.65036 0.481 0.38545 11.77 1.0 0.0006frailty(id, dist = "t") 20.05 13.8 0.1200Iterations: 9 outer, 44 Newton-RaphsonVariance of random effect= 0.807Degrees of freedom for terms= 0.5 0.6 13.8Likelihood ratio test=48.2 on 14.93 df, p=2.24e-05 n= 76> coxph(Surv(time, status) ~ age + sex + frailty(id, dist="gauss", method="aic",+ caic=T), kidney)Call:coxph(formula = Surv(time, status) ~ age + sex + frailty(id, dist = "gauss",method = "aic", caic = T), data = kidney)coef se(coef) se2 Chisq DF page 0.00303 0.0103 0.00895 0.09 1.00 0.7700sex -1.15153 0.3637 0.30556 10.03 1.00 0.0015frailty(id, dist = "gauss 12.36 6.76 0.0800Iterations: 7 outer, 30 Newton-RaphsonVariance of random effect= 0.185Degrees of freedom for terms= 0.8 0.7 6.8Likelihood ratio test=28.4 on 8.22 df, p=0.000476 n= 76> > > # uncorrected aic> coxph(Surv(time, status) ~ age + sex + frailty(id, method="aic", caic=F),+ kidney)Call:coxph(formula = Surv(time, status) ~ age + sex + frailty(id, method = "aic",caic = F), data = kidney)coef se(coef) se2 Chisq DF page 0.00758 0.0146 0.00836 0.27 1.0 0.60000sex -1.86230 0.5503 0.39401 11.45 1.0 0.00071frailty(id, method = "aic 35.99 19.1 0.01100Iterations: 10 outer, 74 Newton-RaphsonVariance of random effect= 0.824 EM likelihood = -182.6Degrees of freedom for terms= 0.3 0.5 19.1Likelihood ratio test=60 on 19.91 df, p=6.83e-06 n= 76Warning messages:Inner loop failed to coverge for iterations 4 in: coxpenal.fit(X, Y, strats,offset, init = init, iter.max = iter.max, outer.max ....> > coxph(Surv(time, status) ~ age + sex + frailty(id, dist="t", caic=F), kidney)Call:coxph(formula = Surv(time, status) ~ age + sex + frailty(id, dist = "t", caic= F), data = kidney)coef se(coef) se2 Chisq DF page 0.00558 0.012 0.00873 0.22 1.0 0.6400sex -1.65036 0.481 0.38545 11.77 1.0 0.0006frailty(id, dist = "t", c 20.05 13.8 0.1200Iterations: 9 outer, 44 Newton-RaphsonVariance of random effect= 0.807Degrees of freedom for terms= 0.5 0.6 13.8Likelihood ratio test=48.2 on 14.93 df, p=2.24e-05 n= 76> #temp <- sas.get("../../../../data/moertel/sasdata", "anal")> #colon <- temp[temp$study==1,]> #rm(temp)> #colon$rx <- factor(colon$rx, levels=1:3, labels=c("Obs", "Lev", "Lev+5FU"))>> data.restore("data.colon")[1] "data.colon"> #> # Fit models to the Colon cancer data used in Lin> #> fitc1 <- coxph(Surv(time, status) ~ rx + extent + node4 + cluster(id)+ + strata(etype), colon)> fitc1Call:coxph(formula = Surv(time, status) ~ rx + extent + node4 + cluster(id) + strata(etype), data = colon)coef exp(coef) se(coef) robust se z prxLev -0.0362 0.964 0.0768 0.1056 -0.343 7.3e-01rxLev+5FU -0.4488 0.638 0.0840 0.1168 -3.842 1.2e-04extent 0.5155 1.674 0.0796 0.1097 4.701 2.6e-06node4 0.8799 2.411 0.0681 0.0961 9.160 0.0e+00Likelihood ratio test=248 on 4 df, p=0 n= 1858> > fitc2 <- coxph(Surv(time, status) ~ rx + extent + node4 ++ frailty(id, dist="gauss", trace=T)+ + strata(etype), colon)theta resid fsum trace[1,] 1 0.5721865 677.2472 498.2323[2,] 3 0.8244916 2430.3958 880.5538new theta= 6theta resid fsum trace[1,] 1 0.5721865 677.2472 498.2323[2,] 3 0.8244916 2430.3958 880.5538[3,] 6 0.3152272 4520.0041 1279.6138new theta= 12theta resid fsum trace[1,] 1 0.5721865 677.2472 498.2323[2,] 3 0.8244916 2430.3958 880.5538[3,] 6 0.3152272 4520.0041 1279.6138[4,] 12 -2.1486199 7550.5646 1950.6313new theta= 7.554873theta resid fsum trace[1,] 1.000000 0.5721865 677.2472 498.2323[2,] 3.000000 0.8244916 2430.3958 880.5538[3,] 6.000000 0.3152272 4520.0041 1279.6138[4,] 12.000000 -2.1486199 7550.5646 1950.6313[5,] 7.554873 -0.1827268 5420.8778 1463.2371new theta= 7.004443theta resid fsum trace[1,] 1.000000 0.57218652 677.2472 498.2323[2,] 3.000000 0.82449159 2430.3958 880.5538[3,] 6.000000 0.31522725 4520.0041 1279.6138[4,] 12.000000 -2.14861992 7550.5646 1950.6313[5,] 7.554873 -0.18272677 5420.8778 1463.2371[6,] 7.004443 0.02102504 5123.0634 1399.3956new theta= 7.06674theta resid fsum trace[1,] 1.000000 0.57218652 677.2472 498.2323[2,] 3.000000 0.82449159 2430.3958 880.5538[3,] 6.000000 0.31522725 4520.0041 1279.6138[4,] 12.000000 -2.14861992 7550.5646 1950.6313[5,] 7.554873 -0.18272677 5420.8778 1463.2371[6,] 7.004443 0.02102504 5123.0634 1399.3956[7,] 7.066740 -0.01293658 5148.9076 1406.6504new theta= 7.04162theta resid fsum trace[1,] 1.000000 0.572186518 677.2472 498.2323[2,] 3.000000 0.824491593 2430.3958 880.5538[3,] 6.000000 0.315227245 4520.0041 1279.6138[4,] 12.000000 -2.148619920 7550.5646 1950.6313[5,] 7.554873 -0.182726773 5420.8778 1463.2371[6,] 7.004443 0.021025043 5123.0634 1399.3956[7,] 7.066740 -0.012936579 5148.9076 1406.6504[8,] 7.041620 0.003463959 5140.3971 1403.7958new theta= 7.047698theta resid fsum trace[1,] 1.000000 0.572186518 677.2472 498.2323[2,] 3.000000 0.824491593 2430.3958 880.5538[3,] 6.000000 0.315227245 4520.0041 1279.6138[4,] 12.000000 -2.148619920 7550.5646 1950.6313[5,] 7.554873 -0.182726773 5420.8778 1463.2371[6,] 7.004443 0.021025043 5123.0634 1399.3956[7,] 7.066740 -0.012936579 5148.9076 1406.6504[8,] 7.041620 0.003463959 5140.3971 1403.7958[9,] 7.047698 -0.001163083 5142.0167 1404.4463new theta= 7.046096theta resid fsum trace[1,] 1.000000 5.721865e-01 677.2472 498.2323[2,] 3.000000 8.244916e-01 2430.3958 880.5538[3,] 6.000000 3.152272e-01 4520.0041 1279.6138[4,] 12.000000 -2.148620e+00 7550.5646 1950.6313[5,] 7.554873 -1.827268e-01 5420.8778 1463.2371[6,] 7.004443 2.102504e-02 5123.0634 1399.3956[7,] 7.066740 -1.293658e-02 5148.9076 1406.6504[8,] 7.041620 3.463959e-03 5140.3971 1403.7958[9,] 7.047698 -1.163083e-03 5142.0167 1404.4463[10,] 7.046096 -1.820316e-05 5141.5307 1404.2792new theta= 7.046071Warning messages:Inner loop failed to coverge for iterations 4 in: coxpenal.fit(X, Y, strats,offset, init = init, iter.max = iter.max, outer.max ....> fitc2Call:coxph(formula = Surv(time, status) ~ rx + extent + node4 + frailty(id, dist ="gauss", trace = T) + strata(etype), data = colon)coef se(coef) se2 Chisq DF prxLev -0.0267 0.241 0.0824 0.01 1 9.1e-01rxLev+5FU -0.7880 0.243 0.1071 10.50 1 1.2e-03extent 1.1305 0.218 0.1068 26.81 1 2.2e-07node4 2.1266 0.210 0.0984 102.56 1 0.0e+00frailty(id, dist = "gauss 5464.64 730 0.0e+00Iterations: 10 outer, 77 Newton-RaphsonVariance of random effect= 7.05Degrees of freedom for terms= 0.3 0.2 0.2 729.7Likelihood ratio test=3544 on 730.49 df, p=0 n= 1858> > fitc3 <- coxph(Surv(time, status) ~ rx + extent + node4 + frailty(id, trace=T)+ + strata(etype), colon)theta loglik c.loglik[1,] 0 -5846.216 -5846.216[2,] 1 -5305.049 -5590.102new theta= 2theta loglik c.loglik[1,] 0 -5846.216 -5846.216[2,] 1 -5305.049 -5590.102[3,] 2 -5036.927 -5479.479new theta= 4theta loglik c.loglik[1,] 0 -5846.216 -5846.216[2,] 1 -5305.049 -5590.102[3,] 2 -5036.927 -5479.479[4,] 4 -4740.394 -5385.887new theta= 8theta loglik c.loglik[1,] 0 -5846.216 -5846.216[2,] 1 -5305.049 -5590.102[3,] 2 -5036.927 -5479.479[4,] 4 -4740.394 -5385.887[5,] 8 -4457.094 -5347.375new theta= 16theta loglik c.loglik[1,] 0 -5846.216 -5846.216[2,] 1 -5305.049 -5590.102[3,] 2 -5036.927 -5479.479[4,] 4 -4740.394 -5385.887[5,] 8 -4457.094 -5347.375[6,] 16 -4223.785 -5393.362new theta= 8.740343theta loglik c.loglik[1,] 0.000000 -5846.216 -5846.216[2,] 1.000000 -5305.049 -5590.102[3,] 2.000000 -5036.927 -5479.479[4,] 4.000000 -4740.394 -5385.887[5,] 8.000000 -4457.094 -5347.375[6,] 16.000000 -4223.785 -5393.362[7,] 8.740343 -4423.925 -5348.128new theta= 8.058theta loglik c.loglik[1,] 0.000000 -5846.216 -5846.216[2,] 1.000000 -5305.049 -5590.102[3,] 2.000000 -5036.927 -5479.479[4,] 4.000000 -4740.394 -5385.887[5,] 8.000000 -4457.094 -5347.375[6,] 16.000000 -4223.785 -5393.362[7,] 8.740343 -4423.925 -5348.128[8,] 8.058000 -4454.347 -5347.375new theta= 8.025556theta loglik c.loglik[1,] 0.000000 -5846.216 -5846.216[2,] 1.000000 -5305.049 -5590.102[3,] 2.000000 -5036.927 -5479.479[4,] 4.000000 -4740.394 -5385.887[5,] 8.000000 -4457.094 -5347.375[6,] 16.000000 -4223.785 -5393.362[7,] 8.740343 -4423.925 -5348.128[8,] 8.058000 -4454.347 -5347.375[9,] 8.025556 -4455.875 -5347.369new theta= 8.028123Warning messages:Inner loop failed to coverge for iterations 5 6 in: coxpenal.fit(X, Y,strats, offset, init = init, iter.max = iter.max, outer.max ....> fitc3Call:coxph(formula = Surv(time, status) ~ rx + extent + node4 + frailty(id, trace =T) + strata(etype), data = colon)coef se(coef) se2 Chisq DF prxLev 0.0434 0.305 0.140 0.02 1 8.9e-01rxLev+5FU -0.5125 0.310 0.170 2.73 1 9.8e-02extent 1.3373 0.251 0.137 28.45 1 9.6e-08node4 2.3381 0.233 0.156 100.81 1 0.0e+00frailty(id, trace = T) 5939.97 867 0.0e+00Iterations: 9 outer, 112 Newton-RaphsonVariance of random effect= 8.03 EM likelihood = -5347.4Degrees of freedom for terms= 0.5 0.3 0.4 866.7Likelihood ratio test=3787 on 868.01 df, p=0 n= 1858> > fitc4 <- coxph(Surv(time, status) ~ rx + extent + node4 + frailty(id, df=30)+ + strata(etype), colon)> fitc4Call:coxph(formula = Surv(time, status) ~ rx + extent + node4 + frailty(id, df = 30) +strata(etype), data = colon)coef se(coef) se2 Chisq DF prxLev -0.0374 0.0789 0.0769 0.22 1 6.4e-01rxLev+5FU -0.4565 0.0859 0.0840 28.27 1 1.1e-07extent 0.5289 0.0815 0.0798 42.13 1 8.5e-11node4 0.9078 0.0701 0.0681 167.85 1 0.0e+00frailty(id, df = 30) 58.56 30 1.4e-03Iterations: 3 outer, 9 Newton-RaphsonVariance of random effect= 0.0337 EM likelihood = -5832.4Degrees of freedom for terms= 1.9 1.0 0.9 30.0Likelihood ratio test=363 on 33.81 df, p=0 n= 1858> > # Do a fit, removing the no-event people> temp <- tapply(colon$status, colon$id, sum)> keep <- !(is.na(match(colon$id, names(temp[temp>0]))))> fitc5 <- coxph(Surv(time, status) ~ rx + extent + node4 +cluster(id)+ + strata(etype), colon, subset=keep)> > #> # Do the factor fit, but first remove the no-event people> #> # Ha! This routine has a factor with 506 levels. It uses all available> # memory, and can't finish in my patience window. Commented out.>> #fitc4 <- coxph(Surv(time, status) ~ rx + extent + node4 + factor(id), colon,> # subset=keep)>>>>>>> #> # The residual methods treat a sparse frailty as a fixed offset with> # no variance> #>> kfit1 <- coxph(Surv(time, status) ~ age + sex ++ frailty(id, dist="gauss"), kidney)> tempf <- predict(kfit1, type="terms")[,3]> temp <- kfit1$frail[match(kidney$id, sort(unique(kidney$id)))]> #all.equal(unclass(tempf), unclass(temp))> all.equal(as.vector(tempf), as.vector(temp))[1] T> > # Now fit a model with explicit offset> kfitx <- coxph(Surv(time, status) ~ age + sex + offset(tempf),kidney,+ eps=1e-7)> > # These are not precisely the same, due to different iteration paths> all.equal(kfitx$coef, kfit1$coef)[1] T> > # This will make them identical> kfitx <- coxph(Surv(time, status) ~ age + sex + offset(temp),kidney,+ iter=0, init=kfit1$coef)> all.equal(resid(kfit1), resid(kfitx))[1] T> all.equal(resid(kfit1, type="score"), resid(kfitx, type="score"))[1] T> all.equal(resid(kfit1, type="schoe"), resid(kfitx, type="schoe"))[1] T> > # These are not the same, due to a different variance matrix> # The frailty model's variance is about 2x the naive "assume an offset" var> # The score residuals are equal, however.> all.equal(resid(kfit1, type="dfbeta"), resid(kfitx, type="dfbeta"))[1] "Mean relative difference: 0.5214642"> zed <- kfitx> zed$var <- kfit1$var> all.equal(resid(kfit1, type="dfbeta"), resid(zed, type="dfbeta"))[1] T> > > temp1 <- resid(kfit1, type="score")> temp2 <- resid(kfitx, type="score")> all.equal(temp1, temp2)[1] T> > #> # Now for some tests of predicted values> #> all.equal(predict(kfit1, type="expected"), predict(kfitx, type="expected"))[1] T> all.equal(predict(kfit1, type="lp"), predict(kfitx, type="lp"))[1] T> > temp1 <- predict(kfit1, type="terms", se.fit=T)> temp2 <- predict(kfitx, type="terms", se.fit=T)> all.equal(temp1$fit[,1:2], temp2$fit)[1] T> all.equal(temp1$se.fit[,1:2], temp2$se.fit) #should be false[1] "Mean relative difference: 0.3023202"> mean(temp1$se.fit[,1:2]/ temp2$se.fit)[1] 1.432742> all.equal(as.vector(temp1$se.fit[,3])^2,+ as.vector(kfit1$fvar[match(kidney$id, sort(unique(kidney$id)))]))[1] T> > print(temp1)$fit:age sex frailty(id, dist = "gauss")1 -0.073958502 1.0394106 0.597861112 -0.073958502 1.0394106 0.597861113 0.020271945 -0.3712181 0.384858324 0.020271945 -0.3712181 0.384858325 -0.055112412 1.0394106 0.202075836 -0.055112412 1.0394106 0.202075837 -0.059823935 -0.3712181 -0.559114858 -0.055112412 -0.3712181 -0.559114859 -0.158765904 1.0394106 0.2854987310 -0.158765904 1.0394106 0.2854987311 -0.130496770 -0.3712181 0.0662606112 -0.125785247 -0.3712181 0.0662606113 0.034406512 1.0394106 0.8045900014 0.034406512 1.0394106 0.8045900015 0.053252601 -0.3712181 -0.4381282316 0.057964123 -0.3712181 -0.4381282317 0.119213914 -0.3712181 -0.0562658218 0.119213914 -0.3712181 -0.0562658219 0.034406512 1.0394106 -0.4995268320 0.039118034 1.0394106 -0.4995268321 0.001425855 -0.3712181 -0.1302046122 0.001425855 -0.3712181 -0.1302046123 -0.045689368 -0.3712181 0.0637408124 -0.045689368 -0.3712181 0.0637408125 -0.040977845 -0.3712181 0.3879628926 -0.040977845 -0.3712181 0.3879628927 -0.007997189 -0.3712181 -0.4762419028 -0.007997189 -0.3712181 -0.4762419029 -0.125785247 -0.3712181 -0.6695487930 -0.125785247 -0.3712181 -0.6695487931 0.076810213 1.0394106 0.1935241432 0.076810213 1.0394106 0.1935241433 0.076810213 -0.3712181 -0.1647446934 0.076810213 -0.3712181 -0.1647446935 -0.003285667 -0.3712181 -0.1578784136 0.001425855 -0.3712181 -0.1578784137 0.043829556 -0.3712181 -0.4620928338 0.043829556 -0.3712181 -0.4620928339 0.001425855 -0.3712181 0.1259611540 0.001425855 -0.3712181 0.1259611541 0.010848900 1.0394106 -1.7424181642 0.015560422 1.0394106 -1.7424181643 -0.064535457 -0.3712181 -0.4519117944 -0.064535457 -0.3712181 -0.4519117945 0.086233257 -0.3712181 0.5154889646 0.090944780 -0.3712181 0.5154889647 -0.007997189 -0.3712181 0.0946934848 -0.003285667 -0.3712181 0.0946934849 -0.003285667 1.0394106 0.0579554850 -0.003285667 1.0394106 0.0579554851 0.062675646 -0.3712181 -0.3791546352 0.067387168 -0.3712181 -0.3791546353 -0.158765904 -0.3712181 0.1124313054 -0.158765904 -0.3712181 0.1124313055 0.039118034 -0.3712181 0.5476257456 0.039118034 -0.3712181 0.5476257457 0.043829556 1.0394106 0.4585691458 0.043829556 1.0394106 0.4585691459 0.048541079 -0.3712181 0.3562396760 0.048541079 -0.3712181 0.3562396761 0.057964123 -0.3712181 0.4877920262 0.057964123 -0.3712181 0.4877920263 0.029694989 -0.3712181 0.2558178364 0.034406512 -0.3712181 0.2558178365 0.062675646 -0.3712181 0.2304640166 0.062675646 -0.3712181 0.2304640167 0.001425855 -0.3712181 -0.1367210868 0.006137378 -0.3712181 -0.1367210869 -0.102227636 -0.3712181 0.5195093070 -0.102227636 -0.3712181 0.5195093071 -0.007997189 -0.3712181 -0.2386267472 -0.007997189 -0.3712181 -0.2386267473 0.039118034 -0.3712181 0.1716482474 0.039118034 -0.3712181 0.1716482475 0.076810213 1.0394106 -0.3579894176 0.076810213 1.0394106 -0.35798941$se.fit:age sex frailty(id, dist = "gauss")1 0.195822035 0.3279661 0.62449192 0.195822035 0.3279661 0.62449193 0.053674606 0.1171308 0.69525954 0.053674606 0.1171308 0.69525955 0.145922707 0.3279661 0.57040616 0.145922707 0.3279661 0.57040617 0.158397539 0.1171308 0.48935548 0.145922707 0.1171308 0.48935549 0.420369012 0.3279661 0.606982210 0.420369012 0.3279661 0.606982211 0.345520020 0.1171308 0.563265912 0.333045188 0.1171308 0.563265913 0.091099103 0.3279661 0.663992314 0.091099103 0.3279661 0.663992315 0.140998431 0.1171308 0.510081516 0.153473263 0.1171308 0.510081517 0.315646080 0.1171308 0.549030718 0.315646080 0.1171308 0.549030719 0.091099103 0.3279661 0.526281320 0.103573935 0.3279661 0.526281321 0.003775278 0.1171308 0.517999222 0.003775278 0.1171308 0.517999223 0.120973042 0.1171308 0.620710624 0.120973042 0.1171308 0.620710625 0.108498210 0.1171308 0.581013426 0.108498210 0.1171308 0.581013427 0.021174386 0.1171308 0.624577628 0.021174386 0.1171308 0.624577629 0.333045188 0.1171308 0.561445330 0.333045188 0.1171308 0.561445331 0.203372591 0.3279661 0.653050832 0.203372591 0.3279661 0.653050833 0.203372591 0.1171308 0.524611734 0.203372591 0.1171308 0.524611735 0.008699554 0.1171308 0.510563336 0.003775278 0.1171308 0.510563337 0.116048767 0.1171308 0.628233838 0.116048767 0.1171308 0.628233839 0.003775278 0.1171308 0.631823040 0.003775278 0.1171308 0.631823041 0.028724942 0.3279661 0.523380542 0.041199774 0.3279661 0.523380543 0.170872371 0.1171308 0.549077344 0.170872371 0.1171308 0.549077345 0.228322255 0.1171308 0.605727746 0.240797087 0.1171308 0.605727747 0.021174386 0.1171308 0.626622448 0.008699554 0.1171308 0.626622449 0.008699554 0.3279661 0.552545450 0.008699554 0.3279661 0.552545451 0.165948095 0.1171308 0.555532852 0.178422927 0.1171308 0.555532853 0.420369012 0.1171308 0.584826054 0.420369012 0.1171308 0.584826055 0.103573935 0.1171308 0.608035756 0.103573935 0.1171308 0.608035757 0.116048767 0.3279661 0.600892358 0.116048767 0.3279661 0.600892359 0.128523599 0.1171308 0.576087960 0.128523599 0.1171308 0.576087961 0.153473263 0.1171308 0.598113862 0.153473263 0.1171308 0.598113863 0.078624270 0.1171308 0.661206564 0.091099103 0.1171308 0.661206565 0.165948095 0.1171308 0.560832566 0.165948095 0.1171308 0.560832567 0.003775278 0.1171308 0.584346868 0.016250110 0.1171308 0.584346869 0.270671027 0.1171308 0.608816870 0.270671027 0.1171308 0.608816871 0.021174386 0.1171308 0.679336572 0.021174386 0.1171308 0.679336573 0.103573935 0.1171308 0.641995274 0.103573935 0.1171308 0.641995275 0.203372591 0.3279661 0.577811576 0.203372591 0.3279661 0.5778115> kfit1Call:coxph(formula = Surv(time, status) ~ age + sex + frailty(id, dist = "gauss"),data = kidney)coef se(coef) se2 Chisq DF page 0.00471 0.0125 0.00856 0.14 1.0 0.7100sex -1.41063 0.4451 0.31503 10.04 1.0 0.0015frailty(id, dist = "gauss 26.54 14.7 0.0290Iterations: 6 outer, 28 Newton-RaphsonVariance of random effect= 0.569Degrees of freedom for terms= 0.5 0.5 14.7Likelihood ratio test=47.5 on 15.69 df, p=4.65e-05 n= 76> kfitxCall:coxph(formula = Surv(time, status) ~ age + sex + offset(temp), data = kidney,init = kfit1$coef, iter.max = 0)coef exp(coef) se(coef) z page 0.00471 1.005 0.00875 0.538 5.9e-01sex -1.41063 0.244 0.30916 -4.563 5.0e-06Likelihood ratio test=0 on 2 df, p=1 n= 76> > rm(temp1, temp2, kfitx, zed, tempf)> #> # The special case of a single sparse frailty> #>> kfit1 <- coxph(Surv(time, status) ~ frailty(id, dist="gauss"), kidney)> tempf <- predict(kfit1, type="terms")> temp <- kfit1$frail[match(kidney$id, sort(unique(kidney$id)))]> all.equal(as.vector(tempf), as.vector(temp))[1] T> > # Now fit a model with explicit offset> kfitx <- coxph(Surv(time, status) ~ offset(tempf),kidney, eps=1e-7)> > all.equal(resid(kfit1), resid(kfitx))[1] T> all.equal(resid(kfit1, type="deviance"), resid(kfitx, type="deviance"))[1] T> > #> # Some tests of predicted values> #> aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))> aeq(predict(kfit1, type="expected"), predict(kfitx, type="expected"))[1] T> aeq(predict(kfit1, type="lp"), predict(kfitx, type="lp"))[1] T> > temp1 <- predict(kfit1, type="terms", se.fit=T)> all.equal(temp1$fit, kfitx$linear)[1] T> all.equal(temp1$se.fit^2,+ kfit1$fvar[match(kidney$id, sort(unique(kidney$id)))])[1] T> > temp1$fit:[1] 0.695322941 0.695322941 0.244363776 0.244363776 0.493777896[6] 0.493777896 -0.658879753 -0.658879753 0.520994314 0.520994314[11] -0.114379760 -0.114379760 0.799261991 0.799261991 -0.487795409[16] -0.487795409 -0.120271844 -0.120271844 0.131081160 0.131081160[21] -0.214822375 -0.214822375 -0.054773695 -0.054773695 0.184575874[26] 0.184575874 -0.509525783 -0.509525783 -0.790241403 -0.790241403[31] 0.324356959 0.324356959 -0.239177390 -0.239177390 -0.264230973[36] -0.264230973 -0.472242362 -0.472242362 0.006350496 0.006350496[41] -0.872904122 -0.872904122 -0.530513459 -0.530513459 0.351179181[46] 0.351179181 -0.037129299 -0.037129299 0.441718240 0.441718240[51] -0.418876639 -0.418876639 -0.107887816 -0.107887816 0.346099359[56] 0.346099359 0.658680102 0.658680102 0.197185714 0.197185714[61] 0.304679320 0.304679320 0.139630271 0.139630271 0.093562534[66] 0.093562534 -0.209483283 -0.209483283 0.301887374 0.301887374[71] -0.278636264 -0.278636264 0.068590872 0.068590872 0.078473259[76] 0.078473259$se.fit:[1] 0.6147113 0.6147113 0.6157024 0.6157024 0.5713312 0.5713312 0.4392271[8] 0.4392271 0.5759022 0.5759022 0.4832712 0.4832712 0.6417752 0.6417752[15] 0.4573402 0.4573402 0.4812024 0.4812024 0.5118065 0.5118065 0.4762681[22] 0.4762681 0.5530110 0.5530110 0.5193604 0.5193604 0.5531610 0.5531610[29] 0.4773814 0.4773814 0.6361045 0.6361045 0.4707460 0.4707460 0.4669424[36] 0.4669424 0.5597930 0.5597930 0.5639378 0.5639378 0.4648879 0.4648879[43] 0.4902904 0.4902904 0.5446419 0.5446419 0.5567662 0.5567662 0.5606017[50] 0.5606017 0.4994076 0.4994076 0.4830103 0.4830103 0.5450207 0.5450207[57] 0.6054682 0.6054682 0.5207580 0.5207580 0.5374623 0.5374623 0.5908553[64] 0.5908553 0.5063650 0.5063650 0.5288166 0.5288166 0.5366511 0.5366511[71] 0.5992889 0.5992889 0.5760201 0.5760201 0.5780104 0.5780104> kfit1Call:coxph(formula = Surv(time, status) ~ frailty(id, dist = "gauss"), data = kidney)coef se(coef) se2 Chisq DF pfrailty(id, dist = "gauss 23 13.8 0.057Iterations: 6 outer, 28 Newton-RaphsonVariance of random effect= 0.457Degrees of freedom for terms= 13.8Likelihood ratio test=33.4 on 13.83 df, p=0.00234 n= 76> > > # From Gail, Sautner and Brown, Biometrics 36, 255-66, 1980>> # 48 rats were injected with a carcinogen, and then randomized to either> # drug or placebo. The number of tumors ranges from 0 to 13; all rats were> # censored at 6 months after randomization.>> # Variables: rat, treatment (1=drug, 0=control), o> # observation # within rat,> # (start, stop] status> # The raw data has some intervals of zero length, i.e., start==stop.> # We add .1 to these times as an approximate solution> #> rat2 <- read.table("data.rat2", col.names=c("id", "rx", "enum", "start",+ "stop", "status"))> temp1 <- rat2$start> temp2 <- rat2$stop> for (i in 1:nrow(rat2)) {+ if (temp1[i] == temp2[i]) {+ temp2[i] <- temp2[i] + .1;+ if (i < nrow(rat2) && rat2$id[i] == rat2$id[i+1]) {+ temp1[i+1] <- temp1[i+1] + .1;+ if (temp2[i+1] <= temp1[i+1]) temp2[i+1] _ temp1[i+1];+ };+ };+ }> rat2$start <- temp1> rat2$stop <- temp2> > r2fit0 <- coxph(Surv(start, stop, status) ~ rx + cluster(id), rat2)> > r2fitg <- coxph(Surv(start, stop, status) ~ rx + frailty(id), rat2)> r2fitm <- coxph(Surv(start, stop, status) ~ rx + frailty.gaussian(id), rat2)> > r2fit0Call:coxph(formula = Surv(start, stop, status) ~ rx + cluster(id), data = rat2)coef exp(coef) se(coef) robust se z prx -0.826 0.438 0.151 0.204 -4.04 5.3e-05Likelihood ratio test=32.9 on 1 df, p=9.89e-09 n= 253> r2fitgCall:coxph(formula = Surv(start, stop, status) ~ rx + frailty(id), data = rat2)coef se(coef) se2 Chisq DF prx -0.838 0.219 0.152 14.6 1.0 0.00013frailty(id) 57.3 26.4 0.00045Iterations: 7 outer, 21 Newton-RaphsonVariance of random effect= 0.317 EM likelihood = -779.1Degrees of freedom for terms= 0.5 26.3Likelihood ratio test=120 on 26.83 df, p=8.43e-14 n= 253> r2fitmCall:coxph(formula = Surv(start, stop, status) ~ rx + frailty.gaussian(id), data =rat2)coef se(coef) se2 Chisq DF prx -0.79 0.22 0.154 12.9 1.0 3.3e-04frailty.gaussian(id) 61.0 24.9 7.3e-05Iterations: 5 outer, 17 Newton-RaphsonVariance of random effect= 0.303Degrees of freedom for terms= 0.5 24.9Likelihood ratio test=118 on 25.37 df, p=6.99e-14 n= 253> > #This example is unusual: the frailties variances end up about the same,> # but the effect on rx differs. Double check it> # Because of different iteration paths, the coef won't be exactly the> # same, but darn close.>> temp <- coxph(Surv(start, stop, status) ~ rx + offset(r2fitm$frail[id]), rat2)> all.equal(temp$coef, r2fitm$coef[1])[1] "Mean relative difference: 0.0005510633"> > temp <- coxph(Surv(start, stop, status) ~ rx + offset(r2fitg$frail[id]), rat2)> all.equal(temp$coef, r2fitg$coef[1])[1] "Mean relative difference: 0.0005360911"> > #> # What do I get with AIC> #> r2fita1 <- coxph(Surv(start, stop, status) ~ rx + frailty(id, method="aic"),+ rat2)> r2fita2 <- coxph(Surv(start, stop, status) ~ rx + frailty(id, method="aic",+ dist="gauss"), rat2)> r2fita3 <- coxph(Surv(start, stop, status) ~ rx + frailty(id, dist="t"),+ rat2)> > r2fita1Call:coxph(formula = Surv(start, stop, status) ~ rx + frailty(id, method = "aic"),data = rat2)coef se(coef) se2 Chisq DF prx -0.838 0.23 0.151 13.3 1.0 0.00026frailty(id, method = "aic 60.4 28.2 0.00039Iterations: 10 outer, 25 Newton-RaphsonVariance of random effect= 0.375 EM likelihood = -779.2Degrees of freedom for terms= 0.4 28.2Likelihood ratio test=124 on 28.62 df, p=7.92e-14 n= 253> r2fita2Call:coxph(formula = Surv(start, stop, status) ~ rx + frailty(id, method = "aic",dist = "gauss"), data = rat2)coef se(coef) se2 Chisq DF prx -0.784 0.255 0.154 9.48 1.0 2.1e-03frailty(id, method = "aic 73.40 29.7 1.4e-05Iterations: 6 outer, 18 Newton-RaphsonVariance of random effect= 0.493Degrees of freedom for terms= 0.4 29.7Likelihood ratio test=127 on 30.03 df, p=6.02e-14 n= 253> r2fita3Call:coxph(formula = Surv(start, stop, status) ~ rx + frailty(id, dist = "t"), data= rat2)coef se(coef) se2 Chisq DF prx -0.789 0.254 0.157 9.67 1 0.00190frailty(id, dist = "t") 64.70 30 0.00024Iterations: 7 outer, 23 Newton-RaphsonVariance of random effect= 0.779Degrees of freedom for terms= 0.4 30.0Likelihood ratio test=126 on 30.37 df, p=1.39e-13 n= 253> q()