I'm trying to add time-varying coefficients using the tt() function with coxph to solve the problem of PH violations in my Cox model.

When I try to do this I encounter the error:

Cannot allocate vector of size...

I guess this is because of what tt() does to the matrix (and my dataset is n=100 000 with 135 variables), however, the problem is I'm using a remote server and I am limited to their capacities.

So, I guess my question is whether there may be other packages that are more memory efficient that would allow me to do this.

Each ID is represented by one record in the data and there are 1627 events.

  • $\begingroup$ How many unique time values are there? How many events among the 100,000 cases? Please edit the question to add that information, as comments are easy to overlook and can be deleted. $\endgroup$
    – EdM
    Aug 10, 2022 at 16:43
  • $\begingroup$ @EdM Thanks for the tip, updated $\endgroup$
    – JED HK
    Aug 11, 2022 at 8:20

1 Answer 1


I would recommend not trying to use tt() to fit continuous functions of time for your coefficients in this case. With smaller data sets that helps describe and adjust for the pattern of associations of predictors with outcome over time. But even after you've done that I don't think that you can use cox.zph() to show that you "fixed" the proportional hazards (PH) problem. And with your large data set it seems to be impossible in practice, anyway, as you need to make up to 1627 copies (one per event time) of hundreds of thousands of data rows.

I'd recommend a few steps instead.

First, make sure that you are modeling continuous predictors flexibly, e.g. with regression splines, rather than as simple linear predictors. A poorly modeled continuous predictor can lead to apparent violation of PH, perhaps even for other predictors in the model.

Second, follow Frank Harrell's advice in a comment on a related question and do some data reduction to cut down on the number of predictors. With 135 predictors you would expect about 7 of them to show "significant" violations of PH at p < 0.05; I understand that you only found 15 violators despite your large number of cases.

As you aren't interested in most of the predictors except for the effects of levels of a treatment, why use such a complicated model if you can adjust for the other predictors more simply? You are about at the limit of reliability for overfitting, with only about 12 events per predictor, so reducing the number of predictors could improve the precision of the parts of the model that you care most about. Section 4.7 of Harrell's course notes or book describes data-reduction methods that can simplify the model without biasing results; all of Chapter 4 deserves close study.

Third, if you still have a PH problem after that, try stratification instead of the tt() function. Stratifying on levels of categorical predictors that fail PH should fix problems associated with them. You don't get regression coefficients for predictors on which you stratify, but if you're just adjusting for them you don't care about their coefficients.

To get hazard ratios you don't want to stratify your categorical treatment levels of primary interest, but you can still use strata by time groups that interact with your multi-level treatment predictor. That's described in Section 4.1 of the R survival time-dependence vignette. With time strata you can use cox.zph() to evaluate PH.

If you use that approach to specify, say, 10 different time intervals you would only have to make 10 copies (at most, versus up to 1627 with tt()) of each of your 100,000 cases. That would probably fit into your computer, fix the PH problem if you choose appropriate time intervals, and would show how hazard ratios for treatment levels change with time, at whatever resolution you choose for the time intervals.

  • $\begingroup$ What a great answer Ed. I'm grateful as always. I maybe should have been clearer - although our main interest is the treatment, we still find it interesting to see coefficients of other covariates. This will limit the data reduction approaches available to us and to which variables we stratify. Oh and btw, with regard to modelling continuous predictors correctly, I of course used Martingale plots to check functional form; how can I check if the use of splines would improve my situation? $\endgroup$
    – JED HK
    Aug 12, 2022 at 9:26
  • 1
    $\begingroup$ @JEDHK include splines in the model. Restricted cubic splines via rcs() in Harrell's rms package provide a linear term and additional nonlinear terms modifying it. The survival package provides penalized smoothing splines via its pspline() function. Both allow for tests on nonlinearity. It's best to decide ahead of time how many degrees of freedom/knots to allocate to each continuous predictor and then stick with that for inference. See Chapters 2 and 4 of Harrell's notes or book. $\endgroup$
    – EdM
    Aug 12, 2022 at 15:12
  • $\begingroup$ Thanks for the suggestion Ed! I used 'pspline()' to add cubic splines of the continuous variables that violated the PH assumption and now the assumption is satisfied. That's great (in theory), but I checked the Schoenfeld plot for one of the covs (as I believe I've seen Harrell suggest doing) and the shape of the curve looks basically the same, though now the Y-axis is zoomed out (i.e., what was 0.1 is now around 1). This surprises me - I expected a straight horizontal line after fitting the spline. Do you have any thoughts on this? $\endgroup$
    – JED HK
    Aug 12, 2022 at 17:04
  • 1
    $\begingroup$ @JEDHK those plots are hard to interpret for spline terms, as they combine all of the multiple associated coefficients into 1 plot. I also see that behavior with pspline() or rcs() splines if you accept the default terms=TRUE in the call to cox.zph(). With an rcs() spline you can examine each spline coefficient individually if you specify terms=FALSE, but the penalization with pspline() doesn't seem to allow for that. Do check how many degrees of freedom you used up with your pspline()terms; if you had a lot of them you might be in danger of overfitting. $\endgroup$
    – EdM
    Aug 12, 2022 at 18:02
  • $\begingroup$ I went for 5 as I believe Harrell suggests (RMS, page 28) that for large datasets (>100) we can afford up to 5 knots. Since I have 100k I thought 5 seemed reasonable. I was also wondering about the individual associated coefficients - what's my interpretation for them then? That is, how can I calculate the contribution to the final model of the continuous variable that I've applied pspline to? $\endgroup$
    – JED HK
    Aug 13, 2022 at 7:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.