jgellar / pcox

Penalized Cox regression models
1 stars 0 forks source link

print.coxph #36

Closed jgellar closed 8 years ago

jgellar commented 9 years ago

Here is the current output when a pcox models is printed, which goes through print.coxph:

d> fit2
Call:
coxph(formula = Surv(time, event) ~ s(x) + male, data = pcoxdata, 
    na.action = function (object, ...) 
    {
        n <- length(object)
        omit <- logical(nrow(object))
        vars <- seq_len(n)
        for (j in vars) {
            x <- object[[j]]
            if (!is.atomic(x)) 
                next
            x <- is.na(x)
            d <- dim(x)
            if (is.null(d) || length(d) != 2L) 
                omit <- omit | x
            else omit <- omit | apply(x, 1, all)
        }
        xx <- object[!omit, , drop = FALSE]
        if (any(omit > 0L)) {
            temp <- setNames(seq(omit)[omit], attr(object, "row.names")[omit])
            attr(temp, "class") <- "omit"
            attr(xx, "na.action") <- temp
        }
        xx
    }, x = TRUE, iter.max = 100, outer.max = 50)

       coef    se(coef) se2    Chisq DF p      
s(x).1 -1.5676 0.343    0.2856 20.92 1  4.8e-06
s(x).2  0.0271 0.831    0.5396  0.00 1  9.7e-01
s(x).3 -0.1831 0.237    0.1453  0.60 1  4.4e-01
s(x).4 -0.0479 0.487    0.2467  0.01 1  9.2e-01
s(x).5 -0.0470 0.201    0.0942  0.05 1  8.1e-01
s(x).6 -0.0840 0.441    0.2132  0.04 1  8.5e-01
s(x).7 -0.0148 0.139    0.0483  0.01 1  9.1e-01
s(x).8 -0.1814 1.156    0.6913  0.02 1  8.8e-01
s(x).9  1.1387 0.500    0.3218  5.18 1  2.3e-02
male    0.5543 0.134    0.1336 17.17 1  3.4e-05

Iterations: 10 outer, 30 Newton-Raphson
Degrees of freedom for terms= 4.1 1.0 
Likelihood ratio test=98.9  on 5.13 df, p=0  n= 500 

There are two things that aren't ideal:

  1. The entire na.action function (na.omit_pcox) gets printed out. Any ideas how to get rid of it? Why does this not happen when one of the standard na.action functions (e.g., na.omit) is used? One thought was because I have na.omit_pcox as internal, but I don't think that would be the only cause of this - and when I exported the function nothing printed differently.
  2. There are 9 rows printed for the smooth function, one for each spline coefficient. Therneau actually built in a pretty nice system of over-riding this for coxph.penalty terms: you can replace those 9 rows with basically whatever you want. This is done by specifying a printfun function as an attribute of the coxph.penalty term. For example, in his pspline function, this is the output:
d> fit1
Call:
coxph(formula = Surv(time, status) ~ ph.ecog + pspline(age, 3), 
    data = cancer)

                        coef   se(coef) se2     Chisq DF   p      
ph.ecog                 0.4480 0.11707  0.11678 14.64 1.00 0.00013
pspline(age, 3), linear 0.0113 0.00928  0.00928  1.47 1.00 0.22000
pspline(age, 3), nonlin                          2.08 2.08 0.37000

Iterations: 4 outer, 12 Newton-Raphson
     Theta= 0.861 
Degrees of freedom for terms= 1.0 3.1 
Likelihood ratio test=21.9  on 4.08 df, p=0.000227
  n=227 (1 observation deleted due to missingness)

So he replaces his 10 spline coefficients with 2 lines, which he calls linear and nonlinear. linear I believe would be the coefficient for a simple linear term for age, and non-linear has something to do with the deviation for this. He adds a test in for the Chisq, DF, and p columns for each row - though I don't really believe this test at all. I'm guessing it's similar to Wood's test in summary.gam.

Anyways, we can replace our 9 rows with anything, with the caveat that it has to fall into these 6 columns. Blanks are allowed. Our options are:

  1. Leave it as is. I don't think this is that bad when there are only 9 basis functions, but for some models (e.g., bivariate smooths), there could be a lot more than that
  2. Take these rows out of the table completely. At least this avoids the problem of too much useless output.
  3. Find some way to summarize the term in a couple lines. We might have to have a different summary depending on the type of term (scalar, functional, historical, etc.).

3 would be ideal if we had a good idea on how to do that. You got anything? We don't have any tests developed yet - I have wanted to think about this at some point down the line.

jgellar commented 9 years ago

Also, it would be nice if we could replace the Call: line at the top with the pcox call, but I'm not too picky on this point.

jgellar commented 9 years ago

So I killed two birds with one stone, and wrote a short print.pcox method that replaces object$call with object$pcox$call and then calls print.coxph. So now the na.action function is not an issue. I'd still like your input on issue number 2 above though (what to include in the table for these terms).