1
$\begingroup$

The question I'm asking is related to this question here and here. And I apologize for asking so many questions here, as I am thoroughly a stats novice and probably my MD thinking is clouding my ability to interpret my findings.

So, because of the input I've received and following some contemplation, I've decided to run a Cox proportional hazards model on this dataset regarding the treatment of patients and how this affects their time of wound healing. As the duration "wound healing" is affected by censoring, following @Glen_b and @Frank Harrell's input, I believe running a normal linear and/or Poisson regression model (as the days are counted and the Poisson's residual qqplots actually have a more linear/normal distribution) isn't quite feasible, because neither of these models actually account for the fact, that when patients are healed, they drop out of the "study" and thus are "censored". Therefore I've decided to explore Cox proportional hazards modelling a bit further.

Although not originally planned from a study endpoint perspective (long story), I've decided to chose the need for operative intervention (because the wound fails to heal) as the time to event here (also because I need an event variable in the coxph and Surv function in R). The study sample size is 4918 patients, of which 575 required operative treatment (11.7%). Now I'm not even sure if this is a high enough event rate as the independent variables are n = 14 and if I go by the "rule of thumb" of 10 events per predictor variable, I'd actually require 1400 events (although I read here, that this rule can sometimes be relaxed), but that's also one of the reasons why I'm posting my question here.

I've run a multivariable Cox proportional hazards model, adjusting for those variables that may influence the decision to "treat" a patient and or/wound healing dynamics as well as the treatment the patients received.

The output looks as follows:

library(survival)
multicoxph<-coxph(Surv(Timetoheal,Operation=="Yes")~FA
+Age+Woundsurface+Mechanism+Sup+Mid+Deep)
summary(multicoxph)
Call:
coxph(formula = Surv(Timetoheal, Operation == "Yes") ~ FA + Age + 
    Woundsurface + Mechanism + Sup + Mid + Deep)

  n= 4877, number of events= 575 
   (41 observations deleted due to missingness)

                         coef exp(coef)  se(coef)      z Pr(>|z|)    
FAInadequate        -0.225084  0.798449  0.089056 -2.527   0.0115 *  
Age                 -0.012107  0.987966  0.002401 -5.042 4.60e-07 ***
Woundsurface         0.042953  1.043889  0.017338  2.477   0.0132 *  
Mechanism1          -0.341706  0.710557  0.168422 -2.029   0.0425 *  
Mechanism2          -0.169782  0.843848  0.312529 -0.543   0.5870    
Mechanism3          -0.421703  0.655929  0.523798 -0.805   0.4208    
Mechanism4          -0.131178  0.877062  0.175603 -0.747   0.4551    
Mechanism5           0.071292  1.073895  0.313708  0.227   0.8202    
Mechanism6          -1.132178  0.322331  0.473741 -2.390   0.0169 *  
Mechanism7          -1.216408  0.296292  0.308236 -3.946 7.94e-05 ***
Mechanism8          -0.011187  0.988876  0.162881 -0.069   0.9452    
SupYes              -1.476173  0.228511  1.121847 -1.316   0.1882    
MidYes              -0.056176  0.945373  1.010415 -0.056   0.9557    
DeepYes              1.568644  4.800135  1.002629  1.565   0.1177    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                    exp(coef) exp(-coef) lower .95 upper .95
FAInadequate           0.7984     1.2524   0.67057    0.9507
Age                    0.9880     1.0122   0.98333    0.9926
Woundsurface           1.0439     0.9580   1.00901    1.0800
Mechanism1             0.7106     1.4073   0.51079    0.9885
Mechanism2             0.8438     1.1850   0.45734    1.5570
Mechanism3             0.6559     1.5246   0.23496    1.8311
Mechanism4             0.8771     1.1402   0.62167    1.2374
Mechanism5             1.0739     0.9312   0.58067    1.9861
Mechanism6             0.3223     3.1024   0.12737    0.8157
Mechanism7             0.2963     3.3750   0.16194    0.5421
Mechanism8             0.9889     1.0112   0.71862    1.3608
SupYes                 0.2285     4.3762   0.02535    2.0598
MidYes                 0.9454     1.0578   0.13048    6.8497
DeepYes                4.8001     0.2083   0.67269   34.2526

Concordance= 0.746  (se = 0.017 )
Rsquare= 0.066   (max possible= 0.769 )
Likelihood ratio test= 331.9  on 14 df,   p=0
Wald test            = 240.3  on 14 df,   p=0
Score (logrank) test = 288.5  on 14 df,   p=0

Now if I understand my results correctly, this is telling me that the variable FAInadequate reduces the hazards for the event which in this case is "requirement for an operation", right? Now this is a bit bizarre, because if I perform simple comparative statistics, I find that FAInadequate patients have longer healing times and more operations as shown by these results:

t.test(log(Timetoheal+0.5)~FA, var.equal=T)

    Two Sample t-test

data:  log(Timetoheal + 0.5) by FA
t = -7.9285, df = 4916, p-value = 2.724e-15
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.2773227 -0.1673666
sample estimates:
  mean in group Adequate mean in group Inadequate 
                2.338534                 2.560879 

and this

require(gmodels)
CrossTable(Operation, FA, format="SPSS", chisq=T)

   Cell Contents
|-------------------------|
|                   Count |
| Chi-square contribution |
|             Row Percent |
|          Column Percent |
|           Total Percent |
|-------------------------|

Total Observations in Table:  4877 

             | FA 
    Operation|   Adequate | Inadequate |  Row Total | 
-------------|------------|------------|------------|
          No |      2583  |      1719  |      4302  | 
             |     2.764  |     3.835  |            | 
             |    60.042% |    39.958% |    88.210% | 
             |    91.143% |    84.141% |            | 
             |    52.963% |    35.247% |            | 
-------------|------------|------------|------------|
         Yes |       251  |       324  |       575  | 
             |    20.682  |    28.690  |            | 
             |    43.652% |    56.348% |    11.790% | 
             |     8.857% |    15.859% |            | 
             |     5.147% |     6.643% |            | 
-------------|------------|------------|------------|
Column Total |      2834  |      2043  |      4877  | 
             |    58.109% |    41.891% |            | 
-------------|------------|------------|------------|


Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  55.971     d.f. =  1     p =  7.354796e-14 

Pearson's Chi-squared test with Yates' continuity correction 
------------------------------------------------------------
Chi^2 =  55.29973     d.f. =  1     p =  1.03483e-13 


       Minimum expected frequency: 240.8704 

The fact of the matter is, that in plain "medical" terms, the longer it takes a wound to heal, the worse the outcome for the patient, so I would've originally expected the inadequate treatment to increase the hazards, especially in light of the simple t-test cross-tabulation results. Or could it be, that in light of this, running a proportional hazards function is incorrect and I should be looking at cumulative hazards? If that is the case how can I run a cumulative hazards function in R? Further, I'm a bit concerned, that the likelihood ratio, logrank and Wald test all resulted in p-values = 0. Finally, how would you ideally chose to cross-validate this model and adjust for the slightly "too low" event per variable rate?

I'm sorry for all the questions, but this has been thoroughly been driving me nuts over the last couple of days, as I can't really wrap my head around the output for some reason.

I appreciate any help greatly.

Thanks.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.