On Jun 28, 2011, at 5:46 PM, Trey Batey wrote:
> 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:
>
I modified your function to have named parameters:
> ############################################################
> siler<-function(a1=0.1, b1=0.5,a2=0.001,a3=0.003,b3=0.05) # model
> Siler parameters
> {
> 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) # 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))
cbind outputs a matrix and the presence of character values means all
your numbers are now character.... bad programming practice. Fix is
here:
> sil.anthro[, 2:6] <- sapply(sil.anthro[, 2:6], function(x)
as.numeric(as.character(x)) )
> str(sil.anthro)
'data.frame': 2 obs. of 6 variables:
$ names: Factor w/ 2 levels "Ache","Hadza": 2 1
$ a1 : num 0.351 0.157
$ b1 : num 0.895 0.721
$ a2 : num 0.011 0.013
$ a3 : num 6.7e-06 4.8e-05
$ b3 : num 0.125 0.103
> with( sil.anthro[1, -1], siler(a1=a1,b1=b1,a2=a2,a3=a3,b3=b3) )
> with( sil.anthro[2, -1], siler(a1=a1,b1=b1,a2=a2,a3=a3,b3=b3) )
Seems to work.
--
David.
>
> 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.
Could use names as index
with( sil.anthro[sil.anthro$names=="Hadza", ],
siler(a1=a1,b1=b1,a2=a2,a3=a3,b3=b3) )
--
David Winsemius, MD
West Hartford, CT