[R] issues with calling predict.coxph.penal (survival) inside a function

[julian.bothe at elitepartner.de]

Hello,
and thanks for the answer.

1) I found a work-around - in  the end it is easier than thought before.
The only thing you have to do is to have the same variable name with the
new values.
So if predict(coxph.penal.fit, newdata[subset,]) does not work inside a
function,  the following works:

pred_function <- function(coxph_model, newdata){
#things to do before
newdata=newdata[the_subset_i_want,]
predict(coxph_model, newdata)

}

I attach a working example.

I am not sure, maybe this is even what is written in the help to
NextMethod ;)
" NextMethod works by creating a special call frame for the next method.
If no new arguments are supplied, the arguments will be the same in
number, order and name as those to the current method but their values
will be promises to evaluate their name in the current method and
environment. Any named arguments matched to ... are handled specially:
they either replace existing arguments of the same name or are appended to
the argument list. They are passed on as the promise that was supplied as
an argument to the current environment. (S does this differently!) *If
they have been evaluated in the current (or a previous environment) they
remain evaluated.* (This is a complex area, and subject to change: see the
draft ‘R Language Definition’.)" (Help to NextMethod )

2) Terry, I am not sure about the work-around you provided in your mail. I
want to do subsetting on newdata, not on the model.
Additionally, when trying the example you provided, I received different
results. Example is attached.

Thanks and all the best
Julian

#---------------------------

test1 <- data.frame(time=c(4,3,1,1,2,2,3),
              status=c(1,1,1,0,1,1,0),
              x=c(0,2,1,1,1,0,0),
              sex=c(0,0,0,0,1,1,1))

# Fit a stratified model
fit1 <- coxph(Surv(time, status) ~ x + strata(sex), test1)
summary(fit1)

#fit stratified wih spline
fit2 <- coxph(Surv(time, status) ~ pspline(x, df=2) + strata(sex), test1)
summary(fit2)

## work-around
predicting_function_which_works<- function(model, newdata ){
  subs <-vector(mode='logical', length=nrow(newdata))
  subs[3:length(subs)]<- TRUE #try with first values set to false

  newdata_alt<-newdata
  newdata<-newdata[subs,]

  ret<-vector(mode='numeric', length=nrow(newdata_alt))
  ret[!subs]<- NA
  ret[subs]<- predict(model,newdata )
  return(ret)
}

predicting_function_which_works(fit1, test1) # works

predicting_function_which_works(fit2,test1)  # works
predicting_function_which_works(fit2,data.frame(time=c(4,3,1,1,2), # works
                                              status=c(1,1,1,0,1),
                                              x=c(0,2,1,1,2),
                                              sex=c(1,1,0,0,1))
)

## How I understood Terry's work-around. Provides different results and
doesn't consider subset

predicting_function_2 <- function(model, newdata){
  subs <-vector(mode='logical', length=nrow(newdata))
  subs[2:length(subs)]<- TRUE

  newX <- model.matrix(model)
  newY <- model$y
  newfit <- coxph(newY ~ newX, iter=0, init=coef(model))
  newfit$var <- model$var

  #print(model)
  #print(newfit)
  #predict(newfit)

  #ret=predict(newfit)

  print("comparison")
  print(paste(" model, original prediction:", paste(predict(model),
collapse=",")))
  print(paste("newfit, original prediction:", paste(predict(newfit),
collapse=",")))

  ret <- predict (newfit, newdata[subs,])
  return(ret)
}

predicting_function_2(fit1, test1)

predicting_function_2(fit2,test1)


-----Ursprüngliche Nachricht-----
Von: Terry Therneau [mailto:therneau at mayo.edu]
Gesendet: Donnerstag, 14. November 2013 16:31
An: r-help at r-project.org; julian.bothe at elitepartner.de
Betreff: Re: issues with calling predict.coxph.penal (survival) inside a
function

Thanks for the reproducable example.  I can confirm that it fails on my
machine using survival 2-37.5, the next soon-to-be-released version,

The issue is with NextMethod, and my assumption that the called routine
inherited everything from the parent, including the environment chain.  A
simple test this AM showed me that the assumption is false.  It might have
been true for Splus.  Working this out may take some time -- every other
one of my wrestling matches with predict inside a function has -- and
there is a reasonable chance that it won't make this already overdue
release.

In the meantime, here is a workaround that I have sometimes used in other
situations.
Inside your function do the following: fit a new coxph model with fixed
coefficients, and do prediction on that.

myfun <- function(oldfit, subset) {
    newX <- model.matrix(oldfit)[subset,]
    newY <- oldfit$y[subset]
    newfit <- coxph(newY ~ newX, iter=0, init=coef(oldfit))
    newfit$var <- oldfit$var
    predict(newfit)
    }

If the subset is all of a particular strata, as you indicated, then all of
the predictions will be correct.  If not, then those that make use of the
the baseline hazard (type=
expect) will be incorrect but all others are ok.

Terry Therneau


On 11/14/2013 05:00 AM, r-help-request at r-project.org wrote:
> Hello everyone,
>
>
>
> I got an issue with calling predict.coxph.penal inside a function.
>
>
>
> Regarding the context: My original problem is that I wrote a function
> that uses predict.coxph and survfit(model) to predict
>
> a lot of survival-curves using only the basis-curves for the strata
> (as delivered by survfit(model) ) and then adapts them with
>
> the predicted risk-scores. Because there are cases where my new data
> has strata which didn't exist in the original model I exclude
>
> them, using a Boolean vector inside the function.
>
> I end up with a call like this: predict (coxph_model,
> newdata[subscript_vector,] )
>
>
>
> This works fine for coxph.model, but when I fit a model with a spline
> (class coxph.penal), I get an error:
>
> "Error in `[.data.frame`(newdata, [subscript_vector, ) : object
> '[subscript_vector ' not found"
>
>
>
> I suppose this is because of NextMethod, but I am not sure how to work
> around it. I also read a little bit about all those
> matching-and-frame-issues,
>
> But must confess I am not really into it.
>