[R] Package "demography" - calculating percentiles of survival probabilities distribution

[David Winsemius]

On Apr 20, 2012, at 4:13 PM, jolo999 wrote:

> Hi,
>
> I am using the package "demography" from Rob Hyndman for the
> Lee-Carter-Model. It is an amazing powerful tool but I am struggling  
> with
> one issue:
>
> I want to compute different percentiles of the survival probability
> distribution derived from the Lee-Carter-Forecast (e.g. the 50%tile,
> 60%tile, 75%tile and 99%tile) for each of the next 10 years. Is  
> there any
> possibility to retrieve these "attachment points"?
>
> I am sure the package possess this functionality but I didn't find  
> any way
> to calculate percentile values.

Even with your revised question I'm having trouble understanding  
exactly what you want. A survival function is a probability that is  
returned from input of a specified duration, i.e. a value between 0  
and 1 (or equivalently between 0 and 100%). Conversely, you could ask  
what time would be associated with a particular survival percentile  
for a particular age. But I am having trouble understanding what it  
means to ask for a "50th percentile for each of the next 10 years"?  
(It seems you are specifying both the input and the output.)

There are very few ages at which the 50th percentile would be reached  
in one year, or even ten years. Are you asking what initial ages would  
be associated with a 50% survival at 10 years when assuming a  
particular mortality structure? Can you offer a particular example for  
further discussion?

Looking at the returned object descriptions, it appears that they are  
in the form of mortality tables. Are you really asking how to produce  
survival estimates from a particular series of mortality predictions?  
Here is one formula for converting m's to S's

exp( -sum(m[i:(1+9)] ) #  for 'i' being the starting age index and m's  
being a series of annual mortalities .

(again, offering an example from the help pages might clarify this  
trans-Atlantic, trans-language exchange.)

france.LC1 <- lca(fr.mort,adjust="e0")
france.fcast <- forecast(france.LC1, 50)
str(france.fcast$rate$total[ , which(france.fcast$year == 2012)])
# num [1:101] 2.37e-03 1.04e-04 5.14e-05 4.19e-05 3.77e-05 ...
france.fcast$rate$total[ 90:101, which(france.fcast$year == 2012)]
# [1] 0.1571447 0.1774625 0.2004623 0.2261028 0.2546273 0.2838347  
0.3061775 0.3229046 0.3398724
# [10] 0.3465093 0.3362461 0.5371380
exp(- france.fcast$rate$total[ 90:101, which(france.fcast$year ==  
2012)] )
  [1] 0.8545804 0.8373924 0.8183523 0.7976361 0.7752053 0.7528911  
0.7362559 0.7240429 0.7118612
[10] 0.7071523 0.7144473 0.5844184

So at no age with assumed starting year of 2012 would there be an  
associated 50% probability.

These are the one year survival probabilities for ages 80-101 in year  
2012 (based on a forcast from 2006 data.)

 > exp(- france.fcast$rate$total[ 80:101, which(france.fcast$year ==  
2012)] )
  [1] 0.9574380 0.9530564 0.9474982 0.9411758 0.9338487 0.9254885  
0.9150893 0.9022344 0.8882844
[10] 0.8719397 0.8545804 0.8373924 0.8183523 0.7976361 0.7752053  
0.7528911 0.7362559 0.7240429
[19] 0.7118612 0.7071523 0.7144473 0.5844184

-- 

David Winsemius, MD
West Hartford, CT