km.R 5.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168
  1. #' Create a KM plot
  2. #'
  3. #' @param x a data frame that contains followup and censored columns
  4. #' @param var a categorical variable to split the data to sub-curves
  5. #' @param comment An additional piece of text to write in the curve
  6. #'
  7. #' @return p probability that curves split by var differ significantly
  8. #'
  9. #' @export
  10. kaplan.meier<-function(x,var,comment=''){
  11. #x should have followup and censored columns
  12. surv.obj<-survival::Surv(x$followup,x$censored)
  13. m=base::max(x$followup)
  14. f<-stats::as.formula(paste('surv.obj',var,sep='~'))
  15. s1<-survival::survfit(f,data=x)
  16. #str(s1)
  17. tit=base::sprintf('Kaplan-Meier plot by %s',var)
  18. cols=base::c('red','blue')
  19. labels=base::c(base::sprintf('%s=true',var),base::sprintf('%s=false',var))
  20. #plot(s1,mark.time=TRUE,col=c('red','blue'),pch=labels,main=tit)
  21. graphics::plot(s1,mark.time=TRUE,col=c('red','blue'),main=tit)
  22. s=survival::survdiff(f,data=x)
  23. #str(s)
  24. p=stats::pchisq(s$chisq, length(s$n)-1, lower.tail = FALSE)
  25. sv=base::sprintf("p=%.3f",p)
  26. nLab=base::sprintf('N=%d',nrow(x))
  27. graphics::text(x=c(0.9*m,0.3*m,0.9*m),y=c(0.2,0.1,0.3),label=c(sv,comment,nLab),cex=1.2)
  28. lLab <- base::gsub("x=","",base::names(s1$strata)) ## legend labels
  29. graphics::legend("top",legend=lLab,col=cols,lty=c(1,1),horiz=FALSE, bty='n')
  30. p
  31. }
  32. set.from.list<-function(var,default,...){
  33. z<-list(...)
  34. set.from.arg.list(var,default,z)
  35. }
  36. set.from.arg.list<-function(var,default,z){
  37. if (var %in% base::names(z)) result<-base::unlist(z[[var]])
  38. else result=default
  39. result
  40. }
  41. #' Evaluate a KM plot
  42. #'
  43. #' @param x a data frame that contains followup and censored columns
  44. #' @param var a categorical variable to split the data to sub-curves
  45. #'
  46. #' @return p probability that curves split by var differ significantly
  47. #'
  48. #' @export
  49. kaplan.meier.stats<-function(x,var){
  50. #x should have followup and censored columns
  51. surv.obj<-survival::Surv(x$followup,x$censored)
  52. m=base::max(x$followup)
  53. f<-stats::as.formula(paste('surv.obj',var,sep='~'))
  54. s1<-survival::survfit(f,data=x)
  55. s=survival::survdiff(f,data=x)
  56. p=stats::pchisq(s$chisq, length(s$n)-1, lower.tail = FALSE)
  57. p
  58. }
  59. #' Plot a Kaplan-Meier curve with ggsurvfit
  60. #'
  61. #'@param x data frame containing followup and censored column
  62. #'@param var name of variable to stratify by
  63. #'@param ... other parameters:
  64. #' * varName name of variable for title (NOT GIVEN)
  65. #' * comment additional text to print on plot (empty string)
  66. #' * reorder whether to change order of categorical variable labels in legend (FALSE)
  67. #' * unit unit for time axis
  68. #' * my.legend.title title to set to legend (varName)
  69. #' * my.title title for the plot (Kaplan-Meier plot by varName)
  70. #' * my.labels labels for cases
  71. #' * draw.axis draw title and axis (FALSE)
  72. #' * my.n number of classes (2)
  73. #' * my.ylab label for y axis (Overall survival probability)
  74. #'@return graphical object
  75. #'
  76. #'@export
  77. kaplan.meier.plot.gg<-function(x,var,...){
  78. if (!requireNamespace('ggsurvfit',quiet=TRUE)){
  79. print('ggsurvfit not available. Use rNIX::kaplan.meier function')
  80. return(NULL)
  81. }
  82. #x should have followup and censored columns
  83. surv.obj<-survival::Surv(x$followup,x$censored)
  84. m=base::max(x$followup)
  85. f<-stats::as.formula(paste('surv.obj',var,sep='~'))
  86. #need survfit2
  87. s1<-ggsurvfit::survfit2(f,data=x)
  88. varName=set.from.list('varName',var,...)
  89. tit=base::sprintf('Kaplan-Meier plot by %s',varName)
  90. ylab='Overall survival probability'
  91. comment=set.from.list('comment','',...)
  92. my.labels=set.from.list('my.labels',c(),...)
  93. reorder=set.from.list('reorder',FALSE,...)
  94. unit=set.from.list('unit','day',...)
  95. my.legend.title=set.from.list('my.legend.title','NONE',...)
  96. draw.axis=set.from.list('draw.axis',FALSE,...)
  97. my.title=set.from.list('my.title',tit,...)
  98. my.n=set.from.list('my.n',2,...)
  99. my.ylab=set.from.list('my.ylab',ylab,...)
  100. base::print(base::sprintf('my.n=%f',my.n))
  101. xlab=base::sprintf('Time (%ss)',unit)
  102. if (my.n==4)
  103. cols=base::c('dodgerblue2', 'orchid2','orange','green')
  104. else
  105. cols=base::c('dodgerblue2', 'orchid2')
  106. nc=base::length(cols)
  107. base::print(base::sprintf('nc=%d',nc))
  108. labs <- base::gsub("^.*=","",base::names(s1$strata)) ## legend labels
  109. if (base::length(my.labels)==0){
  110. my.labels <- labs
  111. }
  112. nl=base::length(labs)
  113. nl1=base::length(my.labels)
  114. base::print(base::sprintf('nl=%d nl1=%d',nl,nl1))
  115. if (my.legend.title=='NONE'){
  116. my.legend.title=varName
  117. }
  118. if (reorder){
  119. perm=base::c(2,1)
  120. cols=cols[base::order(perm)]
  121. labs=labs[base::order(perm)]
  122. my.labels=my.labels[base::order(perm)]
  123. }
  124. q<-ggsurvfit::ggsurvfit(s1,lwd=1.0,censor=TRUE,censor.shape='+',censor.size=10)+
  125. ggplot2::ylim(0,1)+
  126. ggsurvfit::add_legend_title(my.legend.title)+
  127. ggsurvfit::add_pvalue("annotation", size = 5)+
  128. #add_pvalue("caption", size = 5)+
  129. ggsurvfit::add_confidence_interval()+
  130. #add_risktable()+
  131. ggplot2::scale_color_manual(values=cols,breaks = labs,labels=my.labels) +
  132. ggplot2::scale_fill_manual(values=cols, breaks = labs,labels=my.labels) +
  133. ggplot2::labs(x=NULL, y=NULL)
  134. if (draw.axis){
  135. q<-q+ggplot2::ggtitle(my.title)+ggplot2::xlab(xlab)+ggplot2::ylab(my.ylab)
  136. }
  137. q
  138. }