The R Project SVN
R-packages
(root)/trunk/survival/Examples/testfrail/Rtest.out.save @ 1663 – Rev
Rev 123 | Blame | Last modification | View Log | Download | RSS feed
R : Copyright 2004, The R Foundation for Statistical ComputingVersion 2.0.0 Under development (unstable) (2004-08-12), ISBN 3-900051-00-3R 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.R is a collaborative project with many contributors.Type 'contributors()' for more information and'citation()' on how to cite R in publications.Type 'demo()' for some demos, 'help()' for on-line help, or'help.start()' for a HTML browser interface to help.Type 'q()' to quit R.> #> # Set up for the test> #> #dyn.load("../loadmod.o")> #attach("../.Data")> #options(na.action="na.omit", contrasts='contr.treatment')> library(survival)Loading required package: splines> #library(date)> #> # 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] TRUE> all.equal(tfit1[temp], tfit3[temp])[1] TRUE>> 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] TRUE>>>> # 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 I-likelihood = -181.1Degrees of freedom for terms= 1.0 13.9Likelihood ratio test=36.3 on 14.8 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[1] TRUE$historytheta 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 I-likelihood = -181.5Degrees of freedom for terms= 1.0 22.7Likelihood ratio test=50.7 on 23.7 df, p=0.00100 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.0 14.6 0.24Iterations: 6 outer, 17 Newton-RaphsonVariance of random effect= 0.504 I-likelihood = -184.8Degrees of freedom for terms= 14.6Likelihood ratio test=30 on 14.6 df, p=0.0101 n= 150> #> # Here is a test case with multiple smoothing terms> #> data(lung)> 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.69Degrees 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.78Degrees 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> data(ovarian)> 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] TRUE$"pspline(age, df = 2)"$historythetas 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.7485 0.24 1.00 0.6200pspline(age, df = 4), lin 0.139 0.044 0.0440 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>> if (gc()[2,1]>60000){+ set.seed(42)+ 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.372 5.859 7.300 10.520[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) ##FAILS in S-PLUS[1] "Mean relative difference: 4.724822e-08"> all.equal(temp2, kfit1$frail) ##FAILS in S-PLUS[1] "names for target but not for current"[2] "Attributes: < Modes: list, NULL >"[3] "Attributes: < Lengths: 2, 0 >"[4] "Attributes: < names for target but not for current >"[5] "Attributes: < Length mismatch: comparison on first 0 components >"[6] "target is array, current is numeric">> 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 I-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 = 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 I-likelihood = -182.5Degrees of freedom for terms= 0.2 0.5 1.1 18.8Likelihood ratio test=63.8 on 20.6 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.48314 0.227 0.3582 -4.140 3.5e-05diseaseGN 0.08796 1.092 0.4064 0.216 8.3e-01diseaseAN 0.35079 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 I-likelihood = -182.5Degrees of freedom for terms= 0.4 0.5 1.4 18.7Likelihood ratio test=63.8 on 21.0 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.0150 0.0106 0.11 1.0 0.74000sex -1.69703 0.4609 0.3617 13.56 1.0 0.00023diseaseGN 0.17980 0.5447 0.3927 0.11 1.0 0.74000diseaseAN 0.39283 0.5447 0.3982 0.52 1.0 0.47000diseasePKD -1.13630 0.8250 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.9 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.601gauss:1 1.701 0.588 0.5181 5.586gauss:2 1.424 0.702 0.3851 5.266gauss:3 1.159 0.863 0.3828 3.511gauss:4 0.623 1.606 0.2340 1.657gauss:5 1.254 0.797 0.3981 3.953gauss:6 1.135 0.881 0.3834 3.360gauss:7 1.973 0.507 0.5694 6.834gauss:8 0.620 1.614 0.2166 1.772gauss:9 0.823 1.215 0.2888 2.346gauss:10 0.503 1.988 0.1747 1.448gauss:11 0.757 1.322 0.2708 2.113gauss:12 1.105 0.905 0.3343 3.651gauss:13 1.302 0.768 0.4275 3.967gauss:14 0.591 1.691 0.1854 1.885gauss:15 0.545 1.835 0.1858 1.598gauss:16 1.044 0.958 0.3142 3.470gauss:17 0.914 1.095 0.3000 2.782gauss:18 0.918 1.089 0.3248 2.597gauss:19 0.643 1.556 0.1951 2.117gauss:20 1.170 0.855 0.3453 3.963gauss:21 0.334 2.997 0.1020 1.091gauss:22 0.687 1.455 0.2353 2.006gauss:23 1.478 0.677 0.4756 4.592gauss:24 1.017 0.983 0.3156 3.278gauss:25 0.810 1.235 0.2749 2.384gauss:26 0.614 1.627 0.2149 1.757gauss:27 1.088 0.919 0.3282 3.610gauss:28 1.542 0.649 0.4923 4.829gauss:29 1.379 0.725 0.4377 4.342gauss:30 1.375 0.727 0.4444 4.253gauss:31 1.445 0.692 0.4703 4.438gauss:32 1.199 0.834 0.3521 4.085gauss:33 1.945 0.514 0.5523 6.849gauss:34 0.862 1.161 0.2769 2.682gauss:35 1.703 0.587 0.5266 5.508gauss:36 0.827 1.209 0.2281 3.002gauss:37 1.471 0.680 0.3894 5.555gauss: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.1 df, p=0Wald test = 37.4 on 15.1 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 I-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.0120 0.00873 0.22 1.0 0.6400sex -1.65036 0.4810 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.9 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 I-likelihood = -182.6Degrees of freedom for terms= 0.3 0.5 19.1Likelihood ratio test=60 on 19.9 df, p=6.83e-06 n= 76>> 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.0120 0.00873 0.22 1.0 0.6400sex -1.65036 0.4810 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.9 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(colon)> #data.restore('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.106 -0.343 7.3e-01rxLev+5FU -0.4488 0.638 0.0840 0.117 -3.842 1.2e-04extent 0.5155 1.674 0.0796 0.110 4.701 2.6e-06node4 0.8799 2.411 0.0681 0.096 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.046071> 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 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.028123> 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 I-likelihood = -5347.4Degrees of freedom for terms= 0.5 0.3 0.4 866.7Likelihood ratio test=3787 on 868 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 I-likelihood = -5832.4Degrees of freedom for terms= 1.9 1.0 0.9 30.0Likelihood ratio test=363 on 33.8 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))[1] TRUE> all.equal(as.vector(tempf), as.vector(temp))[1] TRUE>> # 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] TRUE>> # 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] TRUE> all.equal(resid(kfit1, type='score'), resid(kfitx, type='score'))[1] TRUE> all.equal(resid(kfit1, type='schoe'), resid(kfitx, type='schoe'))[1] TRUE>> # 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] TRUE>>> temp1 <- resid(kfit1, type='score')> temp2 <- resid(kfitx, type='score')> all.equal(temp1, temp2)[1] TRUE>> #> # Now for some tests of predicted values> #> all.equal(predict(kfit1, type='expected'), predict(kfitx, type='expected'))[1] TRUE> all.equal(predict(kfit1, type='lp'), predict(kfitx, type='lp'))[1] TRUE>> 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] TRUE> 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] TRUE>> print(temp1)$fitage 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.fitage 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.7 df, p=4.65e-05 n= 76> kfitxCall:coxph(formula = Surv(time, status) ~ age + sex + offset(temp),data = kidney, init = kfit1$coef, iter = 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] TRUE>> # 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] TRUE> all.equal(resid(kfit1, type='deviance'), resid(kfitx, type='deviance'))[1] TRUE>> #> # 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] TRUE> aeq(predict(kfit1, type='lp'), predict(kfitx, type='lp'))[1] TRUE>> temp1 <- predict(kfit1, type='terms', se.fit=T)> all.equal(temp1$fit, kfitx$linear)[1] TRUE> all.equal(temp1$se.fit^2,+ kfit1$fvar[match(kidney$id, sort(unique(kidney$id)))])[1] TRUE>> 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.0 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.8 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.827 0.438 0.151 0.204 -4.05 5.2e-05Likelihood ratio test=32.9 on 1 df, p=9.9e-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 I-likelihood = -779.1Degrees of freedom for terms= 0.5 26.3Likelihood ratio test=120 on 26.8 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.220 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.4 df, p=7e-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]) ##not quite[1] TRUE>> temp <- coxph(Surv(start, stop, status) ~ rx + offset(r2fitg$frail[id]), rat2)> all.equal(temp$coef, r2fitg$coef[1]) ##not quite[1] TRUE>> #> # 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.230 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 I-likelihood = -779.2Degrees of freedom for terms= 0.4 28.2Likelihood ratio test=124 on 28.6 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 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.79 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.4 df, p=1.39e-13 n= 253> q()