options(width = 100) knitr::opts_chunk$set(tidy = FALSE, size = "small")
In order to assess the accuracy of life-table summary measures over varying levels of data aggregation I distinguish between 3 levels of age discreteness:
The better the discretization of the summary measure the closer its level 1 and 2 values are to the value calculated directly from level 0 data.
library(pash) library(tidyverse) radix = 10000 # create pash objects lt_lvl1 <- Inputlx(x = australia_1y$x, lx = australia_1y$lx, nax = australia_1y$nax, nx = australia_1y$nx, last_open = TRUE) lt_lvl2 <- Inputlx(x = australia_10y$x, lx = australia_10y$lx, nax = australia_10y$nax, nx = australia_10y$nx, last_open = TRUE)
# not a shape measure, just as a test, should be exactly equal e0 <- c( lvl0 = mean(australia_surv_times), lvl1 = GetPace(lt_lvl1, type = "e0"), lvl2 = GetPace(lt_lvl1, type = "e0") ) e0
Formulation of individual level Gini coefficient: Damgaard, C., & Weiner, J. (2000). Describing inequality in plant size or fecundity. Ecology, 81(4), 1139-1142.
X = data.frame(i = 1:length(australia_surv_times), x = sort(australia_surv_times)) G_lvl0 = sum((2*X$i-nrow(X)-1)*X$x)/(nrow(X)^2*mean(australia_surv_times)) gini <- c( lvl0 = G_lvl0, lvl1 = GetShape(lt_lvl1, type = "gini", harmonized = FALSE), lvl2 = GetShape(lt_lvl2, type = "gini", harmonized = FALSE) ) gini
cv <- c( lvl0 = sd(australia_surv_times) / mean(australia_surv_times), lvl1 = GetShape(lt_lvl1, type = "cv", harmonized = FALSE), lvl2 = GetShape(lt_lvl2, type = "cv", harmonized = FALSE) ) cv
eDaggerx = vector(mode = "double", length = length(australia_surv_times)) for (i in seq_along(australia_surv_times)) { if (i < length(australia_surv_times)) { eDaggerx[i] = mean(australia_surv_times[seq(i+1, length(australia_surv_times))]) - australia_surv_times[i] } } entropy <- c( lvl0 = sum(eDaggerx)/radix/mean(australia_surv_times), lvl1 = GetShape(lt_lvl1, type = "entropy", harmonized = FALSE), lvl2 = GetShape(lt_lvl2, type = "entropy", harmonized = FALSE) ) entropy
psmad <- c( lvl0 = sum(australia_surv_times > mean(australia_surv_times)) / radix, lvl1 = GetShape(lt_lvl1, type = "psmad", harmonized = FALSE), lvl2 = GetShape(lt_lvl2, type = "psmad", harmonized = FALSE) ) psmad
# remaining life-expectancy at mean age of death (expressed as share of total life-expectancy) ler <- c( lvl0 = (mean(australia_surv_times[australia_surv_times > mean(australia_surv_times)])-mean(australia_surv_times))/mean(australia_surv_times), lvl1 = GetShape(lt_lvl1, type = "ler", harmonized = FALSE), lvl2 = GetShape(lt_lvl2, type = "ler", harmonized = FALSE) ) ler
data.frame(cv, gini, entropy, ler, psmad) %>% mutate(lvl = rownames(.)) %>% gather(key = measure, value = value, -lvl) %>% group_by(measure) %>% do({ rbind( data.frame(lvl = 1, error = (.$value[2]-.$value[1])/.$value[1]), data.frame(lvl = 2, error = (.$value[3]-.$value[1])/.$value[1]) ) }) %>% ggplot() + geom_col(aes(x = interaction(lvl, measure), y = error, fill = measure)) + scale_x_discrete("", labels = c("CV-1Y.", "CV-10Y", "Entropy-1Y", "Entropy-10Y", "Gini-1Y", "Gini-10Y", "LER-1Y", "LER-10Y", "PSMAD-1Y", "PSMAD-10Y")) + scale_y_continuous(labels = scales::percent_format(), limits = c(-0.05, 0.05), breaks = seq(-0.05, 0.05, 0.01)) + theme_minimal() + theme(panel.grid.minor = element_blank(), panel.grid.major.y = element_blank()) + guides(fill = FALSE) + coord_flip()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.