Simulando o Tempo de Vida Futuro

Documento que objetiva ilustrar o uso da função rLife() do pacote {lifecontingencies} para simular o tempo de vida futuro para um conjunto de indivíduos,

Author

Marcos F. Silva

Published

January 27, 2025

Introdução

Este documento tem por objetivo mostrar como usar a função rLife() do pacote {lifecontingencies} para simular o tempo de vida futuro (\(T(X)\)) até a morte dos indivíduos de um conjunto de dados que poderiam ser, por exemplo, os servidores aposentados de um RPPS.

Artigo publicado na Revista de Estatística UFOP, Vol VI, 2017 - ISSN 2237-8111

A inspiração para a elaboração do documento surgiu da leitura do artigo Estimação do Tempo de Vida Futuro Aleatório Individual para Grupos de Indivíduos de autoria da Prof^a Cristiane Silva Correa, no qual ela apresenta uma função desenvolvida em R para estimar tempo de vida futuro individual para grupos de indivíduos.

A autora menciona que a necessidade de criar a função decorreu do fato de que a função rLife() do pacote {lifecontingencies} não é aplicável a “uma matriz de vetores com idades diferentes ou sexos diferentes de uma só vez, de forma que, se utilizada nesse caso, teria que ser aplicado para cada célula da matriz individualmente, aumentando o tempo de processamento”. (grifo nosso)

Considerando o avanço ocorrido na linguagem R desde então, em especial o surgimento dos pacotes {purrr} e {furrr} pretendemos mostrar que a limitação apontada pela autora pode ser contornada com o uso dos pacotes em referência.

Dados utilizados

Para ilustrar o procedimento, vamos utilizar uma base de dados fictícia utilizada pela Prof^a Andréia Vanzilota em seu curso “Avaliação Atuarial em Planos de Previdência (2ª edição)” ministrado na Pós-Graduação em Atuária da UFRJ cujos vídeos e materiais estão disponíveis no youtube.

Vamos também utilizar a tábua biométrica IBGE de 2022.

A função `rLife()` e o tempo de vida futuro (\(T(X)\))

O tempo de vida futuro é um valor aleatório que representa a quantidade de anos que um indivíduo irá viver a partir da idade atual \(x\).

Essa função do pacote {lifecontingencies} simula uma quantidade \(n\) de tempo de vida futuro \(T(X)\) para um indivíduo de idade \(x\) segundo as taxas de mortalidade de uma tábua biométrica especificada.

Conforme apresentado no artigo da Prof^a Cristiane, o tempo de vida futuro (\(T(X)\)) pode ser obtido a partir dos valores de \(\ell_x\) de uma tábua biométrica e, dada a limitação então identificada na função rLife(), a autora desenvolveu três funções (EstimaTX(), procNDif()e GeraTX()) em R que, usadas em conjunto, nos fornecem a distribuição do tempo de vida futuro para um conjunto de indivíduos com idades e sexos distintos.

A função rLife() fornece a mesma funcionalidade, mas tem o inconveniente de aceitar apenas um único valor para o argumento relativo à idade para a qual se deseja simular o tempo de vida futuro.

Ocorre que, na prática, deseja-se simular o tempo de vida futuro não apenas para um único indivíduo de idade \(x\), mas para um grupo de indivíduos de diversas idades tanto do sexo masculino como do sexo feminino.

Inicialmente vamos mostrar como funciona a função rLife() para depois vermos como fazer com que ela funcione para diversas idades e para isso vamos precisar criar dois objetos da classe lifetable para representar as tábuas IBGE 2022 Masculina e Feminina. Vamos fazer isso a seguir mas, primeiro, vamos carregar os pacotes que serão necessários:

suppressPackageStartupMessages({
  library(lifecontingencies)
  library(readxl)
  library(tidyverse)
  library(knitr)
  library(tictoc)
  library(furrr)
})

Agora vamos criar os objetos da classe lifetable.

Primeiro devemos importar os dados das tábuas de mortalidade:

# MASCULINO

ibge2022M <- read_excel("tabua-ibge-2022-extrapolada-drpsp_19_12_2023.xlsx",
                        range = "A4:D115",
                        col_names = c("x", "lx", "qx", "ex"))

# FEMININO

ibge2022F <- read_excel("tabua-ibge-2022-extrapolada-drpsp_19_12_2023.xlsx",
                        range = "E4:G115",
                        col_names = c("lx", "qx", "ex")) %>% 
             mutate(x = 0:111, .before = 1)

Agora podemos criar os objetos da classe lifetable que serão usados como um dos argumentos da função rLife():

ibge2022M_ltb <- new("lifetable", lx=ibge2022M$lx, x=ibge2022M$x,  name="IBGE2022M")
ibge2022F_ltb <- new("lifetable", lx=ibge2022F$lx, x=ibge2022F$x,  name="IBGE2022F")

Com isso já estamos em condições de usar a função rLife(). Vamos gerar 10.000 simulações para o tempo de vida futuro de uma pessoa de 40 anos de idade do sexo masculino e com os resultados obtidos vamos gerar algumas estatísticas descritivas e um histograma para “vermos” a distribuição dos valores obtidos.

sim_Tx <- rLife(n = 10000, object = ibge2022M_ltb, x = 40)
summary(sim_Tx)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.50   27.50   38.50   36.16   46.50   65.50

Gerados os 10.000 valores que estão agora armazenados no vetor sim_Tx, vamos fazer um histograma para visualizarmos a distribuição dos valores. Nesse histograma vamos marcar com linhas verticais três quantidades: a expectativa de vida \(e_x\) fornecida pela tábua para a idade de 40 anos, a expectiva de vida calculada a partir dos valores simulados e a mediana dos valores simulados.

data.frame(Tx = sim_Tx) %>% 
  ggplot(aes(x = Tx)) +
  geom_histogram(fill = "lightblue", color="white") +
  geom_vline(xintercept = c(36.34, 36.54, 38.50), linewidth=0.7, color=c("tomato", "black", "orange")) +
  theme_bw() +
  labs(x = "Tempo de Vida Futuro - T(X)", y = "Qtd. de Pessoas")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

A simulação realizada produziu resultados que variam de um mínimo de 0.5 a um máximo de 65.5, significando que uma pessoa que tenha 40 anos pode sobreviver uma quantidade de anos contida nesse intervalo, mas alguns tempos de vida futuro são mais prováveis do que outros conforme nos mostra o histograma.

Por exemplo, é menos provável morrer aos 40 anos (tempo de vida futuro de 0 anos) do que aos 80 anos (tempo de vida futuro de 40 anos) e é mais provável morrer aos 80 anos do que aos 100 anos (tempo de vida futuro de 60 anos).

No gráfico marcamos em vermelho e preto, respectivamente, e expectativa de vida teórica (\(e_{40}\)) para um indivído do sexo masculino de 40 anos dada pela tábua de mortalidade (36,34) e a média dos 10.000 valores obtidos pela simulação (36,43). Os dois valores são praticamente os mesmos.

Cálculo Atuarial Aplicado. Teoria e Aplicações. Exercícios Resolvidos e Propostos, 2^a Edição. Atlas, 2019 (pág. 13).

Marcamos também no gráfico em laranja o valor referente à mediana dos 10.000 valores obtidos pela simulação (38,50), quantidade também conhecida por “vida provável” conforme nos ensina Antonio Cordeiro Filho.

A função rLife() nos permite obter essa distribuição de valores apenas para uma única idade (40 anos no nosso exemplo), ou seja, o argumento “x =” aceita apenas um único valor.

Para contornar este inconveniente vamos mostrar a seguir como obter, por exemplo, a vida média, ou expectativa de vida (\(e_x\)), e a “vida provável” a partir da simulação do tempo de vida futuro para todas as idades da tábua ou seja, para \(x = 0, 1, 2, \dots, 111\), tanto para o sexo masculino como para o feminino.

Simulação da Vida Média

O comando a seguir irá produzir as simulações com o uso da função map() e, a partir delas fará o cálculo da média dos tempos de vida futuro simulados com a função map_dbl().

# média do tempo de vida futuro

vida_media_masc <- map(0:getOmega(ibge2022M_ltb), \(k) rLife(n=10000, object=ibge2022M_ltb, x = k)) %>% 
                   map_dbl(.f = mean)

vida_media_fem <- map(0:getOmega(ibge2022F_ltb), \(k) rLife(n=10000, object=ibge2022F_ltb, x = k)) %>% 
                  map_dbl(.f = mean)

Os comandos acima simulam 10.000 realizações aleatórias do tempo de vida futuro para cada uma das idades da tábua e obtém a média dos valores simulados para cada idade.

A vida média, ou expectativa de vida, é uma função biométrica fornecida pelas tábuas de mortalidade, de forma que é possível comparar os valores obtidos acima, via simulação, com os valores teóricos constantes da tábua, o que faremos a seguir por intermédio de um gráfico de dispersão:

# Masculino

ibge2022M %>% 
  mutate(vm_masc = vida_media_masc) %>% 
  ggplot(aes(x = ex, y = vm_masc)) +
  geom_point(shape = 1, size = 2) +
  geom_abline(slope = 1, intercept = 0, color="blue") +
  labs(title = "Comparação das expectativas de vida\n(Masculino)",
       x = "ex (Tábua)",
       y = "ex (Simulação)") +
  theme_bw()

# Feminino

ibge2022F %>% 
  mutate(vm_fem = vida_media_fem) %>% 
  ggplot(aes(x = ex, y = vm_fem)) +
  geom_point(shape = 1, size = 2) +
  geom_abline(slope = 1, intercept = 0, color="blue") +
   labs(title = "Comparação das expectativas de vida\n(Feminino)",
       x = "ex (Tábua)",
       y = "ex (Simulação)") +
  theme_bw()

Como pode ser visto, os dois valores são praticamente iguais para todas as idades, o que nos mostra que usando simulação podemos conseguir excelentes aproximações dos valores teóricos.

Simulação da Vida Provável

O mesmo que fizemos acima pode ser feito para calcular a mediana dos valores simulados, o que é feito a seguir:

# mediana do tempo de vida futuro
vida_mediana_masc <- map(0:getOmega(ibge2022M_ltb), \(k) rLife(n=10000, object=ibge2022M_ltb, x = k)) %>% 
                     map_dbl(.f = median)
  
vida_mediana_fem <- map(0:getOmega(ibge2022F_ltb), \(k) rLife(n=10000, object=ibge2022F_ltb, x = k)) %>% 
                    map_dbl(.f = median)

Note

Devemos observar que, para fins didáticos, refizemos novamente as 10.000 simulações e calculamos a mediana dos valores simulados para cada idade, mas poderíamos ter reutilizado os valores anteriormente simulados para calcular a vida média.

Vamos comparar a vida provável com a vida média para os homens.

ibge2022M %>% 
  mutate(vida_provavel = vida_mediana_masc) %>%
  select(x, ex, vida_provavel) %>% 
  pivot_longer(cols = c(ex, vida_provavel), names_to = "tipo_vida", values_to = "vlr") %>% 
  ggplot(aes(x = x, y = vlr, color=tipo_vida)) +
  geom_line(linewidth = 1.5) +
  labs(x = "Idades (x)", y="Vida Média | Vida provável")

Agora o mesmo para as mulheres:

ibge2022F %>% 
  mutate(vida_provavel = vida_mediana_fem) %>%
  select(x, ex, vida_provavel) %>% 
  pivot_longer(cols = c(ex, vida_provavel), names_to = "tipo_vida", values_to = "vlr") %>% 
  ggplot(aes(x = x, y = vlr, color=tipo_vida)) +
  geom_line(linewidth = 1.5) +
  labs(x = "Idades (x)", y = "Vida Média | Vida provável")

Poderíamos calcular qualquer outra medida resumo para os valores simulados do tempo de vida futuro para cada idade.

Uma vantagem da simulação é que podemos estipular intervalos que definam uma certa probabilidade de conter as quantidades calculadas (vida média, vida provável, etc), em vez de fornecer apenas uma estimativa pontual da quantidade desejada.

Uma aplicação prática

dados <- read_excel("BaseDadosVanzzilotta.xlsx") %>% 
             mutate(IDADE = round(time_length(difftime(ymd("2019-12-31"), as.Date(`DT NASC`)), "years")))  

head(dados) %>% kable()

IDENTIF	DT NASC	DT ADMISSÃO	DT ADESÃO	SEXO	SALARIO	RESGATE	IDADE
1	1962-10-15	1980-02-15	1983-03-15	M	18844.5	289620.5	57
2	1955-02-17	1981-05-15	1981-06-06	M	20604.0	453101.0	65
3	1963-05-17	1983-11-18	1983-11-18	M	12665.0	177650.0	57
4	1954-03-12	1986-02-22	1986-02-22	M	9401.0	178211.0	66
5	1968-02-05	1990-03-18	1990-03-18	F	5652.5	44591.0	52
6	1970-04-05	1992-03-28	1992-03-28	M	19764.0	331423.5	50

Importada a base de dados, fizemos o cálculo das idades dos indivíduos na data de 31/12/2019.

Essa base de dados possui apenas 100 indivíduos e está sendo utilizada para fins didáticos. Na prática, os arquivos contendo a base cadastral para os cálculos atuariais conterão muito mais do que 100 observações o que implicará uma maior quantidade de tempo para a realização das simulações.

Para realizar as simulações é necessário separar as idades dos indivíduos do sexo masculino das idades dos indivíduos do sexo feminino:

idades_masc <- dados %>% filter(SEXO == "M") %>% pull(IDADE)
idades_fem  <- dados %>% filter(SEXO == "F") %>% pull(IDADE)

Separadas as idades das pessoas por sexo, podemos realizar as simulações. Vamos simular 10.000 realizações do tempo de vida futuro para cada pessoa de ambos os sexos.

sim_Tx_masc <- map(idades_masc, \(k) rLife(n=10000, object=ibge2022M_ltb, x = k)) 
sim_Tx_fem  <- map(idades_fem,  \(k) rLife(n=10000, object=ibge2022F_ltb, x = k))

Os objetos sim_Tx_masc e sim_Tx_fem são listas em que cada componente possui os resultados das 10.000 simulações realizadas do tempo de vida futuro para cada idade.

Quantidade de pessoas em cada lista:

c("Feminino"   = length(sim_Tx_fem),
  "Masculino"  = length(sim_Tx_masc))

 Feminino Masculino 
       32        68

Obtidos os valores simulados podemos, como já visto acima, calcular quaisquer medidas resumo a partir das 10.000 simulações realizadas para cada indivíduo da base de dados.

O percentil de ordem 95 dos tempos de vida futuro é o valor para o qual 95% dos valores simulados serão menores ou iguais a ele e 5% superiores, indicando que a probabilidade de se obter um tempo de vida futuro superior é de apenas 5%.

Para exemplificar, vamos calcular o percentil de ordem 95, a vida média e a vida mediana e reunir estes valores como colunas em um data frame contendo também as idades e o sexo dos indivíduos.

medidas_resumo <- bind_rows(
  
  tibble(x    = idades_fem,
         sexo = "F",
         ex   = map_dbl(sim_Tx_fem, mean),
         p50  = map_dbl(sim_Tx_fem, median),
         p95  = map_dbl(sim_Tx_fem, \(x) quantile(x, probs = 0.95))), 
         
  tibble(x   = idades_masc,
        sexo = "M",
        ex   =  map_dbl(sim_Tx_masc, mean),
        p50  =  map_dbl(sim_Tx_masc, median),
        p95  =  map_dbl(sim_Tx_masc, \(x) quantile(x, probs = 0.95)))
  
)

head(medidas_resumo) %>% kable()

x	sexo	ex	p50	p95
52	F	30.3091	31.5	46.5
48	F	33.9491	35.5	50.5
43	F	38.5783	40.5	55.5
52	F	30.4776	32.0	46.5
58	F	25.1753	26.5	40.5
45	F	36.5411	38.5	53.5

Obtivemos as medidas resumo desejadas a partir das listas sim_Tx_fem e sim_Tx_masc que contém os 10.000 valores simulados do tempo de vida futuro para cada indivíduo na base de dados.

Eficiência computacional

Uma outra questão levantada no artigo, e talvez até mais importante que a inadequação da função rLife(), diz respeito à velocidade com que as simulações devem ser feitas para um conjunto de dados de tamanho compatível com o que vamos encontrar na prática.

[Os RPPS de pequeno porte são os que mais se beneficiariam da utilização de simulações para o cálculo das quantidades atuariais.

A esse respeito recomendamos a leitura da tese de doutoramento da Prof^a Cristine Correa intitulada “Tamanho populacional e aleatoriedade de eventos demográficos na solvência em RPPS municipais capitalizados” disponível para download em https://repositorio.ufmg.br/handle/1843/AMSA-9TNH49]{.aside}

O porte de um RPPS municipal é dado pela quantidade de beneficiários (ativos, aposentados e pensionistas) que o mesmo possui, sendo classificados em RPPS de grande porte, médio porte e pequeno porte.

Os RPPS estaduais e distrital integram a categoria de porte especial e independe da quantidade de beneficiários.

Com base nos dados do ISP de 2024 elaboramos o quadro a seguir que nos indica os quantitativos mínimo e máximo de beneficiários em RPPS de cada um dos quatro portes acima mencionados.

Porte	Qtd. de RPPS	Qtd. Mín. Benef.	Qtd. Máx. Benef.
Especial	27	18.561	808.143
Grande	105	6.051	233.302
Médio	943	641	5.966
Pequeno	1.047	1	640

Com base nos dados acima, parece razoável admitir que seria factível a realização de simulações sem maiores preocupações com eficiência computacional para todos os RPPS de pequeno e médio portes, que representam mais de 90% dos RPPS.

Vamos a seguir calcular o tempo de execução de algumas simulações e ver como esses tempos podem ser melhorados.

Note

A velocidade de processamento depende fortemente da configuração do computador utilizado, de forma que é possível que o leitor possa obter valores bem distintos dos que são aqui apresentados.

Vamos iniciar calculando o tempo para o exemplo que acabamos de mostrar acima. Para isso vamos usar o pacote {tictoc}.

tic()
sim_Tx_masc <- map(idades_masc, \(k) rLife(n=10000, object=ibge2022M_ltb, x = k)) 
sim_Tx_fem  <- map(idades_fem,  \(k) rLife(n=10000, object=ibge2022F_ltb, x = k)) 
t1 <- toc()

2.23 sec elapsed

As 10.000 simulações para cada indivíduo dessa base de dados bem pequena demorou apenas 2.23 segundos.E se aumentarmos o tamanho dessa base de dados de 100 para 10.000 pessoas? Quanto tempo demoraria para realizar as 10.000 simulações para cada pessoa?

Vamos criar uma base de dados “sintética” a partir da base de dados utilizada acima por intermédio de reamostragem:

n <- 10000 
dadosextendido <- tibble(idade = sample(x = dados$IDADE, size = n, replace = TRUE),
                         sexo  = sample(x = dados$SEXO,  size = n, replace = TRUE))
  
nrow(dadosextendido)

[1] 10000

Como se pode ver, a nova base de dados chamada dadosextendido possui 10.000 pessoas.

Separação das idades das pessoas do sexo masculino e feminino:

x_masc <- dadosextendido %>% filter(sexo == "M") %>% pull(idade)
x_fem  <- dadosextendido %>% filter(sexo == "F") %>% pull(idade)

c("Masculino" = length(x_masc),
  "Feminino"  = length(x_fem))

Masculino  Feminino 
     6854      3146

Dessas 10.000 pessoas, 6854 são do sexo masculino e 3146 são do sexo feminino.

Vamos à medição do tempo necessário para realizar a simulação do tempo de vida futuro para essa quantidade de pessoas:

tic()
sim_Tx_masc <- map(x_masc, \(k) rLife(n=10000, object=ibge2022M_ltb, x = k)) 
sim_Tx_fem  <- map(x_fem,  \(k) rLife(n=10000, object=ibge2022F_ltb, x = k)) 
t2 <- toc()

253.25 sec elapsed

Realizamos as simulações para alguns tamanhos de bases de dados (\(n\)) e os resultados são apresentados na tabela a seguir:

n	tempo (aprox.)
500	13 seg.
1000	27 seg.
5000	125 seg. (2 min.)
10000	253.25 seg. (4 min.)

Apenas com esses recursos, sem buscar qualquer melhorias de eficiência, poderíamos realizar simulações para RPPS de portes pequeno e médio e ainda alguns de grande porte conforme quadro acima.

Paralelização

Embora os tempos de execução das simulaçõs não sejam impraticáveis, uma pergunta importante é: será que teria alguma forma de melhorar esses tempos?

A resposta é sim, usando o recurso da paralelização que diversos pacotes no R implementam, e também
usando pacotes que sabidamente são mais rápidos na manipulação de dados como o {data.table}.

Nessa sessão vamos tratar apenas da paralelização, que diz respeito ao aspecto da eficiência computacional relacionado à velocidade de processamento, mas existe também o aspecto relacionado ao consumo de memória.

Abaixo vamos deixar alguns links para materiais que possibilitarão um aprofundamento na temática da computação em paralelo que será necessária caso se deseje de fato implementar a realização de calculos atuariais via simulação na prática dos RPPS:

https://nceas.github.io/oss-lessons/parallel-computing-in-r/parallel-computing-in-r.html https://dept.stat.lsa.umich.edu/~jerrick/courses/stat701/notes/parallel.html https://bookdown.org/rdpeng/rprogdatascience/parallel-computation.html https://www.appsilon.com/post/r-doparallel https://github.com/fastverse/fastverse https://cran.r-project.org/web/views/HighPerformanceComputing.html

Para mais informações sobre esse pacote recomendamos consultar a documentação disponível em https://furrr.futureverse.org/

Embora não seja nosso objetivo discutir em maior profundidade ganhos de eficiência computacional relacionados à velocidade de processamento, gostaríamos apenas de comentar que o R disponibiliza vários recursos que possibilitam essa abordagem, sendo uma delas o pacote {furrr} que vamos usar para ilustrar, de forma introdutória, a implementação de paralelização no R.

Para utilizar esse pacote, precisamos identificar e definir a quantidade de “núcleos” que desejamos utilizar, o que fazemos a seguir:

No windows é possível consultar a quantidade de núcleos disponíveis a partir do menu “Desempenho” do “Gerenciador de Tarefas” que é acessado com “Ctrl+Shift+Esc.

plan(multisession, workers = 10)

Agora vamos repetir a simulação feita para a base de dados sintética de 10.000 indivíduos e medir o tempo necessário, realizando a simulação de forma paralizada usando a função future_map() do pacote {furrr}

tic()
simTx_masc <- future_map(x_masc, \(k) rLife(n=10000, object=ibge2022M_ltb, x = k),
                          .options = furrr_options(seed = 123))

simTx_fem  <- future_map(x_fem,  \(k) rLife(n=10000, object=ibge2022F_ltb, x = k),
                          .options = furrr_options(seed = 123)) 
t3 <- toc()

138.92 sec elapsed

Como pode ser visto, o tempo passou de 253.25 seg. para 138.92 seg. Um ganho de aproximadamente 45% no tempo de processamento.

Considerações finais

Ao longo do tutorial fizemos algumas escolhas que, embora corretas e didaticamente úteis, talvez não sejam as melhores do ponto de vista prático e da eficiência computacional.

Por exemplo, usamos sempre a quantidade de 10.000 simulações por indivíduo, mas será que para efeitos práticos apenas 1.000 simulações por indivíduo não nos forneceriam estimativas suficientemente precisas?

Outro ponto é que optamos por usar o vetor de idades para realizar as simulações exatamene da forma como estavam na base de dados. Assim, por exemplo, se na base de dados existir 10 pessoas de 35 anos do sexo masculino, realizamos 10.000 simulações para cada uma delas, totalizando 100.000 simulações. Mas será que não seria suficiente realizar apenas 10.000 simulações (ou 1.000) e usar os resultados da simulação para todos os indivíduos de 35 anos do sexo masculino existentes na base de dados?

Essas são questões que merecem uma maior reflexão e que podem fazer muita diferença na prática.

Até a próxima elucubração atuarial!

Introdução

Dados utilizados

A função rLife() e o tempo de vida futuro (\(T(X)\))

Uma aplicação prática

Eficiência computacional

Paralelização

Considerações finais

A função `rLife()` e o tempo de vida futuro (\(T(X)\))