2
$\begingroup$

We made a cox regression and ended up with huge HR for some of our variables. One of them (which was an interaction) gave us a HR=2747,093. Our dataset consists of only 74 observations. What is wrong?

$\endgroup$
2
  • $\begingroup$ How many explanatory variables (independent variables, features) do you have? $\endgroup$ – Todd D Jan 6 '19 at 17:27
  • $\begingroup$ Hi Todd. Glad to hear from you! We have 6 explanatory variables (3 categorical and 3 continuous) and three interactions. Furthermore we found that two of our continuous variables was described best cubic. Thank you for your reply! $\endgroup$ – Katrine Jan 6 '19 at 17:40
6
$\begingroup$

Most likely you have a situation where, for a particular level of a categorical variable, >0 people have an event, but for the other levels, no one does. That makes it so that the HR is infinite--the extremely large HR you find is the computer's attempt to tell you the MLE is infinity. In logistic regression, this is called "complete separation"; I assume this term still applies for survival analysis but, regardless of the name, you get the idea.

The root problem is when there is complete separation is that your model is overfitting the data. Evidently, for some level of your categorical predictor, there's a very low chance of an event, so you'd need a much larger sample size to be able to estimate the HR. It's analogous to a contingency table with some very small expected counts.

$\endgroup$
2
  • $\begingroup$ Thank you for your answer :-) Very helpful! We removed our interactions and now we have nice HR. $\endgroup$ – Katrine Jan 9 '19 at 14:00
  • $\begingroup$ Katrine - but are your HR now valid? If you removed all interactions without a research on that, and they were important, your main effects are now masked by those unhandled interactions. Maybe you should verify which interactions should remain by simple graphical analysis of predictions using these variables. Very often main effects are both non-significant and meaningless alone (itself) because there is a strong interaction between them (or them - and other covariates), which truly reveals the reality. Ignoring this leads to looking at distorted picture of it. $\endgroup$ – GibbsSampler10 Oct 7 '20 at 17:22
5
$\begingroup$

Based on the information provided in the question and the comment, your model clearly suffers strong overfitting, and possibly as a consequence, potentially strong multicollinearity.

Regarding overfitting, it is not the number of covariables that count but rather the number of parameters estimated. In your case, we have

  • 7 parameters for the continuous variables (2 cubic, 1 linear)

  • $\ge$ 3 parameters for the categorical variables. A factor with 10 levels is represented by 9 dummy variables resp. 9 parameters. So in your case, we have at least 3, but it could be a lot more.

  • $\ge$ 3 parameters representing the three interactions. Could be much more parameters, depending on whether they involve categoricals with > 2 factor levels or the cubic polynomials.

13+ parameters are too much for such a small data set. Following [1], the number of parameters should not be much larger than the number of events divided by 10. In my view, since you seem to interprete the coefficients (rather than just performing some tests), this simple rule is even too optimistic.

Now, if you fit a too complex model to a too small data set, one consequence is that some regressors are strongly correlated by chance, which can cause very large and unstable estimates of the coefficients. The same effect happens e.g. if you would use $x, x^2, x^2$ instead of orthogonal polynomials as regressors. They will be so strongly correlated that some coefficients would be very large and positive, while the other would be very large and negative.

Let's study an example with around 120 events.

library(survival)
library(car) # for function vif

options(scipen = 4)
head(veteran)

# 6 covariables with a total of 8 parameters. Should be ok on a data set with 120 events
fit <- coxph(Surv(time, status) ~ ., data = veteran)
round(exp(coef(fit)), 3)

# Output: a lot of normal sized relative effects
          trt celltypesmallcell     celltypeadeno     celltypelarge 
         1.34              2.37              3.31              1.49 
        karno          diagtime               age             prior 
         0.97              1.00              0.99              1.01 

# Variable inflation factors. All values are close to 1, so no covariable is strongly correlated with the others.
vif(fit)

# 6 covariables with a total of 17 parameters.
overfit <- coxph(Surv(time, status) ~ trt + karno + I(karno^2) + I(karno^3) +
                   (diagtime + age + celltype)^2 + prior, 
                 data = veteran)
round(exp(coef(overfit)), 2) # The coefficient of celltype adeno is huge: 381.12
vif(overfit) # many values are above 4, indicating that covariables are very strongly correlated, which basically makes it impossible to study individual effects.

[1] Regression Modeling Strategies. Frank Harrell

$\endgroup$
1
  • $\begingroup$ Thank you for your answer!! We removed our interactions and now we have nice HR :-) $\endgroup$ – Katrine Jan 9 '19 at 14:00
3
$\begingroup$

Your data likely exhibit the "curse of high dimensionality." Essentially, you have a table with 7 dimensions and additional dimensions for the interaction. Each cell in this table has very few members. If all members survived or died, the model generates unbelievable estimates. Ultimately, you have too few observations for the number of variables and interactions used to model the data.

Generating some random data, we can see that estimation of unbeleivable hazard ratios occurs very quickly (it only took 4 random seeds of sys.time() to get these estimates).

require(survival)

set.seed("2019-01-06 18:46:39 GMT")

time.days <- round(rnorm(100, 100, 20),0)
event <- rbinom(length(time.days), 1, 0.33)

x1 <- rbinom(length(time.days), 1, runif(1))
x2 <- rbinom(length(time.days), 1, runif(1))
x3 <- rbinom(length(time.days), 1, runif(1))
x4 <- rbinom(length(time.days), 1, runif(1))
x5 <- rbinom(length(time.days), 1, runif(1))
x6 <- rbinom(length(time.days), 1, runif(1))
x7 <- rbinom(length(time.days), 1, runif(1))
x8 <- rbinom(length(time.days), 1, runif(1))
x9 <- rbinom(length(time.days), 1, runif(1))

df <- data.frame(time.days, event, x1, x2,x3,x4,x5,x6,x7,x8,x9)

surv.set <- Surv(df$time.days, df$event)

summary(coxph(surv.set ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x1:x2))

      exp(coef) exp(-coef) lower .95 upper .95
x1    9.709e-01  1.030e+00    0.2676     3.523
x2    1.017e+00  9.837e-01    0.4181     2.472
x3    1.101e+00  9.085e-01    0.5443     2.226
x4    6.747e-01  1.482e+00    0.3211     1.417
x5    1.814e+00  5.512e-01    0.8956     3.675
x6    1.459e+00  6.855e-01    0.5833     3.649
x7    1.699e+00  5.885e-01    0.8583     3.364
x8    1.085e+00  9.217e-01    0.5154     2.284
x9    9.956e-01  1.004e+00    0.4903     2.022
x1:x2 8.289e-06  1.206e+05    0.0000       Inf
$\endgroup$
5
  • $\begingroup$ Thank you for a very helpful answer! We removed all our interactions and now our HR looks better :-) $\endgroup$ – Katrine Jan 9 '19 at 13:54
  • $\begingroup$ That may be a solution. Consider using cross validation to further confirm the stability of your HRs. $\endgroup$ – Todd D Jan 9 '19 at 15:03
  • $\begingroup$ We made a Log Likelihood test - would that be a way to validate the stability of our HRs? $\endgroup$ – Katrine Jan 9 '19 at 15:05
  • $\begingroup$ That assesses the condition with all data present and is a test of significance. Cross validation will assess whether the estimates remain the same when the sample is changed. $\endgroup$ – Todd D Jan 9 '19 at 15:06
  • $\begingroup$ @Katrine How did you remove the interactions? $\endgroup$ – Snorrlaxxx Mar 12 at 17:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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