km.R 5.5 KB

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