[R] plotting survival curves with model parameters

[Trey Batey]

Hello.

I am trying to write an R function to plot the survival function (and
associated hazard and density) for a Siler competing hazards model.
This model is similar to the Gompertz-Makeham, with the addition of a
juvenile component that includes two parameters---one that describes
the initial infant mortality rate, and a negative exponential that
describes typical mortality decline over the juvenile period.  The
entire hazard is expressed as


h(x) = a1*exp(-b1*x)+a2+a3*exp(b3*x)


I've had success in plotting the curves using the following function:

############################################################
siler<-function(x=c(0.1,0.5,0.001,0.003,0.05))   # model Siler parameters
  {
    sil=function(t)
      {
        a1 = x[1]
        b1 = x[2]
        a2 = x[3]
        a3 = x[4]
        b3 = x[5]

        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)    # returns the survival function
        #return(h.t)
      }

    t=seq(0,100,0.01)
    plot(t,sil(t),ylim=c(0,1),type='l',main="Siler model of mortality:
 Wood et al. (2002, Figure
7.4)",cex.main=0.9,cex.lab=0.8,cex.axis=0.75,ylab='S(t)',xlab='Age
(years)')    # reproduces Figure 7.4 from Wood et al. (2002)

  }

siler()
#########################################################################

How can I modify the function so that I can plot curves using
published Siler parameters I have already compiled into a dataframe
(below)?

names<-c("Hadza","Ache")
a1<-c(0.351,0.157)
b1<-c(0.895,0.721)
a2<-c(0.011,0.013)
a3<-c(0.0000067,0.000048)
b3<-c(0.125,0.103)
sil.anthro<-data.frame(cbind(names,a1,b1,a2,a3,b3))

For example, how can I modify the function above to produce a curve
for the a specific group (e.g., Hadza, Ache...) or multiple groups on
one graph?  Thanks.

--Trey