0
$\begingroup$

I have data that I analyzed using a stratified PWP-GT extension of the Cox PH model and I am having trouble figuring out what I need to report in the upcome manuscript and how to interpret the stratified components. Briefly, the data contains information about how long it took individual insects to reach certain developmental points ("event": 1, 2, 3, 4) after being treated with a pesticide. Including the control group, there were four levels of pesticide ("tx": C, L, M, and H). After some preliminary data exploration, I found that when the insects arrived ("ship": 1, 2, 3, 4) and started the experiment had a significant effect on development time.

From all of this, I made the following model in R:

> coxph(formula = Surv(rep(0, dim(sub1)[1]), tstop - tstart, status) ~ 
tx * strata(event) + ship * strata(event) + cluster(id) + 
    strata(event), data = sub1, method = "breslow")

Where "id" equals the individual insects and "tstop - tstart" indicates the gap time between events.

This outputs: (EDITED TO INCLUDE "ship" AS A FACTOR)

Call:
coxph(formula = Surv(rep(0, dim(sub1)[1]), tstop - tstart, status) ~ 
tx * strata(event) + ship * strata(event) + cluster(id) + 
    strata(event), data = sub1, method = "breslow")

  n= 492, number of events= 478 

                                    coef         exp(coef)          se(coef)         robust se         z       Pr(>|z|)    
txh                        -0.26370516438318  0.76819999899147  0.26233252301991  0.22689202991138  -1.16225  0.2451341    
txl                        -0.11987649941476  0.88702997867426  0.25672484677695  0.18618575422004  -0.64385  0.5196699    
txm                        -0.51362257507084  0.59832417391156  0.27195280894235  0.23832152919253  -2.15517  0.0311488 *  
ship2                       1.65646739004232  5.24076452832491  1.03399484661244  0.21696246528283   7.63481 2.2649e-14 ***
ship3                       1.33578200803071  3.80296873598721  0.23058142160552  0.22861869093861   5.84284 5.1319e-09 ***
txh:strata(event)event=2    0.09842786823565  1.10343480882414  0.37065411069586  0.31825353892903   0.30928  0.7571123    
txl:strata(event)event=2   -0.07427724599606  0.92841425902055  0.36269132158140  0.24381582455421  -0.30464  0.7606366    
txm:strata(event)event=2    0.49030567879075  1.63281526067281  0.38153511309810  0.31897747833434   1.53712  0.1242647    
txh:strata(event)event=3    0.08024387819237  1.08355128998456  0.38622480003047  0.34278826007442   0.23409  0.8149138    
txl:strata(event)event=3   -0.05197138144960  0.94935603564767  0.37101662106664  0.32942786091746  -0.15776  0.8746439    
txm:strata(event)event=3    0.48005016860266  1.61615548042121  0.38211881930279  0.35384318204038   1.35667  0.1748845    
txh:strata(event)event=4    0.25828662043757  1.29470985575366  0.37678415307640  0.29040930717579   0.88939  0.3737945    
txl:strata(event)event=4    0.01633743201781  1.01647161761367  0.37635750078046  0.23765893060307   0.06874  0.9451940    
txm:strata(event)event=4    0.79156021665709  2.20683688529518  0.39001352229631  0.25824225217817   3.06518  0.0021754 ** 
strata(event)event=2:ship2 -2.87880676955391  0.05620178452044  1.45849731918900  0.28690571886754 -10.03398 < 2.22e-16 ***
strata(event)event=3:ship2  0.59127090608651  1.80628257277547  1.46984114123398  0.34311762620912   1.72323  0.0848468 .  
strata(event)event=4:ship2 -1.80811450547487  0.16396299686714  1.45261480585070  0.21783939462209  -8.30022 < 2.22e-16 ***
strata(event)event=2:ship3 -2.15379915896782  0.11604245552070  0.31838628524726  0.28307156268510  -7.60867 2.7645e-14 ***
strata(event)event=3:ship3  0.22972491271199  1.25825383268272  0.34233273407819  0.36139360103116   0.63566  0.5249954    
strata(event)event=4:ship3 -1.15889323548630  0.31383332833294  0.32161353666032  0.26768619263227  -4.32930 1.4959e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                              exp(coef)       exp(-coef)        lower .95        upper .95
txh                        0.76819999899147  1.3017443391211 0.49242881278360 1.19840923831125
txl                        0.88702997867426  1.1273576136565 0.61582412220287 1.27767353485976
txm                        0.59832417391156  1.6713347773707 0.37503940106341 0.95454455204408
ship2                      5.24076452832491  0.1908118547581 3.42543603154333 8.01813625723261
ship3                      3.80296873598721  0.2629524640939 2.42952947665736 5.95282804586274
txh:strata(event)event=2   1.10343480882414  0.9062610604659 0.59135568542433 2.05894423159404
txl:strata(event)event=2   0.92841425902055  1.0771053872600 0.57571215103501 1.49719444830733
txm:strata(event)event=2   1.63281526067281  0.6124391559079 0.87382190804916 3.05106298082884
txh:strata(event)event=3   1.08355128998456  0.9228912458904 0.55343621886470 2.12144301006479
txl:strata(event)event=3   0.94935603564767  1.0533455968580 0.49775956270853 1.81066713719455
txm:strata(event)event=3   1.61615548042121  0.6187523490867 0.80777661572505 3.23351590779965
txh:strata(event)event=4   1.29470985575366  0.7723738222552 0.73278336798836 2.28754320009662
txl:strata(event)event=4   1.01647161761367  0.9837953000081 0.63796904694880 1.61953711446611
txm:strata(event)event=4   2.20683688529518  0.4531372511776 1.33031359891976 3.66088796074397
strata(event)event=2:ship2 0.05620178452044 17.7930293234080 0.03202841999485 0.09861993141683
strata(event)event=3:ship2 1.80628257277547  0.5536232343002 0.92198411562698 3.53873421180760
strata(event)event=4:ship2 0.16396299686714  6.0989370718218 0.10698444024686 0.25128761041906
strata(event)event=2:ship3 0.11604245552070  8.6175356727230 0.06662940869666 0.20210072018766
strata(event)event=3:ship3 1.25825383268272  0.7947521986624 0.61965430467327 2.55497734062470
strata(event)event=4:ship3 0.31383332833294  3.1864047241633 0.18571379401475 0.53033948552421

Concordance= 0.653  (se = 0.043 )
Rsquare= 0.164   (max possible= 1 )
Likelihood ratio test= 88.11  on 20 df,   p=1.581164088549e-10
Wald test            = 404.6  on 20 df,   p=0
Score (logrank) test = 110.98  on 20 df,   p=1.298960938811e-14,   Robust = 52.59  p=9.340668759683e-05

  (Note: the likelihood ratio and score tests assume independence of
 observations within a cluster, the Wald and robust score tests do not).

I understand the individual components of the output (coef, exp(coef), se(coef), z, Pr(>|z|), etc.) and that the numbers are compared to a baseline group (I think... txC, event=1, ship=1). However, I'm still unsure about two things:

  1. What do I report in a manuscript? I have a couple other of these outputs and it would quickly get overwhelming for a reader to make their way through all of them. Would it be adequate just to explain the effects of the covariate of interest (tx)?

  2. What about comparing the other treatments to each other? For example, I'd like to know the difference between L and H. My output doesn't contain that information. Would I just rerun the model without C to set a new baseline group? I'm hesitant because I don't want to inflate my type I error.

$\endgroup$
2
  • $\begingroup$ It looks like your coding evaluated ship as a standard numeric variable rather than as a factor (which is, I think, what you would have intended). You probably want to fix that first. $\endgroup$
    – EdM
    Dec 4, 2017 at 19:17
  • $\begingroup$ Thanks @EdM! I've changed that and have updated the output in the question post. $\endgroup$ Dec 6, 2017 at 21:10

1 Answer 1

1
$\begingroup$

You need to go beyond the standard output from coxph() and analyze the data in a way that accomplishes two things. For your first question, what you presumably want is a more general idea whether any of your treatments affect developmental times significantly, including both their direct effects and their interactions with shipments and types of events. For your second question, you need flexibility to examine particular combinations (contrasts) of factor levels other than the comparisons of each individually against the reference level (as you correctly take the coxph() output to represent.)

Although it's possible to do this yourself by hand, the rms package in R provides these capabilities directly. There is a bit of learning curve in terms of pre-processing the data (datadist is important) and using the cph() function with appropriate parameter values to allow useful downstream processing, but it's well worth it. Then the anova() function applied to the cph() output provides information about whether the treatments differ among themselves (as main effects or including interaction terms) and the contrast() function allows for great flexibility in terms of examining pre-specified combinations of factor levels.

The rms package also provides ways to calibrate and validate your model. You need to know that your data adequately meet the proportional hazards assumption before you report results that depend on the assumption, and you want to make sure that your results aren't too dependent on the particular data sample you have.

If you are publishing in a journal that provides for supplemental material outside of the main text, showing the anova tables and specific contrast results in the main text and putting more detailed results in the supplement might work. In clinical survival analysis, usual practice is to present hazard ratios (exp(coef) in the output) along with their confidence intervals. As you might be more interested in whether the treatments delay rather than accelerate the developmental steps, you might consider using exp(-coef) instead as a measure of delay and recalculate the confidence intervals from -coef and robust se.

In this particular application, you will of course have to deal with the biological problem of the differences among shipments being the greatest of all influences. That makes it very hard to go from your specific data to generalizations about these types of insects and treatments. A more than 5-fold difference in rates of developmental progression between 2 shipments of insects (ship2 versus ship1) is pretty remarkable. My guess is that, in the background of those shipment differences and their striking interactions with specific developmental stages, any effects of the treatments themselves will be lost.

$\endgroup$
4
  • $\begingroup$ Thanks for your response @EdM! I've come across rms before, but I was not sure how to implement it for a multi-event model. For coxph() I had come across a tutorial by Amarim and Cai 2015 (International Journal of Epidemiology). I'll have to look into rms more and see if I can figure it out. For some of my single-event models that do not meet the PH assumption, I have used the package MRH and its estimateMRH() function. That seemed to work well, but, again, I'm not sure how to expand it to multi-event models or if it is even possible. $\endgroup$ Dec 11, 2017 at 21:44
  • $\begingroup$ Briefly, regarding ship, differences were expected and can be explained fairly satisfactorily. The insect is univoltine (1 generation per year) and individuals from eggs laid later in the year develop more quickly. $\endgroup$ Dec 11, 2017 at 21:45
  • $\begingroup$ cph in rms pretty much works like coxph except that you use strat() instead of strata(); cluster() constructs are allowed in the formula. Worth the effort to learn rms if you will continue to do regression analysis of any type. $\endgroup$
    – EdM
    Dec 12, 2017 at 1:43
  • $\begingroup$ I've located Frank Harrell Jr.'s book on regression modeling strategies. Hopefully that will give me a good grasp on the rms package. It's a little intimidating at nearly 600 pages. The R documentation, in comparison, is a svelte 246 pages. Thanks again for your help @EdM. I've had a heck of a time finding help at my own institution, and yours is greatly appreciated. $\endgroup$ Dec 12, 2017 at 2:59

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.