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!
sequential treatment of a vector for formula
9 messages · Frostygoat, David Winsemius, Wu Gong +1 more
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.
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.
______________________________________________ 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.
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[]))
-----
A R learner.
View this message in context: http://r.789695.n4.nabble.com/sequential-treatment-of-a-vector-for-formula-tp2232517p2232598.html Sent from the R help mailing list archive at Nabble.com.
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
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.
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
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: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of David Winsemius
Sent: Thursday, May 27, 2010 10:04 AM
To: Frostygoat
Cc: R-help at r-project.org 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 ______________________________________________ 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.
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")
-----
A R learner.
View this message in context: http://r.789695.n4.nabble.com/sequential-treatment-of-a-vector-for-formula-tp2232517p2233491.html Sent from the R help mailing list archive at Nabble.com.