km.R 4.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
  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. if (var %in% base::names(z)) result<-base::unlist(z[[var]])
  35. else result=default
  36. result
  37. }
  38. #' Plot a Kaplan-Meier curve with ggsurvfit
  39. #'
  40. #'@param x data frame containing followup and censored column
  41. #'@param var name of variable to stratify by
  42. #'@param ... other parameters:
  43. #' * varName name of variable for title (NOT GIVEN)
  44. #' * comment additional text to print on plot (empty string)
  45. #' * reorder whether to change order of categorical variable labels in legend (FALSE)
  46. #' * unit unit for time axis
  47. #' * my.legend.title title to set to legend (varName)
  48. #' * my.title title for the plot (Kaplan-Meier plot by varName)
  49. #' * my.labels labels for cases
  50. #' * draw.axis draw title and axis (FALSE)
  51. #' * my.n number of classes (2)
  52. #'@return graphical object
  53. #'
  54. #'@export
  55. kaplan.meier.plot.gg<-function(x,var,...){
  56. if (!requireNamespace('ggsurvfit',quiet=TRUE)){
  57. print('ggsurvfit not available. Use rNIX::kaplan.meier function')
  58. return(NULL)
  59. }
  60. #x should have followup and censored columns
  61. surv.obj<-survival::Surv(x$.followup,x$.censored)
  62. m=base::max(x$.followup)
  63. f<-stats::as.formula(paste('surv.obj',var,sep='~'))
  64. #need survfit2
  65. s1<-ggsurvfit::survfit2(f,data=x)
  66. varName=set.from.list('varName',var,...)
  67. tit=base::sprintf('Kaplan-Meier plot by %s',varName)
  68. comment=set.from.list('comment','',...)
  69. my.labels=set.from.list('my.labels',c(),...)
  70. reorder=set.from.list('reorder',FALSE,...)
  71. unit=set.from.list('unit','day',...)
  72. my.legend.title=set.from.list('my.legend.title','NONE',...)
  73. draw.axis=set.from.list('draw.axis',FALSE,...)
  74. my.title=set.from.list('my.title',tit,...)
  75. my.n=set.from.list('my.n',2,...)
  76. base::print(base::sprintf('my.n=%f',my.n))
  77. xlab=base::sprintf('Time (%ss)',unit)
  78. ylab='Overall survival probability'
  79. if (my.n==4)
  80. cols=base::c('dodgerblue2', 'orchid2','orange','green')
  81. else
  82. cols=base::c('dodgerblue2', 'orchid2')
  83. nc=base::length(cols)
  84. base::print(base::sprintf('nc=%d',nc))
  85. labs <- base::gsub("^.*=","",base::names(s1$strata)) ## legend labels
  86. if (base::length(my.labels)==0){
  87. my.labels <- labs
  88. }
  89. nl=base::length(labs)
  90. nl1=base::length(my.labels)
  91. base::print(base::sprintf('nl=%d nl1=%d',nl,nl1))
  92. if (my.legend.title=='NONE'){
  93. my.legend.title=varName
  94. }
  95. if (reorder){
  96. perm=base::c(2,1)
  97. cols=cols[base::order(perm)]
  98. labs=labs[base::order(perm)]
  99. my.labels=my.labels[base::order(perm)]
  100. }
  101. q<-ggsurvfit::ggsurvfit(s1,lwd=1.0,censor=TRUE,censor.shape='+',censor.size=10)+
  102. ggplot2::ylim(0,1)+
  103. ggsurvfit::add_legend_title(my.legend.title)+
  104. ggsurvfit::add_pvalue("annotation", size = 5)+
  105. #add_pvalue("caption", size = 5)+
  106. ggsurvfit::add_confidence_interval()+
  107. #add_risktable()+
  108. ggplot2::scale_color_manual(values=cols,breaks = labs,labels=my.labels) +
  109. ggplot2::scale_fill_manual(values=cols, breaks = labs,labels=my.labels) +
  110. ggplot2::labs(x=NULL, y=NULL)
  111. if (draw.axis){
  112. q<-q+ggplot2::ggtitle(my.title)+ggplot2::xlab(xlab)+ggplot2::ylab(ylab)
  113. }
  114. q
  115. }