Ramon Diaz-Uriarte asks a hard question about nested models containing a frailty term. Assume that factorA is binary. > fit1<-coxph(Surv(time,censor)~other.stuff+factorA+ frailty(id)) > fit2<-coxph(Surv(time,censor)~other.stuff+ frailty(id)) > chi1 <- 2*diff(fit1$loglik) > chi2 <- 2*diff(fit2$loglik) > 1-pchisq(chi1-chi2, df=sum(fit1$df)-sum(fit2$df)) The above can give a negative chisquare and/or degrees of freedom, which is certainly disturbing. There are two ways to look at the frailty term: as a single extra variance parameter (integrated likelihood view), or as a large collection of coefficients subject to a shrinkage constraint (penalized likelihood view). 1. If we take the first view, then the two models should be compared based on the integrated likelihood, the term currently labeled "EM likelihood" on the printout. This has integrated out the individual frailty coefficients. The test comparing the two models will have 1 degree of freedom. In the example below, this will give 2*(181.7-179.1)= 5.1 on 1 df: disease is a significant covariate. a. Pro: the theory is worked out, & seems reliable b. Con: theory is only present for a single term, and that term must be a factor variable. (At least for the distributions supported in Splus). 2. The second view treats id as a factor variable, with a penalty term to keep the number of degrees of freedom `sensible'. This is the way that smoothing spline (pspline) inference is done. But if this approach is taken, the above code may not be valid because the frailty term has not necessarily retained the same number of degrees of freedom for the two runs. For example: > coxph(Surv(time, status) ~ sex + disease + frailty(id), data=kidney) coef se(coef) se2 Chisq DF p sex -1.477 0.357 0.357 17.14 1 3.5e-05 diseaseGN 0.139 0.364 0.364 0.15 1 7.0e-01 diseaseAN 0.413 0.336 0.336 1.51 1 2.2e-01 diseasePKD -1.367 0.589 0.589 5.39 1 2.0e-02 frailty(id) 0.00 0 9.3e-01 Iterations: 6 outer, 28 Newton-Raphson Variance of random effect= 5e-07 EM likelihood = -179.1 Degrees of freedom for terms= 1 3 0 Likelihood ratio test=17.6 on 4 df, p=0.0015 n= 76 > coxph(Surv(time, status) ~ sex + frailty(id), data=kidney) coef se(coef) se2 Chisq DF p sex -1.56 0.454 0.351 11.8 1.0 0.00058 frailty(id) 23.1 13.2 0.04400 Iterations: 6 outer, 36 Newton-Raphson Variance of random effect= 0.399 EM likelihood = -181.7 Degrees of freedom for terms= 0.6 13.2 Likelihood ratio test=45.8 on 13.78 df, p=2.63e-05 n= 76 The second model has dropped one term, but the degrees of freedom has actually gone UP, due to a different frailty solution. To get an "honest" test of the disease term, compare the second model instead to > coxph(Surv(time, status) ~ sex + disease + frailty(id, theta=.399), data=kidney) coef se(coef) se2 Chisq DF p sex -1.731 0.465 0.369 13.86 1.0 0.0002 diseaseGN 0.268 0.496 0.369 0.29 1.0 0.5900 diseaseAN 0.493 0.451 0.343 1.20 1.0 0.2700 diseasePKD -0.788 0.812 0.620 0.94 1.0 0.3300 frailty(id, theta = 0.399 15.86 11.8 0.1900 Iterations: 1 outer, 7 Newton-Raphson Variance of random effect= 0.399 EM likelihood = -180.1 Degrees of freedom for terms= 0.6 1.7 11.8 Likelihood ratio test=46.7 on 14.22 df, p=2.5e-05 n= 76 Now the test for addition of disease is (46.7 - 45.8)=0.9 on (14.22 - 13.78) = 0.44 df. (The kidney data set is an extreme case, since both the frailty and disease effects are dominated by a single large outlier. Rather innocent looking changes in the model can lead to very different solutions.) a. Pro: The approach works for multiple frailty terms, and for all frailty distributions. Related to the Wald test(s) beta/se(beta) Integrates nicely with AIC and BIC selection criteria b. Con: The asymptotic framework is shaky. The number of "parameters" increases with n. --- The jury is still out on which outlook is the better one. In fact, I'm not all that sure that the jury has begun deliberations. But I am very interested in the answer.