Hi: Here's one way to put all the plots in one graph using ggplot2 and a couple of tricks using the plyr package. You could take the data frame I generate below and use it as input to lattice graphics if you prefer. For groupwise plots, either as an ensemble or as separate panels, these packages are usually less work than an equivalent base graphic. OTOH, if you want a graphic to be 'just so', it may be better to use base graphics since both ggplot2 and lattice make certain design decisions that are not always trivial to override. I concur with Bert's comments re lattice and ggplot2, though.
The first step is to produce a suitable data frame. To generate the Siler curves for each group, one can use the mapply() function in base R or the mdply() function in the plyr package. I chose the latter, but some post-processing is necessary. Once the curves are obtained, the next step is to generate a data frame that contains the curves which will subsequently be used as input to render the curves by group. library(ggplot2) # also loads package plyr time <- 0L:90L # Function to generate a Siler survival curve siler <- function(a1, b1, a2, a3, b3) { exp(-a1/b1 * (1 - exp(-b1 * time))- a2 * time + a3/b3 * (1 - exp(b3 * time))) } # See explanation below silerCurves <- as.vector(t(mdply(hazanth[, 3:7], siler)[, 6:96])) # Generate a data frame with the populations, times and curve values silerdat <- data.frame(pop = rep(hazanth[['pop']], each = length(time)), time = rep(time, nrow(hazanth)), surv = silerCurves) # Basic form of graph: ggplot(silerdat, aes(x = time, y = surv, colour = pop)) + geom_line(size = 1) + xlab('Years') + ylab('S(t)') To get the Siler curves, I used that part of the hazanth data frame containing the parameters that determine a curve. The mdply() function operates rowwise, and substitutes the values of the parameters into the calling function (siler) to generate each curve. (The order of the parameters in the input data frame have to be exactly the same as those in the calling function.) The result is a 6 x 96 data frame, but the first five columns are the input parameters which are no longer needed; moreover, we want to transpose the result (t) and then stack the columns one below the other (i.e., the first set of curve values is on top, the last on the bottom), which is what as.vector() does. The new data frame silerdat then associates a population with each curve and matches the times corresponding to the curve values. This data frame is then input into ggplot(). If ggplot2 is new to you, the basic idea is that a plot is assembled incrementally in 'layers'. The '+' sign indicates that a new layer is being added to the existing plot. 'Geoms' represent different geometric elements of a plot; in this case, geom_smooth() is a reasonable choice to fit a smooth curve through the (time, surv) pairs. Color is an 'aesthetic' that is used here to distinguish among the various populations. A default legend is produced with the plot, which can be optionally repositioned or modified. The on-line help pages for ggplot2 are found at http://had.co.nz/ggplot2/ Scroll down to 'Reference manual'. A freely available chapter of the ggplot2 book is a useful way to get started with the package; see the 'ggplot book' heading on the same page and click on the book website to gain access to it. Similar graphs using the lattice package are obtained as follows, with a couple of choices of legend: # Legend inside the plot region mykey <- list( corner = c(1, 1), title = 'Population', cex.title = 1.1, text = list(levels(silerdat$pop), cex = 0.8), lines = list(col = 1:6, lwd = 2)) xyplot(surv ~ time, data = silerdat, groups = pop, type = 'l', lwd = 2, col.line = 1:6, key = mykey, xlab = 'Years', ylab = 'S(t)') # Legend outside the plot region (similar to the ggplot() above) mykey2 <- list( space = 'right', title = 'Population', cex.title = 1.1, text = list(levels(silerdat$pop), cex = 0.8), lines = list(col = 1:6, lwd = 2)) xyplot(surv ~ time, data = silerdat, groups = pop, type = 'l', lwd = 2, col.line = 1:6, key = mykey2, xlab = 'Years', ylab = 'S(t)') HTH, Dennis On Tue, Jul 5, 2011 at 3:08 PM, Trey Batey <ekt.ba...@gmail.com> 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. > > --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@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. > ______________________________________________ R-help@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.