On Jul 5, 2011, at 6:24 PM, David Winsemius wrote:
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()
Er, ... make that par(new=TRUE)
?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@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
______________________________________________
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.
David Winsemius, MD
West Hartford, CT
______________________________________________
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.
