Please pardon the simplicity of this question of biological nature.
I'm trying to calculate a statistic, px, the proportion of a cohort that survives through the interval x:x+1. I have the vector from which the calc is to be made but I can't figure out how to tell R to take the current value and divide it by the next value. The formula is P0=L1/LO The following is an example of the lifetable I'm constructing: example=as.vector(c(8,2,1,5,6,7,7,0,8,10,13,8,11,11,11,2,7,1,5,6,8,6)) #prime k=1 #Deaths# deaths=numeric(k) for(k in 0:max(example)) { deaths[k]=sum(example==k)} #adjust for no zero!!# deaths=c(0,deaths) #Alive, Kx# alive=sum(deaths)-cumsum(deaths) #Day, or age class,x# day=seq(from=0,to=length(deaths)-1) #mortality, lx# lx=alive/sum(deaths) #proportion surviving to next interval, px# #####Here's where I'm having trouble!!### p=1 px=numeric(p) for(p in 0:length(lx)) { px[p]=lx/lx} #create data frame# iC=data.frame("x"=day,"Kx"=alive,"Dx"=deaths,"Lx"=lx,"Px"=px) print(iC) Thanks for any suggestions! ______________________________________________ [hidden email] 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. |
On May 26, 2010, at 6:24 PM, Frostygoat wrote: > Please pardon the simplicity of this question of biological nature. > I'm trying to calculate a statistic, px, the proportion of a cohort > that survives through the interval x:x+1. I have the vector from > which the calc is to be made but I can't figure out how to tell R to > take the current value and divide it by the next value. > > The formula is P0=L1/LO > > The following is an example of the lifetable I'm constructing: > > example=as.vector(c(8,2,1,5,6,7,7,0,8,10,13,8,11,11,11,2,7,1,5,6,8,6)) > #prime > k=1 > #Deaths# > deaths=numeric(k) > for(k in 0:max(example)) > { > deaths[k]=sum(example==k)} > #adjust for no zero!!# > deaths=c(0,deaths) > #Alive, Kx# > alive=sum(deaths)-cumsum(deaths) > #Day, or age class,x# > day=seq(from=0,to=length(deaths)-1) > #mortality, lx# > lx=alive/sum(deaths) > #proportion surviving to next interval, px# > #####Here's where I'm having trouble!!### > p=1 > px=numeric(p) > for(p in 0:length(lx)) > { > px[p]=lx/lx} > > #create data frame# > iC=data.frame("x"=day,"Kx"=alive,"Dx"=deaths,"Lx"=lx,"Px"=px) > print(iC) iC$Px <- with(iC, c(1, Lx[2:nrow(iC)]/Lx[1:(nrow(iC)-1)] ) ) -- David. ______________________________________________ [hidden email] 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. |
On May 26, 2010, at 10:50 PM, David Winsemius wrote: > > On May 26, 2010, at 6:24 PM, Frostygoat wrote: > >> Please pardon the simplicity of this question of biological nature. >> I'm trying to calculate a statistic, px, the proportion of a cohort >> that survives through the interval x:x+1. I have the vector from >> which the calc is to be made but I can't figure out how to tell R to >> take the current value and divide it by the next value. >> >> The formula is P0=L1/LO >> >> The following is an example of the lifetable I'm constructing: >> >> example >> =as.vector(c(8,2,1,5,6,7,7,0,8,10,13,8,11,11,11,2,7,1,5,6,8,6)) >> #prime >> k=1 >> #Deaths# >> deaths=numeric(k) >> for(k in 0:max(example)) >> { >> deaths[k]=sum(example==k)} >> #adjust for no zero!!# >> deaths=c(0,deaths) >> #Alive, Kx# >> alive=sum(deaths)-cumsum(deaths) >> #Day, or age class,x# >> day=seq(from=0,to=length(deaths)-1) >> #mortality, lx# >> lx=alive/sum(deaths) >> #proportion surviving to next interval, px# >> #####Here's where I'm having trouble!!### >> p=1 >> px=numeric(p) >> for(p in 0:length(lx)) >> { >> px[p]=lx/lx} >> >> #create data frame# >> iC=data.frame("x"=day,"Kx"=alive,"Dx"=deaths,"Lx"=lx,"Px"=px) >> print(iC) > > iC$Px <- with(iC, c(1, Lx[2:nrow(iC)]/Lx[1:(nrow(iC)-1)] ) ) > This does the same but is more compact: iC$Px2 <- with(iC, c(1, Lx[-1]/Lx[-nrow(iC)] ) ) iC > -- > David. > > ______________________________________________ > [hidden email] 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. ______________________________________________ [hidden email] 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. |
In reply to this post by Frostygoat
I don't know if my understanding of P is right.
P ?= (the number of lives at the end of the interval)/(the number of lives at the beginning of the interval) ### Compute proportion of a cohort that survives through the interval ### The formula is P0=L1/LO ## Original data is a vector of death days example=as.vector(c(8,2,1,5,6,7,7,0,8,10,13,8,11,11,11,2,7,1,5,6,8,6)) ## Count the frenquency of each death day ## Add more levels to fulfill days range tb <- table(factor(example,levels=seq(0, max(example)))) ## Output iC <- data.frame("x"=names(tb), "Kx"=sum(tb)-cumsum(tb), "Dx"=tb[], "Lx"=iC$Kx/sum(tb), "Px"=iC$Kx/(iC$Kx+tb[])) |
On May 27, 2010, at 12:42 AM, Wu Gong wrote: > > I don't know if my understanding of P is right. > > P ?= (the number of lives at the end of the interval)/(the number of > lives > at the beginning of the interval) In standard actuarial parlance one would use lower case p_x for that quantity and would reserve upper case P_x for the cumulative survival fraction. The cumulative product of the p_x', P_x, with appropriate corrections for censoring, was rebranded in biostatistical circles as the Kaplan-Meier estimator. > > ### Compute proportion of a cohort that survives through the interval > ### The formula is P0=L1/LO > ## Original data is a vector of death days > example=as.vector(c(8,2,1,5,6,7,7,0,8,10,13,8,11,11,11,2,7,1,5,6,8,6)) > > ## Count the frenquency of each death day > ## Add more levels to fulfill days range > tb <- table(factor(example,levels=seq(0, max(example)))) > > ## Output > iC <- data.frame("x"=names(tb), > "Kx"=sum(tb)-cumsum(tb), > "Dx"=tb[], > "Lx"=iC$Kx/sum(tb), > "Px"=iC$Kx/(iC$Kx+tb[])) In a fresh session that R code does not succeed, since the data.frame function does not handle recursive references. > > ----- > A R learner. David Winsemius, MD West Hartford, CT ______________________________________________ [hidden email] 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. |
Thanks very much David. Bert, I've seen the Survival package and it
does not do what I need, which is to fit Gompertz curves to survival data and compute lifetable statistics for other functions. There is the SSgompertz package for growth curves, not what I need either. ______________________________________________ [hidden email] 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. |
On May 27, 2010, at 12:32 PM, Frostygoat wrote: > Thanks very much David. Bert, I've seen the Survival package and it survival ... correct capitalization is important in R. > does not do what I need, which is to fit Gompertz curves to survival > data and compute lifetable statistics for other functions. There is > the SSgompertz package for growth curves, not what I need either. Have you reviewed the capabilities of the eha package? http://cran.r-project.org/web/packages/eha/eha.pdf I thought that I had seen examples of the use of the Design package (and presumably this would have applied to the rms package) for fitting Gompertz models in an accelerated failure time framework, but checking Harrell's RMS text fails to confirm this memory. -- David Winsemius, MD West Hartford, CT ______________________________________________ [hidden email] 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. |
Are you using R's search capabilities (to say nothing of Google's)?
RSiteSearch("Gompertz", restr="func") brings up several functions/packages that might do what you need. Bert Gunter Genentech Nonclinical Biostatistics -----Original Message----- From: [hidden email] [mailto:[hidden email]] On Behalf Of David Winsemius Sent: Thursday, May 27, 2010 10:04 AM To: Frostygoat Cc: [hidden email] Help; Bert Gunter Subject: Re: [R] sequential treatment of a vector for formula On May 27, 2010, at 12:32 PM, Frostygoat wrote: > Thanks very much David. Bert, I've seen the Survival package and it survival ... correct capitalization is important in R. > does not do what I need, which is to fit Gompertz curves to survival > data and compute lifetable statistics for other functions. There is > the SSgompertz package for growth curves, not what I need either. Have you reviewed the capabilities of the eha package? http://cran.r-project.org/web/packages/eha/eha.pdf I thought that I had seen examples of the use of the Design package (and presumably this would have applied to the rms package) for fitting Gompertz models in an accelerated failure time framework, but checking Harrell's RMS text fails to confirm this memory. -- David Winsemius, MD West Hartford, CT ______________________________________________ [hidden email] 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. ______________________________________________ [hidden email] 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. |
In reply to this post by David Winsemius
Thank you very much David. I'm sorry for this fault , hope it has not confused Frostygoat.
I was clueless of recursive reference and I didn't meet any error when I test the code. So I wonder if there are some useful tips to prevent making this kind of faults:) The revised code is followed. ### Kaplan-Meier estimate ### Compute surviving proportion of a cohort ## Original data is a vector of death days example <- as.vector(c(8,2,1,5,6,7,7,0,8,10,13,8,11,11,11,2,7,1,5,6,8,6)) ## Count the frequency of each death day tb <- table(example) ## Output ## Interval(Start-End) "Interval" ## At risk at start of interval "RSI" ## Deaths during interval "Deaths" ## At risk at end of interval "REI" ## Propoition surviving this interval "PSI" = REI/RSI ## Cumulative surviving at end of interval "CSI" = cumprod(PSI) iC <- data.frame("Interval"=as.numeric(names(tb)), "RSI"=sum(tb)-cumsum(tb)+tb[], "Deaths"=tb[], "REI"=sum(tb)-cumsum(tb), "PSI"=(sum(tb)-cumsum(tb))/(sum(tb)-cumsum(tb)+tb[]), "CSI"=cumprod((sum(tb)-cumsum(tb))/(sum(tb)-cumsum(tb)+tb[]))) ## Plot survival curve plot(iC$Interval,iC$CSI,"s") |
Free forum by Nabble | Edit this page |