4
$\begingroup$

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!

Improve this question
$\endgroup$

3 Answers 3

Reset to default
3
$\begingroup$

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.)

Improve this answer
$\endgroup$
4
  • $\begingroup$ DWin, would the specific predict.coxph() and anova.coxph() functions accomplish what you mean? $\endgroup$
    – user21068
    Commented Feb 26, 2013 at 22:32
  • $\begingroup$ 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. $\endgroup$
    – DWin
    Commented Feb 26, 2013 at 22:45
  • $\begingroup$ yes, not all possible permutations were actually performed as these experiments use live animals $\endgroup$
    – user21068
    Commented Feb 26, 2013 at 23:10
  • $\begingroup$ another set of questions! $\endgroup$
    – user21068
    Commented Feb 28, 2013 at 18:41
1
$\begingroup$

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.

Improve this answer
$\endgroup$
1
$\begingroup$

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!!

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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