1

On a study is described which evaluates a protocol change in disinfectant practices in a large midwestern university medical center. Of primary interest in the study is a comparison of two methods of body cleansing. The first method, used exclusively from January 1983 to June 1984, consisted of a routine bathing care method (initial surface decontamination with 10% povidone-iodine followed with regular bathing with Dial soap). From June 1984 to the end of the study period in December1985, body cleansing was initially performedusing 4% chlorhexidine gluconate. Eighty-four patients were in the group who received the new bathing solution, chlorhexidine, and 70 patients served as the control group who received routine bathing care, povidoneiodine. Included in the data set is a covariate that measures the total surface area burned. The data is (burn). I want to test for: 1- any difference in survival functions for the two groups. 2- any difference in survival functions for the two groups adjusting for total area burned.

library(KMsurv)
data()
data(burn)
burn

library(survival)

I know that the function that would be used is coxph(), but I'm not sure which groups that I should test for (from the above information). Are they T1 and D2? so that for 1, Coxfit1<-coxph(Surv(T1,D2)~group, data = burn)? and for 2, Coxfit2<-coxph(Surv(T1,D2)~Z4, data = burn)?

What is this code doing?

for(i in 1:154){
  if (burn$??[i]==2)
    burn$Z1[i]<-1
  else burn$Z1[i]<-0
}

for(i in 1:154){
  if (burn$??[i]==3)
    burn$Z2[i]<-1
  else burn$Z2[i]<-0
}
Dana
  • 75
  • 5
  • See this: https://stackoverflow.com/questions/60593369/using-2-sided-log-rank-test-for-the-data-burn-from-kmsurv-package/60593870#60593870 and this: https://stackoverflow.com/questions/60605574/using-gehan-s-test-and-tarone-and-ware-weights/60610020#60610020 – Edward Apr 01 '20 at 02:30
  • I think you shouldn't be mixing questions. It's quite confusing now. "What is this code doing?" --- What is `burn$??[i]==2` ? That code returns an error. – Edward Apr 01 '20 at 04:16
  • There should be something here. That's why I wrote ` ??` there instead. Might be Z1? – Dana Apr 01 '20 at 04:24
  • I'm sorry but you asked what the code is doing and then you ask how to write the code. But without knowing the answer to either question, you'll just keep going around in circles. What do you _want_ the code to do? Can we go to chat? – Edward Apr 01 '20 at 04:34

2 Answers2

2

For question 1, you want to test the survival distributions* between the levels of the Z1 variable. There is no variable called group in the dataset. Z1=0 means routine bathing and Z1=1 means body cleansing. You may want to convert all Z variables to factors before continuing further (except Z4). 

library(survival)
library(KMsurv)
library (dplyr)

burn$Z1 <- factor(burn$Z1, label=c("Routine bathing", "Body cleansing"))

* The word survival needs some clarification. Presumably it is time until first straphylocous aureaus infection (D3) or on study time if no event occurred. The time is in variable T3.

The command to perform the test is:

coxph(Surv(T3,D3) ~ Z1, data=burn)
                    coef exp(coef) se(coef)      z      p
Z1Body cleansing -0.5614    0.5704   0.2934 -1.914 0.0557

For question 2, Z4 contains the percentage of total surface area burned, the variable to adjust for.

coxph(Surv(T3,D3)~Z1+Z4, data=burn)

                      coef exp(coef)  se(coef)      z     p
Z1Body cleansing -0.524764  0.591695  0.295769 -1.774 0.076
Z4                0.007248  1.007275  0.007145  1.015 0.310

So there appears to be no difference in time until first infection between those who were given routine bathing vs body cleansing.

Edward
  • 10,360
  • 2
  • 11
  • 26
  • Thanks for the clarification and being always helpful. There is a code that is used in some cases for the cox proportional hazard model. Do you know what it is this used for? Please see the code in my question. It would not be clear in the comment. – Dana Apr 01 '20 at 03:12
  • Let me check.... That's not entirely clear. `group` is not defined and neither is `myData`. Are you trying to decipher someone else's code? – Edward Apr 01 '20 at 03:16
  • I used this code for another data, which is bone marrow transplant study at Ohio State University. But I don't know what is doing. I did not modify the variables since I am not sure about it. – Dana Apr 01 '20 at 03:28
  • see this https://stackoverflow.com/questions/60919837/cox-proportional-hazard-model for clarification – Dana Apr 01 '20 at 03:32
  • Is the data in that other link the same as the burn dataset from the __KMsurv__ package? If not, why are we mixing datasets? It's confusing! – Edward Apr 01 '20 at 04:01
  • It is not the same dataset. I'm sorry if I confused you. I just gave an example to know what that code is doing. – Dana Apr 01 '20 at 04:12
  • OK - but it is fairly obvious what those code chunks do. The thing is, it has nothing to do with the burn data! Oh, you made another edit... >. – Edward Apr 01 '20 at 04:13
  • Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/210712/discussion-between-edward-and-dana). – Edward Apr 01 '20 at 04:49
1

??burn tells you what the variables mean; Z1 and Z4 seem to be what you are after:

This data frame contains the following columns:

Obs Observation number

Z1 Treatment: 0-routine bathing 1-Body cleansing

Z2 Gender (0=male 1=female)

Z3 Race: 0=nonwhite 1=white

Z4 Percentage of total surface area burned

Z5 Burn site indicator: head 1=yes, 0=no

Z6 Burn site indicator: buttock 1=yes, 0=no

Z7 Burn site indicator: trunk 1=yes, 0=no

Z8 Burn site indicator: upper leg 1=yes, 0=no

Z9 Burn site indicator: lower leg 1=yes, 0=no

Z10 Burn site indicator: respiratory tract 1=yes, 0=no

Z11 Type of burn: 1=chemical, 2=scald, 3=electric, 4=flame

T1 Time to excision or on study time

D1 Excision indicator: 1=yes 0=no

T2 Time to prophylactic antibiotic treatment or on study time

D2 Prophylactic antibiotic treatment: 1=yes 0=no

T3 Time to straphylocous aureaus infection or on study time

D3 Straphylocous aureaus infection: 1=yes 0=no

Source Klein and Moeschberger (1997) Survival Analysis Techniques for Censored and truncated data, Springer. Ichida et al. Stat. Med. 12 (1993): 301-310.

Edit: In your case, there is a significant difference between routine bathing and body cleansing (Z1), but Percentage of total surface area burned (Z4) is not significant in a univariate analysis.

library(KMsurv)
library(survival)
library(survminer)
#> Loading required package: ggplot2
#> Loading required package: ggpubr
#> Loading required package: magrittr
data(burn)

## Univariate Cox regression analysis to see whether Z1 and Z4 are significant:
res.cox <- coxph(Surv(T1, D1) ~ Z1, data = burn)
summary(res.cox)
#> Call:
#> coxph(formula = Surv(T1, D1) ~ Z1, data = burn)
#> 
#>   n= 154, number of events= 99 
#> 
#>      coef exp(coef) se(coef)     z Pr(>|z|)   
#> Z1 0.5504    1.7339   0.2072 2.656   0.0079 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>    exp(coef) exp(-coef) lower .95 upper .95
#> Z1     1.734     0.5767     1.155     2.602
#> 
#> Concordance= 0.599  (se = 0.027 )
#> Likelihood ratio test= 7.24  on 1 df,   p=0.007
#> Wald test            = 7.06  on 1 df,   p=0.008
#> Score (logrank) test = 7.23  on 1 df,   p=0.007
ggsurvplot(surv_fit(Surv(T1, D1) ~ Z1, data = burn), data = burn,
           conf.int = TRUE, pval = TRUE)



res.cox <- coxph(Surv(T1, D1) ~ Z4, data = burn)
summary(res.cox)
#> Call:
#> coxph(formula = Surv(T1, D1) ~ Z4, data = burn)
#> 
#>   n= 154, number of events= 99 
#> 
#>         coef exp(coef)  se(coef)      z Pr(>|z|)
#> Z4 -0.005108  0.994905  0.005408 -0.945    0.345
#> 
#>    exp(coef) exp(-coef) lower .95 upper .95
#> Z4    0.9949      1.005    0.9844     1.006
#> 
#> Concordance= 0.529  (se = 0.034 )
#> Likelihood ratio test= 0.94  on 1 df,   p=0.3
#> Wald test            = 0.89  on 1 df,   p=0.3
#> Score (logrank) test = 0.89  on 1 df,   p=0.3

## Multivariate Cox regression analysis to see whether Z1 and Z4 remain significant
## here, univariate Z4 was n.s., so not that relevant...
res.cox <- coxph(Surv(T1, D1) ~ Z1 + Z4, data = burn)
summary(res.cox)
#> Call:
#> coxph(formula = Surv(T1, D1) ~ Z1 + Z4, data = burn)
#> 
#>   n= 154, number of events= 99 
#> 
#>         coef exp(coef)  se(coef)      z Pr(>|z|)  
#> Z1  0.534232  1.706138  0.208651  2.560   0.0105 *
#> Z4 -0.003458  0.996548  0.005435 -0.636   0.5246  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>    exp(coef) exp(-coef) lower .95 upper .95
#> Z1    1.7061     0.5861     1.133     2.568
#> Z4    0.9965     1.0035     0.986     1.007
#> 
#> Concordance= 0.606  (se = 0.033 )
#> Likelihood ratio test= 7.66  on 2 df,   p=0neither.02
#> Wald test            = 7.44  on 2 df,   p=0.02
#> Score (logrank) test = 7.61  on 2 df,   p=0.02
user12728748
  • 8,106
  • 2
  • 9
  • 14
  • I have looked into that before, but I do not know what the 2 groups are to be tested. – Dana Apr 01 '20 at 02:16