[R] plotting survival curves with model parameters

[David Winsemius]

On Jun 28, 2011, at 6:26 PM, David Winsemius wrote:

>
> 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)
>>
>
> 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) )

And maybe even clearer would be:
 > row.names(sil.anthro) <- sil.anthro$names
 > with( sil.anthro["Hadza", ], siler(a1=a1,b1=b1,a2=a2,a3=a3,b3=b3) )

You probably should have done that originally as:

sil.anthro<-data.frame(a1=a1,b1=b1,a2=a2,a3=a3,b3=b3)
row.names(sil.anthro) <- names

(But even better woulbe to choose a 'name' other than 'names'.)


>
> -- 
> David Winsemius, MD
> West Hartford, CT
>
> ______________________________________________
> 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