Hazard Rate Modelling in R

Question

I am writing my thesis about bankruptcy prediction Pr(Y=1) through a sentiment score I have calculated (x1) and a control variable called Z-Score (x2). However I am very unsure about how to integrate my data in the R formula I have found. This is my code:

library(survival)
## Add survival object. status = 1 is bankruptcy
WRDS$SurvObj <- with(WRDS, Surv(as.numeric(DEL == 1)))

## model
res.cox1 <- coxph(SurvObj ~ SCORE10K+TIME+Z, data =  WRDS)
summary(res.cox1)

The following is my dataset, x1,x2,x3 ,x4 & x5 are sub variables that I used to create Z, so dont take those into consideration:

    CIK TIME DEL SCORE10K SCOREMDA AbsDiff   x1    x2   x3        x4   x5         Z
1 10254    1   0     0.69     0.13    0.56 0.24 -0.36 0.15 20.578713 0.39 13.016228
2 10254    2   0     0.66     0.13    0.53 0.25 -0.16 0.20 10.676150 0.45  7.591690
3 10254    3   0     0.65     0.18    0.47 0.02 -0.07 0.12  5.063745 0.33  3.690247
4 10254    4   0     0.62     0.19    0.43 0.06  0.03 0.20  6.476520 0.38  5.039912
5 20629    1   0     0.70     0.26    0.44 0.57  1.15 0.12  4.859852 1.34  6.945911
6 20629    2   0     0.74     0.30    0.44 0.61  1.17 0.13  6.950391 1.26  8.229235

Is it correct to plug in my continuous SCORE10K and Z variables together with the TIME variable? Or is that already integrated in the function through the Survival Object? (DEL is my binary variable that shows if there is bankruptcy Y=1).

TIME is coded as 1,2,3,4 for each company, and each company has a score for a duration of 4 years before they go bankrupt, OR NOT. My sample includes 50 companies that go default in year 4 and a matching healthy sample of 100 that also has scores for the same 4 years. So in total I have 600 data points but only 150 "company-specific data points".

I also gave random effects logit model a thought but in this case we are really confronted with a survival analysis if I am correct.

EDIT after input:

I have discovered that due to the structure of my dataset, the Cox model cox.zph function estimates my variance to be equal to 0, since all my default events Y=1 happen in the last period t=4. This invalidates the model technically, even though theoretically it makes sense to use a hazard model.

If you take a look at my dataset, that just means that DEL=1 happens only at TIME=4, IF it happens.

I cannot add data points anymore due to time constraints, so I am wondering whether I need to change the structure of the data points or if there is any way around the Cox model not being able to estimate my survival rates.

Answer 1

First, as a more convenient way:

model.coxph <- coxph(Surv(TIME, DEL) ~ SCORE10K + Z, ties="exact", data=WRDS)

Note several important things:

1. Your TIME variable cannot be used within the cox model. It is the underlying property of it. If done correctly, I think R will throw an error if you do.
2. Use the ties="exact" option to properly handle the fact that you use a small number of discrete time intervals. Generally, coxph deals with continuous time. See the coxph documentation for more information on this.
3. It is unclear to me if you have time-dependent covariates in your data. If you do, you need to add start_time and end_time to the Surv function.

Second, just to make sure that the data is build properly, you can see here for examples of time-dependent datasets and analysis, or here otherwise. The later link has excellent information on general cox hazard analysis and importantly - checking the proportionality assumptions.

Edit: after seeing your dataset -

So you DO have time-dependent covariates. Make sure to see Therneau's text (first link above). At any rate, you need to create a second TIME variable to encompass the time-frame in each row. If each unique ID (company) has 4 rows (or less if event occurred) than just add a 1-lagged of 1-added time variable like: WRDS$TIME_2 <- WRDS$TIME + 1. Then enter both starting and ending times to the Surv function:

model.coxph <- coxph(Surv(TIME, TIME_2, DEL) ~ SCORE10K + Z, ties="exact", data=WRDS)

Lastly, check for proportional assumption violations using cox.zph. You can see how to deal with testing and handling them here(with some shameless self promotion): [Extended Cox model and cox.zph