I'm running a coxph right now. My setup: I have a reference (no treatment), and then three different treatments (A, B, and C). I also have the interactions of A, B, and C, (e.g. samples treated with both treatment A and B, or A and C, etc...). I created dummy variables for these treatments coded as 1 or 2 (1 = received treatment, 2 = did not receive treatment). I use as.factor() to load these variables.
example:
A<-as.factor(Data$A)
I can run this as follows, and get a result showing that receiving treatment B (aka B = 1) is beneficial to Lifespan (coef is positive). All three are significant in some way:
> coxph1<-coxph(Surv(Lifespan,Status)~A+B+C
> summary(coxph1)
Call:
coxph(formula = Surv(Life, Status) ~ A + B + C, data = FlyData,
method = "efron")
n= 162, number of events= 140
coef exp(coef) se(coef) z Pr(>|z|)
A -0.3486 0.7057 0.1761 -1.980 0.047753 *
B 0.5911 1.8059 0.1787 3.307 0.000944 ***
C -0.6956 0.4988 0.1815 -3.832 0.000127 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
A 0.7057 1.4170 0.4997 0.9966
B 1.8059 0.5537 1.2722 2.5635
C 0.4988 2.0050 0.3494 0.7119
Concordance= 0.822 (se = 0.095 )
Rsquare= 0.227 (max possible= 1 )
Likelihood ratio test= 41.75 on 3 df, p=5e-09
Wald test = 41.35 on 3 df, p=6e-09
Score (logrank) test = 43.6 on 3 df, p=2e-09
But when I run a coxph with interaction terms, where I want to know if A:B or A:C etc... have some interaction different from just A or B, I get the following:
> int.coxph <- coxph(Surv(Life, Status)~A*B*C, data=FlyData, method='efron')
Warning message: In fitter(X, Y, strats, offset, init, control, weights = weights, : Loglik converged before variable 1,2,3,4,5,6,7 ; beta may be infinite.
> summary(int.coxph)
Call:
coxph(formula = Surv(Life, Status) ~ A * B * C, data = FlyData,
method = "efron")
n= 162, number of events= 140
coef exp(coef) se(coef) z Pr(>|z|)
A 3.987e+01 2.066e+17 4.945e+03 0.008 0.994
B 1.856e+01 1.148e+08 2.472e+03 0.008 0.994
C 3.799e+01 3.144e+16 4.945e+03 0.008 0.994
A:B -1.964e+01 2.967e-09 2.472e+03 -0.008 0.994
A:C -3.954e+01 6.737e-18 4.945e+03 -0.008 0.994
B:C -1.874e+01 7.241e-09 2.472e+03 -0.008 0.994
A:B:C 1.962e+01 3.318e+08 2.472e+03 0.008 0.994
exp(coef) exp(-coef) lower .95 upper .95
A 2.066e+17 4.841e-18 0 Inf
B 1.148e+08 8.714e-09 0 Inf
C 3.144e+16 3.180e-17 0 Inf
A:B 2.967e-09 3.370e+08 0 Inf
A:C 6.737e-18 1.484e+17 0 Inf
B:C 7.241e-09 1.381e+08 0 Inf
A:B:C 3.318e+08 3.014e-09 0 Inf
Concordance= 0.869 (se = 0.095 )
Rsquare= 0.51 (max possible= 1 )
Likelihood ratio test= 115.6 on 7 df, p=<2e-16
Wald test = 9.24 on 7 df, p=0.2
Score (logrank) test = 73.69 on 7 df, p=3e-13
So... this is similar to some other questions... but why does beta approach infinite? The added twist I have for this question is that if I recode the variables as 0 or 1 (instead of 1 and 2), then I can change the output in the interaction coxph(). This recoding for the coxph:
coxph2<-coxph(Surv(Lifespan, Status)~A2+B2+C2))
summary(coxph2)
Call:
coxph(formula = Surv(Life, Status) ~ A2 + B2 + C2, data = FlyData,
method = "efron")
n= 162, number of events= 140
coef exp(coef) se(coef) z Pr(>|z|)
A2 0.3486 1.4170 0.1761 1.980 0.047753 *
B2 -0.5911 0.5537 0.1787 -3.307 0.000944 ***
C2 0.6956 2.0050 0.1815 3.832 0.000127 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
A2 1.4170 0.7057 1.0035 2.001
B2 0.5537 1.8059 0.3901 0.786
C2 2.0050 0.4988 1.4048 2.862
Concordance= 0.822 (se = 0.095 )
Rsquare= 0.227 (max possible= 1 )
Likelihood ratio test= 41.75 on 3 df, p=5e-09
Wald test = 41.35 on 3 df, p=6e-09
Score (logrank) test = 43.6 on 3 df, p=2e-09
is just the inverse, but the interaction coxph is different...
> full.coxph <- coxph(Surv(Life, Status)~A2*B2*C2, data=FlyData, method='efron')
Warning message:
In fitter(X, Y, strats, offset, init, control, weights = weights, :
Loglik converged before variable 2,4,6,7 ; beta may be infinite.
> summary(full.coxph)
Call:
coxph(formula = Surv(Life, Status) ~ A2 * B2 * C2, data = FlyData,
method = "efron")
n= 162, number of events= 140
coef exp(coef) se(coef) z Pr(>|z|)
A2 -7.067e-15 1.000e+00 3.204e-01 0.000 1.000
B2 -2.028e+01 1.558e-09 2.472e+03 -0.008 0.993
C2 9.821e-02 1.103e+00 3.204e-01 0.307 0.759
A2:B2 1.960e+01 3.266e+08 2.472e+03 0.008 0.994
A2:C2 -2.991e-01 7.415e-01 4.475e-01 -0.668 0.504
B2:C2 2.050e+01 7.970e+08 2.472e+03 0.008 0.993
A2:B2:C2 -1.962e+01 3.014e-09 2.472e+03 -0.008 0.994
exp(coef) exp(-coef) lower .95 upper .95
A2 1.000e+00 1.000e+00 0.5337 1.874
B2 1.558e-09 6.417e+08 0.0000 Inf
C2 1.103e+00 9.065e-01 0.5888 2.067
A2:B2 3.266e+08 3.062e-09 0.0000 Inf
A2:C2 7.415e-01 1.349e+00 0.3085 1.782
B2:C2 7.970e+08 1.255e-09 0.0000 Inf
A2:B2:C2 3.014e-09 3.318e+08 0.0000 Inf
Concordance= 0.869 (se = 0.095 )
Rsquare= 0.51 (max possible= 1 )
Likelihood ratio test= 115.6 on 7 df, p=<2e-16
Wald test = 9.24 on 7 df, p=0.2
Score (logrank) test = 73.69 on 7 df, p=3e-13
Why should changing the numerical value of a categorical variable matter? :S What am I missing here... Re-trying this with non-numeric variables ("no" and "yes") gives the same result as using 0 and 1. e.g. upper .95 for A is "1.874", for B is "inf". Similarly, coxph(Surv()~A+B+C) gives a negative coef for B, just like the above.
You probably (almost certainly in fact) have a nearly degenerate "hat matrix" which is what is formed from the model matrix with that interaction. You have all the second order interactions as well as the third order interactions. Depending on the number of levels in the factors the numbers of terms required to fully populate the model matrix might be very large. What I would try next is a model with slightly fewer terms in the model. You can use R's formula interface to remove the third order terms and only leave the first and second terms in one of two ways:
int.coxph <- coxph(Surv(Life, Status)~( A+B+C)^2, data=FlyData, method='efron')
Or:
int.coxph <- coxph(Surv(Life, Status)~ A*B*C - A:B:C, data=FlyData, method='efron')
It's not certain you will get satisfaction this way. It's possible that you don't have sufficient data to avoid the degeneracy in constructing the XX^t matrix but if your results do not blow up in as obvious a manner as is seen above, then perhaps the results will be meaningful. Another safer method would be to look first at the reduced model and then add back in specific interactions:
int.coxph.base <- coxph(Surv(Life, Status)~A+B+C, data=FlyData, method='efron')
int.coxph.intAB <- coxph(Surv(Life, Status)~A+B+C +A:B, data=FlyData, method='efron')
This second option would have the added advantage of allowing you to easily construct tests based on the change in log-likelihood rather than depending on the less reliable Wald-type tests that you see in the default printouts for print.coxph or summary.coxph.