# Two significant covariates in coxph analysis: exploring whether the two covariates interact significantly

I am using R 2.15.2 GUI1.53 64-bit on MacOSX 10.5.8 to perform these analyses. I am a molecular biologist by training, not a statistician, so all this analysis represents my own bootstrapped learning and I may have lots of follow-up questions.

I am performing a cox proportional hazards survival analysis on a large dataset of binomial type 1 right censored data. These data come from controlled exposures fish to different doses of virus at time zero, followed by a set time of observation to monitor mortality. All experiments revealed a strong dose-response in mortality, so I am using virus dose as a stratification term in the coxph analysis. I am interested in determining if two covariates significantly contribute to the proportional hazard of mortality: fish size(age) at time of exposure and the virus type used at exposure. A priori, I have good reason to think that there IS a significant effect from virus type, based on field observations of epidemic disease. The question of whether fish age (proxy measure being size) is a significant predictor of mortality is much less clear from field observations, but my hypothesis is that younger fish are more susceptible than older fish.

when I analyse the effect of each covariate independently via coxph, I see that both have some level of significant effect:

for virus type, which is six different viruses compared to sham groups:

    coxph(formula = surv.ALL ~ survivals$virus.type + strata(survivals$virus.dose >
2000))
coef exp(coef) se(coef)      z Pr(>|z|)
survivals$virus.type002 1.5724 4.8184 0.3415 4.605 4.13e-06 *** survivals$virus.type007  2.7217   15.2067   0.3016  9.024  < 2e-16 ***
survivals$virus.type009 3.1833 24.1272 0.3046 10.453 < 2e-16 *** survivals$virus.type110  3.5779   35.7998   0.2913 12.281  < 2e-16 ***
survivals$virus.type111 3.2039 24.6284 0.2981 10.747 < 2e-16 *** survivals$virus.type139  2.8502   17.2911   0.3013  9.460  < 2e-16 ***


AND for fish size, which is a factor with levels = sm, m, lg:

coxph(formula = surv.ALL ~ survivals$fish.size + strata(survivals$virus.dose <
2000))

coef exp(coef) se(coef)     z Pr(>|z|)
survivals$fish.sizem 0.02582 1.02615 0.08565 0.301 0.763098 survivals$fish.sizesm 0.33225   1.39411  0.08914 3.727 0.000194 ***


The effects of the virus type fit my a priori assumptions in terms of which type is most virulent (e.g. 110>111>002). My a priori hypothesis about fish size was that the smallest would suffer the greatest proportional hazard, and that appears to be true. (please correct me if these interpretations are incorrect?)

Now I would like to determine if fish size and virus stock interact in their impact on the proportional hazard of mortality. The output in R has me a bit baffled in that each factor of fish size and virus type are listed singly (where their values differ from the analyses above) and then listed in combination (please note that not all size:type combinations were tested, hence the NA values).

coxph(formula = surv.ALL ~ survivals$virus.type * survivals$fish.size +
strata(survivals$virus.dose < 2000)) coef exp(coef) survivals$virus.type002                        1.7585    5.8038
survivals$virus.type007 1.9316 6.9004 survivals$virus.type009                        3.3127   27.4600
survivals$virus.type110 2.9438 18.9882 survivals$virus.type111                        2.4122   11.1587
survivals$virus.type139 2.3408 10.3895 survivals$fish.sizem                          -0.6003    0.5486
survivals$fish.sizesm 0.6900 1.9938 survivals$virus.type002:survivals$fish.sizem NA NA survivals$virus.type007:survivals$fish.sizem 0.7983 2.2218 survivals$virus.type009:survivals$fish.sizem NA NA survivals$virus.type110:survivals$fish.sizem 0.9628 2.6189 survivals$virus.type111:survivals$fish.sizem 0.8156 2.2606 survivals$virus.type139:survivals$fish.sizem 0.2687 1.3083 survivals$virus.type002:survivals$fish.sizesm NA NA survivals$virus.type007:survivals$fish.sizesm NA NA survivals$virus.type009:survivals$fish.sizesm NA NA survivals$virus.type110:survivals$fish.sizesm -0.3135 0.7309 survivals$virus.type111:survivals$fish.sizesm NA NA survivals$virus.type139:survivals$fish.sizesm NA NA se(coef) z Pr(>|z|) survivals$virus.type002                         0.5363  3.279  0.00104
survivals$virus.type007 0.4288 4.505 6.65e-06 survivals$virus.type009                         0.5137  6.448 1.13e-10
survivals$virus.type110 0.7158 4.112 3.91e-05 survivals$virus.type111                         0.4241  5.687 1.29e-08
survivals$virus.type139 0.4244 5.515 3.49e-08 survivals$fish.sizem                            0.8662 -0.693  0.48824
survivals$fish.sizesm 0.8165 0.845 0.39807 survivals$virus.type002:survivals$fish.sizem 0.0000 NA NA survivals$virus.type007:survivals$fish.sizem 0.6661 1.199 0.23072 survivals$virus.type009:survivals$fish.sizem 0.0000 NA NA survivals$virus.type110:survivals$fish.sizem 0.8706 1.106 0.26881 survivals$virus.type111:survivals$fish.sizem 0.6598 1.236 0.21637 survivals$virus.type139:survivals$fish.sizem 0.6688 0.402 0.68783 survivals$virus.type002:survivals$fish.sizesm 0.0000 NA NA survivals$virus.type007:survivals$fish.sizesm 0.0000 NA NA survivals$virus.type009:survivals$fish.sizesm 0.0000 NA NA survivals$virus.type110:survivals$fish.sizesm 0.8221 -0.381 0.70296 survivals$virus.type111:survivals$fish.sizesm 0.0000 NA NA survivals$virus.type139:survivals$fish.sizesm 0.0000 NA NA survivals$virus.type002                       **
survivals$virus.type007 *** survivals$virus.type009                       ***
survivals$virus.type110 *** survivals$virus.type111                       ***
survivals$virus.type139 *** survivals$fish.sizem
survivals$fish.sizesm survivals$virus.type002:survivals$fish.sizem survivals$virus.type007:survivals$fish.sizem survivals$virus.type009:survivals$fish.sizem survivals$virus.type110:survivals$fish.sizem survivals$virus.type111:survivals$fish.sizem survivals$virus.type139:survivals$fish.sizem survivals$virus.type002:survivals$fish.sizesm survivals$virus.type007:survivals$fish.sizesm survivals$virus.type009:survivals$fish.sizesm survivals$virus.type110:survivals$fish.sizesm survivals$virus.type111:survivals$fish.sizesm survivals$virus.type139:survivals$fish.sizesm  so why are the single covariate lines returning different values compared to the respective analyses above? how do I interpret the exp(coef) values in the interaction pairs- since none have significant p values, are there no significant interactions? beyond this question about interaction, I would like to look for signs that the underlying assumption of coxph are not being violated. Is there a specific function in the survival package for this (I really struggle to understand the R documentation language, as it assumes a much greater statistical knowledge than I have)? also, I came across source('http://www.stat.ucla.edu/~david/teac/surv/local-coxph-test.R') as a script to perform a local test for significant explanatory variables, but it is defunct as far as I can tell. is there a more up to date function for this? thank you for your attention! - add comment ## 3 Answers Q: "how do I interpret the exp(coef) values in the interaction pairs" A: You don't. You use predict(). And you don't separately examine the p-values, either. The separate p-values are pretty much meaningless. You use anova()-type analyses to compare nested models, ideally with some method to penalize the process so you are not susceptible to multiple comparisons artifacts. (Note: it is possible to compare some of the coefficients. The single level "main-effects" fishsize coefficients in an interaction model would be the expected log(hazard ratios) for the other fishsizes compared with the baseline levels but only for the reference value of virus. And similarly for the "main-effects" virus coefficients.) - thank you both. – user21068 Feb 26 '13 at 22:29 DWin, would the specific predict.coxph() and anova.coxph() functions accomplish what you mean? – user21068 Feb 26 '13 at 22:32 Yes, those functions should do quite nicely. You will notice that you have some NA's in your coefficient combination. You seem to be missing events in some of those category crosses. – DWin Feb 26 '13 at 22:45 yes, not all possible permutations were actually performed as these experiments use live animals – user21068 Feb 26 '13 at 23:10 another set of questions! – user21068 Feb 28 '13 at 18:41 add comment Only a partial answer to your questions: why are the single covariate lines returning different values compared to the respective analyses above For exactly the same reason that this happens in multiple regression. If you fit a main effects model then fit an interaction model, then in general the main-effects coefficients (and their standard errors!) will change. This change in coefficients occurs because the main effects and the interactions are not independent. - add comment DWin: can you explain this 'compared with the baseline levels but only for the reference value of virus' a little bit more for me? I was following your description to this point, then lost it also have a question about the way factor covariates are handled in coxph: In the virus.type only block above, there are actually 7 factors of virus type and the one that is not shown is the mock virus exposure control group (I think of this as the comparator). the interpretation of coef for a continuous covariate (proportional hazard for a single unit increase in the covariate) would not make sense for a factor covariate, so the individual proportional hazards are calculated for each factor value. is the comparator factor value, the one not printed, assumed to have a exp(coef) = 1? and therefore the p-value is only for the pairwise comparison of the printed factor and the non-printed comparator, right? what role does the ordering or level of the factor play in the calculation? when I alter their order, in the naive hope that I could explore the proportional hazards of one virus directly compared to another (e.g. virus type 002 is the comparator) the coef values go all wacky:  coef exp(coef) se(coef) z factor(survivals$virus.type)009  2.1170    8.3064   0.3400  6.227
factor(survivals$virus.type)007 3.2342 25.3861 0.3001 10.778 factor(survivals$virus.type)111  3.6765   39.5067   0.3031 12.128
factor(survivals$virus.type)110 4.0690 58.4997 0.2898 14.039 factor(survivals$virus.type)139  3.7135   40.9957   0.2966 12.521
factor(survivals\$virus.type)0    3.3732   29.1729   0.2997 11.254


this suggests that the mock exposure carries a 29-fold greater risk of death than exposure to the 002 virus? that seems nonsense. I suspect that I am doing something fundamentally improper here, but I would like to understand what/why, and how I might think about how the different factors compare to one another instead of just the mock control. would I have to do that with anova on the coef/se values?

thanks!!

-