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?
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.
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 , 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.
 Regression Modeling Strategies. Frank Harrell
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