sequential treatment of a vector for formula

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

sequential treatment of a vector for formula

Frostygoat
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.
Reply | Threaded
Open this post in threaded view
|

Re: sequential treatment of a vector for formula

David Winsemius

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.
Reply | Threaded
Open this post in threaded view
|

Re: sequential treatment of a vector for formula

David Winsemius

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.
Reply | Threaded
Open this post in threaded view
|

Re: sequential treatment of a vector for formula

Wu Gong
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[]))
Reply | Threaded
Open this post in threaded view
|

Re: sequential treatment of a vector for formula

David Winsemius

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.
Reply | Threaded
Open this post in threaded view
|

Re: sequential treatment of a vector for formula

Frostygoat
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.
Reply | Threaded
Open this post in threaded view
|

Re: sequential treatment of a vector for formula

David Winsemius

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.
Reply | Threaded
Open this post in threaded view
|

Re: sequential treatment of a vector for formula

Bert Gunter
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.
Reply | Threaded
Open this post in threaded view
|

Re: sequential treatment of a vector for formula

Wu Gong
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")