[R] plotting survival curves (multiple curves on single graph)

[David Winsemius]

On Jul 5, 2011, at 6:08 PM, Trey Batey wrote:

> Hello.
>
> This is a follow-up to a question I posted last week.  With some
> previous suggestions from the R-help community, I have been able to
> plot survival (, hazard, and density) curves using published data for
> Siler hazard parameters from a number of ethnographic populations.
> Can the function below be modified, perhaps with a "for" statement, so
> that multiple curves (different line types---one for each population)
> are plotted on a single graph for comparison?  Thanks so much.
>

There are (at least) three methods to plot multiple curves in base  
plotting:
-- plot() then lines()

?lines

--plot(); par(add=TRUE); plot()

?par

# There is also matplot()
?matplot

After extracting the sil function to exist on its own, you could try:

  matplot(x=t, y=apply(hazanth[ ,3:7], 1, sil)

My first choice would be to make a modified version of your silsurv  
that uses the lines function rather than plot and then you can just  
use the lines of code you already have.

-- 
David
> --Trey
>
> The function and calls below use the data in this Excel file (feel
> free to access):
> https://docs.google.com/leaf?id=0B5zZGW2utJN0ZDk1NjA0ZjUtMWU0ZS00ZGQ3LWIxZTUtOWE0NGVmYWMxODJl&hl=en_US
>
> ## - plot Siler survival curve
> ##############################
> silsurv<-function(a1,b1,a2,a3,b3)
>  {
>    sil=function(t)
>      {
>        h.t<-a1*exp(-b1*t)+a2+a3*exp(b3*t)
>        S.t<-exp(-a1/b1*(1-exp(-b1*t))-a2*t+a3/b3*(1-exp(b3*t)))
>        d.t<-S.t*h.t
>
>        #return(d.t)
>        return(S.t)
>        #return(h.t)
>      }
>    t<-seq(0,90,1)
> plot 
> (t 
> ,sil 
> (t 
> ),ylim=c(0,1),type='l',cex.lab=0.8,cex.axis=0.75,ylab='S(t)',xlab='Age
> (years)')
>  }
>
> with 
> (hazanth 
> [1,3 
> : 
> 7 
> ],silsurv 
> (a1=a1,b1=b1,a2=a2,a3=a3,b3=b3));title(main=hazanth[1,1],cex.main=0.9)
> # plot for Hadza
> with 
> (hazanth 
> [2,3 
> : 
> 7 
> ],silsurv 
> (a1=a1,b1=b1,a2=a2,a3=a3,b3=b3));title(main=hazanth[2,1],cex.main=0.9)
> # plot for Ache
> with 
> (hazanth 
> [3,3 
> : 
> 7 
> ],silsurv 
> (a1=a1,b1=b1,a2=a2,a3=a3,b3=b3));title(main=hazanth[3,1],cex.main=0.9)
> # plot for Hiwi
> with 
> (hazanth 
> [4,3 
> : 
> 7 
> ],silsurv 
> (a1=a1,b1=b1,a2=a2,a3=a3,b3=b3));title(main=hazanth[4,1],cex.main=0.9)
> # plot for !Kung
> with 
> (hazanth 
> [5,3 
> : 
> 7 
> ],silsurv 
> (a1=a1,b1=b1,a2=a2,a3=a3,b3=b3));title(main=hazanth[5,1],cex.main=0.9)
> # plot for Yanomamo
> with 
> (hazanth 
> [6,3 
> : 
> 7 
> ],silsurv 
> (a1=a1,b1=b1,a2=a2,a3=a3,b3=b3));title(main=hazanth[6,1],cex.main=0.9)
> # plot for Tsimane
>
> ###############################
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius, MD
West Hartford, CT