Source of Data

Data Input

Graunt’s Life Table

rm(list = ls())
graunt <- data.frame(x = c(0, seq(6, 76, by = 10)), lx.17th = c(100, 64, 40, 25, 16, 10, 6, 3, 1))

More data

US 1993 life table for the same age group

us.93 <- data.frame(x = graunt$x, lx.93 = c(100, 99, 99, 98, 97, 95, 92, 84, 70))

Data Extraction

There are many ways to extract part of us.93 data frame.

us.93["lx.93"]

##   lx.93
## 1   100
## 2    99
## 3    99
## 4    98
## 5    97
## 6    95
## 7    92
## 8    84
## 9    70

us.93["lx.93"][[1]]

## [1] 100  99  99  98  97  95  92  84  70

us.93["lx.93"]$lx.93

## [1] 100  99  99  98  97  95  92  84  70

us.93["lx.93"]$lx

## [1] 100  99  99  98  97  95  92  84  70

us.93[2]

##   lx.93
## 1   100
## 2    99
## 3    99
## 4    98
## 5    97
## 6    95
## 7    92
## 8    84
## 9    70

us.93[2][[1]]

## [1] 100  99  99  98  97  95  92  84  70

us.93[2]$lx.93

## [1] 100  99  99  98  97  95  92  84  70

us.93[, "lx.93"]

## [1] 100  99  99  98  97  95  92  84  70

us.93[, 2]

## [1] 100  99  99  98  97  95  92  84  70

us.93$lx.93

## [1] 100  99  99  98  97  95  92  84  70

us.93$lx

## [1] 100  99  99  98  97  95  92  84  70

Into one single data frame

Combine two data frames into one single data frame

(graunt.us <- data.frame(graunt, lx.93 = us.93$lx))

##    x lx.17th lx.93
## 1  0     100   100
## 2  6      64    99
## 3 16      40    99
## 4 26      25    98
## 5 36      16    97
## 6 46      10    95
## 7 56       6    92
## 8 66       3    84
## 9 76       1    70

Life Expectancy

The basic principle is that the area under the survival function is the life expectancy.

\(X \ge 0\), \(X \sim F(x)\) => \(X \equiv F^{-1}(U), U \sim U(0,1)\), therefore,

\(E(X) = E\{F^{-1}(U)\} = \int_{0}^{1} F^{-1}(u)du = \int_0^{\infty} 1-F(x) dx = \int_{0}^{\infty} S(x) dx\)

Step by step approach to draw survival function plot

Basic plot with points and lines, compare the following threes methods

# library(extrafont)
par(mfrow = c(2, 2))
plot(graunt$x, graunt$lx)
plot(lx.17th ~ x, data = graunt)
plot(graunt)
plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")

Denote the ages and observed survival rates on the axes

plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")
axis(side = 1, at = graunt$x, labels = graunt$x)
axis(side = 2, at = graunt$lx.17th, labels = graunt$lx.17th)

Denote the age 0 and 76 by dotted lines

plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")
axis(side = 1, at=graunt$x, labels=graunt$x)
axis(side = 2, at = graunt$lx.17th, labels = graunt$lx.17th)
abline(v = c(0, 76), lty = 2)

Setting up coordinates for `polygon()` (Clockwise)

(graunt.x <- c(graunt$x, 0))

##  [1]  0  6 16 26 36 46 56 66 76  0

(graunt.y <- c(graunt$lx.17th, 0))

##  [1] 100  64  40  25  16  10   6   3   1   0

graunt.poly <- data.frame(x = graunt.x, y = graunt.y)

Shading

plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")
axis(side = 1, at = graunt$x, labels = graunt$x)
axis(side = 2, at = graunt$lx.17th, labels = graunt$lx.17th)
abline(v = c(0, 76), lty = 4)
polygon(graunt.poly, density = 15, angle = 135)
# polygon(graunt.x, graunt.y, density = 15, angle = 135)
points(graunt, pch = 21, col = "black", bg = "white")

Grids

plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")
axis(side = 1, at = graunt$x, labels = graunt$x)
axis(side = 2, at = graunt$lx.17th, labels = graunt$lx.17th)
abline(v = c(0, 76), lty = 2)
polygon(graunt.poly, density = 15)
# polygon(graunt.x, graunt.y, density = 15)
abline(v = graunt$x, lty = 2)
points(graunt, pch = 21, col = "black", bg = "white")

Title, x-axis label, and y-axis label

plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")
axis(side = 1, at = graunt$x, labels = graunt$x)
axis(side = 2, at = graunt$lx.17th, labels = graunt$lx.17th)
abline(v = c(0, 76), lty = 2)
polygon(graunt.poly, density = 15)
# polygon(graunt.x, graunt.y, density = 15)
abline(v = graunt$x, lty = 2)
points(graunt, pch = 21, col = "black", bg = "white")
main.title <- "Graunt's Survival Function"
x.lab <- "Age (years)"
y.lab <- "Survival Rates (%)"
title(main = main.title, xlab = x.lab, ylab = y.lab)

Area under the curve

The area under the curve can be approximated by the sum of the areas of trapezoids, therefore the area is \(\sum_{i=1}^{n-1} (x_{i+1}-x_i)\times\frac{1}{2}(y_i + y_{i+1})\).

diff(), head(), and tail() can be used to write a function to compute the area easily.

area.R <- function(x, y) {
  sum(diff(x) * (head(y, -1) + tail(y, -1))/2)
  }
area.R(graunt$x, graunt$lx.17th)/100

## [1] 18.17

Comparison with US 1993 life table

The shaded area between the survival functions of Graunt’s and US 1993 represents the difference of life expectancies.

Draw Graunt’s first with axes, lower and upper limits

plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")
axis(side = 1, at = graunt$x, labels = graunt$x)
axis(side = 2, at = graunt$lx, labels = graunt$lx.17th)
abline(v = c(0, 76), lty = 2)

Add US 1993 survival function

plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")
axis(side = 1, at = graunt$x, labels = graunt$x)
axis(side = 2, at = graunt$lx, labels = graunt$lx.17th)
abline(v = c(0, 76), lty = 2)
lines(us.93, type = "b")

Actually, US 1993 life table is truncated at the age 76. Specify that point.

plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")
axis(side = 1, at = graunt$x, labels = graunt$x)
axis(side = 2, at = graunt$lx, labels = graunt$lx.17th)
abline(v = c(0, 76), lty = 2)
lines(us.93, type = "b")
abline(h = 70, lty = 2)

Using las = 1 to specify 70%.

plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")
axis(side = 1, at = graunt$x, labels = graunt$x)
axis(side = 2, at = graunt$lx, labels = graunt$lx.17th)
abline(v = c(0, 76), lty = 2)
lines(us.93, type = "b")
abline(h = 70, lty = 2)
axis(side = 2, at = 70, labels = 70, las = 1)

Setting coordinates for `polygon()`

us.graunt.x <- c(us.93$x, rev(graunt$x))
us.graunt.y <- c(us.93$lx.93, rev(graunt$lx.17th))
us.graunt <- data.frame(x = us.graunt.x, y = us.graunt.y)

Shading

plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")
axis(side = 1, at = graunt$x, labels = graunt$x)
axis(side = 2, at = graunt$lx, labels = graunt$lx.17th)
abline(v = c(0, 76), lty = 2)
lines(us.93, type = "b")
abline(h = 70, lty = 2)
axis(side = 2, at = 70, labels = 70, las = 1)
polygon(us.graunt, density = 15, col = "red", border = NA)
# polygon(us.graunt.x, us.graunt.y, density = 15, col = "red", border = NA)
points(us.graunt, pch = 21, col = "black", bg = "white")

Grids

plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")
axis(side = 1, at = graunt$x, labels = graunt$x)
axis(side = 2, at = graunt$lx, labels = graunt$lx.17th)
abline(v = c(0, 76), lty = 2)
lines(us.93, type = "b")
abline(h = 70, lty = 2)
axis(side = 2, at = 70, labels = 70, las = 1)
polygon(us.graunt, density = 15, col = "red", border = NA)
abline(v = graunt$x, lty = 2)
points(us.graunt, pch = 21, col = "black", bg = "white")

Title, x-axis and y-axis labels

plot(graunt, ann = FALSE, xaxt = "n", yaxt = "n", type = "b")
axis(side = 1, at = graunt$x, labels = graunt$x)
axis(side = 2, at = graunt$lx, labels = graunt$lx.17th)
abline(v = c(0, 76), lty = 2)
lines(us.93, type = "b")
abline(h = 70, lty = 2)
axis(side = 2, at = 70, labels = 70, las = 1)
polygon(us.graunt, density = 15, col = "red", border = NA)
abline(v = graunt$x, lty = 2)
points(us.graunt, pch = 21, col = "black", bg = "white")
main.title.g.us <- "Survival Function of Graunt and US 1993"
title(main = main.title.g.us, xlab = x.lab, ylab = y.lab)

dev.copy(device = png, file = "../pics/graunt_us93.png")

## quartz_off_screen 
##                 3

Life expectancy

The area under the US 1993 survival function is

area.R(us.93$x, us.93$lx.93)/100

## [1] 70.92

The area of shaded region is

area.R(us.93$x, us.93$lx.93)/100 - area.R(graunt$x, graunt$lx.17th)/100

## [1] 52.75

Comparison with Halley’s life table

Halley’s life table

age <- 0:84
lx <- c(1238, 1000, 855, 798, 760, 732, 710, 692, 680, 670, 661, 653, 646, 640, 634, 628, 622, 616, 610, 604, 598, 592, 586, 579, 573, 567, 560, 553, 546, 539, 531, 523, 515, 507, 499, 490, 481, 472, 463, 454, 445, 436, 427, 417, 407, 397, 387, 377, 367, 357, 346, 335, 324, 313, 302, 292, 282, 272, 262, 252, 242, 232, 222, 212, 202, 192, 182, 172, 162, 152, 142, 131, 120, 109, 98, 88, 78, 68, 58, 50, 41, 34, 28, 23, 20)
length(lx)

## [1] 85

halley <- data.frame(age, lx)
halley$px <- round(halley$lx/1238*100, digits = 1)
head(halley)

##   age   lx    px
## 1   0 1238 100.0
## 2   1 1000  80.8
## 3   2  855  69.1
## 4   3  798  64.5
## 5   4  760  61.4
## 6   5  732  59.1

tail(halley)

##    age lx  px
## 80  79 50 4.0
## 81  80 41 3.3
## 82  81 34 2.7
## 83  82 28 2.3
## 84  83 23 1.9
## 85  84 20 1.6

halley.lx <- halley[-3]
halley <- halley[-2]
head(halley)

##   age    px
## 1   0 100.0
## 2   1  80.8
## 3   2  69.1
## 4   3  64.5
## 5   4  61.4
## 6   5  59.1

tail(halley)

##    age  px
## 80  79 4.0
## 81  80 3.3
## 82  81 2.7
## 83  82 2.3
## 84  83 1.9
## 85  84 1.6

R base graphics

To make the comparison easy, plot the points at the same age group of Graunt’s, 0, 6, 16, 26, 36, 46, 56, 66, 76. Step by step approach

Halley’s survival function first

plot(halley, ann = FALSE, xaxt = "n", yaxt = "n", type = "l")

Denote the age at 0, 76, and 84 by vertical dotted lines

plot(halley, ann = FALSE, xaxt = "n", yaxt = "n", type = "l")
abline(v = c(0, 76, 84), lty = 2)

Mark the points at 0, 6, 16, 26, 36, 46, 56, 66, 76 on Halley’s survival function.

age.graunt <- age %in% graunt$x
plot(px ~ age, data = halley, ann = FALSE, xaxt = "n", yaxt = "n", type = "l")
abline(v = c(0, 76, 84), lty = 2)
points(px[age.graunt] ~ age[age.graunt], data = halley, pch = 21, col = "black", bg = "white")

Using `subset()`

halley.graunt <- subset(halley, age.graunt)
plot(halley, ann = FALSE, xaxt = "n", yaxt = "n", type = "l")
abline(v = c(0, 76, 84), lty = 2)
points(halley.graunt, pch = 21, col = "black", bg = "white")

Add Graunt’s survival function

plot(halley, ann = FALSE, xaxt = "n", yaxt = "n", type = "l")
abline(v = c(0, 76, 84), lty = 2)
points(halley.graunt, pch = 21, col = "black", bg = "white")
lines(graunt, type = "b", pch = 21, col = "black", bg = "white")

x-axis label and y-axis label with las = 1

plot(halley, ann = FALSE, xaxt = "n", yaxt = "n", type = "l")
abline(v = c(0, 76, 84), lty = 2)
points(halley.graunt, pch = 21, col = "black", bg = "white")
lines(graunt, type = "b", pch = 21, col = "black", bg = "white")
axis(side = 1, at = c(graunt$x, 84), labels = c(graunt$x, 84))
axis(side = 2, at = graunt$lx.17th, labels = graunt$lx.17th, las = 1)

Specify the developers at proper coordinates with text()

plot(halley, ann = FALSE, xaxt = "n", yaxt = "n", type = "l")
abline(v = c(0, 76, 84), lty = 2)
points(halley.graunt, pch = 21, col = "black", bg = "white")
lines(graunt, type = "b", pch = 21, col = "black", bg = "white")
axis(side = 1, at = c(graunt$x, 84), labels = c(graunt$x, 84))
axis(side = 2, at = graunt$lx.17th, labels = graunt$lx.17th, las = 1)
text(x = c(16, 36), y = c(20, 50), label = c("Graunt", "Halley"))

Main title, x-axis label, and y-axis label

plot(halley, ann = FALSE, xaxt = "n", yaxt = "n", type = "l")
abline(v = c(0, 76, 84), lty = 2)
points(halley.graunt, pch = 21, col = "black", bg = "white")
lines(graunt, type = "b", pch = 21, col = "black", bg = "white")
axis(side = 1, at = c(graunt$x, 84), labels = c(graunt$x, 84))
axis(side = 2, at = graunt$lx.17th, labels = graunt$lx.17th, las = 1)
text(x = c(16, 36), y = c(20, 50), label = c("Graunt", "Halley"))
main.title.2 <- "Survival Function of Graunt and Halley"
title(main = main.title.2, xlab = x.lab, ylab = y.lab)

Polygon

Setting the coordinates for polygon(). The intersection is found at x = 10.8, y = 52.8 with locator(1) and couple of trial and errors.

Upper region

poly.1.x <- c(graunt$x[1:2], 10.8, halley$age[11:1])
poly.1.y <- c(graunt$lx.17th[1:2], 52.8, halley$px[11:1])
poly.upper <- data.frame(x = poly.1.x, y = poly.1.y)

Lower region

poly.2.x <- c(10.8, halley$age[12:85], graunt$x[9:3])
poly.2.y <- c(52.8, halley$px[12:85], graunt$lx.17th[9:3])
poly.lower <- data.frame(x = poly.2.x, y = poly.2.y)

Shading upper region first

plot(halley, ann = FALSE, xaxt = "n", yaxt = "n", type = "l")
abline(v = c(0, 76, 84), lty = 2)
points(halley.graunt, pch = 21, col = "black", bg = "white")
lines(graunt, type = "b", pch = 21, col = "black", bg = "white")
axis(side = 1, at = c(graunt$x, 84), labels = c(graunt$x, 84))
axis(side = 2, at = graunt$lx.17th, labels = graunt$lx.17th, las = 1)
text(x = c(16, 36), y = c(20, 50), label = c("Graunt", "Halley"))
title(main = main.title.2, xlab = x.lab, ylab = y.lab)
polygon(poly.upper, angle = 45, density = 15, col = "blue")

Shading lower region next

plot(halley, ann = FALSE, xaxt = "n", yaxt = "n", type = "l")
abline(v = c(0, 76, 84), lty = 2)
points(halley.graunt, pch = 21, col = "black", bg = "white")
lines(graunt, type = "b", pch = 21, col = "black", bg = "white")
axis(side = 1, at = c(graunt$x, 84), labels = c(graunt$x, 84))
axis(side = 2, at = graunt$lx.17th, labels = graunt$lx.17th, las = 1)
text(x = c(16, 36), y = c(20, 50), label = c("Graunt", "Halley"))
title(main = main.title.2, xlab = x.lab, ylab = y.lab)
polygon(poly.upper, angle = 45, density = 15, col = "blue")
polygon(poly.lower, angle = 45, density = 15, col = "red")

Fill the points. Extra points at the 84.

plot(halley, ann = FALSE, xaxt = "n", yaxt = "n", type = "l")
abline(v = c(0, 76, 84), lty = 2)
points(halley.graunt, pch = 21, col = "black", bg = "white")
lines(graunt, type = "b", pch = 21, col = "black", bg = "white")
axis(side = 1, at = c(graunt$x, 84), labels = c(graunt$x, 84))
axis(side = 2, at = graunt$lx.17th, labels = graunt$lx.17th, las = 1)
text(x = c(16, 36), y = c(20, 50), label = c("Graunt", "Halley"))
title(main = main.title.2, xlab = x.lab, ylab = y.lab)
polygon(poly.upper, angle = 45, density = 15, col = "blue")
polygon(poly.lower, angle = 45, density = 15, col = "red")
points(graunt, pch = 21, col = "black", bg = "white")
points(halley.graunt, pch = 21, col = "black", bg = "white")
points(x = 84, y = halley$px[85], pch = 21, col = "black", bg = "white")

dev.copy(device = png, file = "../pics/graunt_halley.png")

## quartz_off_screen 
##                 4

Life expectancy

Compute the difference of life expectancies

(life.exp.halley <- area.R(halley$age, halley$px)/100)

## [1] 27.872

(life.exp.graunt <- area.R(graunt$x, graunt$lx.17th)/100)

## [1] 18.17

ggplot

library(ggplot2)

Data Reshape

Attach reshape2 package to change wide format to long format

library(reshape2)

How melt() works

graunt.us.melt <- melt(graunt.us, id.vars = "x", measure.vars = c("lx.17th", "lx.93"), value.name = "lx", variable.name = "times")
graunt.us.melt

##     x   times  lx
## 1   0 lx.17th 100
## 2   6 lx.17th  64
## 3  16 lx.17th  40
## 4  26 lx.17th  25
## 5  36 lx.17th  16
## 6  46 lx.17th  10
## 7  56 lx.17th   6
## 8  66 lx.17th   3
## 9  76 lx.17th   1
## 10  0   lx.93 100
## 11  6   lx.93  99
## 12 16   lx.93  99
## 13 26   lx.93  98
## 14 36   lx.93  97
## 15 46   lx.93  95
## 16 56   lx.93  92
## 17 66   lx.93  84
## 18 76   lx.93  70

str(graunt.us.melt)

## 'data.frame':    18 obs. of  3 variables:
##  $ x    : num  0 6 16 26 36 46 56 66 76 0 ...
##  $ times: Factor w/ 2 levels "lx.17th","lx.93": 1 1 1 1 1 1 1 1 1 2 ...
##  $ lx   : num  100 64 40 25 16 10 6 3 1 100 ...

Change factor levels of times

levels(graunt.us.melt$times) <- c("17th", "1993")
str(graunt.us.melt)

## 'data.frame':    18 obs. of  3 variables:
##  $ x    : num  0 6 16 26 36 46 56 66 76 0 ...
##  $ times: Factor w/ 2 levels "17th","1993": 1 1 1 1 1 1 1 1 1 2 ...
##  $ lx   : num  100 64 40 25 16 10 6 3 1 100 ...

Plot

Points and Lines

Step by step approach to understand the grammar of ggplot

We set ggplot() to accept varying data.frame() and aes()in geom_polygon

(g1 <- ggplot() + 
  geom_line(data = graunt.us.melt, aes(x = x, y = lx, colour = times)))

(g2 <- g1 + 
  geom_point(data = graunt.us.melt, aes(x = x, y = lx, colour = times), shape = 21, fill = "white"))

(g3 <- g2 + 
  theme_bw())

Polygon

Reuse us.graunt which contains x = us.graunt.x and y = us.graunt.y for polygon(). Note that we start with g3, and also note how to remove default legends.

(p3 <- g3 + 
  geom_polygon(data = us.graunt, aes(x = x, y = y), alpha = 0.3, fill = "red"))

(p4 <- p3 + 
  guides(colour = "none"))

Change default annotations

Points and Lines

Change the x-axis and y-axis labels

(g4 <- g3 + 
   xlab(x.lab) + 
   ylab(y.lab))

Add main title

(g4 <- g3 + 
   xlab(x.lab) + 
   ylab(y.lab) + 
   ggtitle(main.title.g.us))

Change legend title

(g4 <- g3 + 
   xlab(x.lab) + 
   ylab(y.lab) + 
   ggtitle(main.title.g.us) +
   labs(colour = "Era"))

Change legends.

(g4 <- g3 + 
   xlab(x.lab) + 
   ylab(y.lab) +
   ggtitle(main.title.g.us) +
   labs(colour = "Era") +
   scale_colour_discrete(labels = c("Graunt Era", "US 1993")))

Place legends inside the plot

(g5 <- g4 + 
   theme(legend.position = c(0.8, 0.5)))

Change x-axis and y-axis tick marks

(g6 <- g5 + 
   scale_x_continuous(breaks = graunt$x) + scale_y_continuous(breaks = graunt$lx.17th))

ggsave("../pics/graunt_us_plot.png", g6)

## Saving 6 x 6 in image

Polygon

Add information to the plot drawn with polygon()

Start with p4

p4

Main title, x-axis and y-axis labels

(p5 <- p4 + 
   xlab(x.lab) + 
   ylab(y.lab) +
   ggtitle(main.title.g.us))

"Graunt Era", "US 1993", "Difference of Life Expectancies" at proper positions

(p6 <- p5 + 
   annotate("text", x = c(20, 40, 70), y = c(20, 60, 90), label = c("Graunt Era", "Difference of\nLife Expectancies", "US 1993"), family = "Helvetica"))

x-axis and y-axis tick marks

(p7 <- p6 + 
   scale_x_continuous(breaks = graunt$x) + scale_y_continuous(breaks = graunt$lx.17th))

ggsave("../pics/graunt_us_poly.png", p7)

## Saving 6 x 6 in image

Graunt and Halley

Data Reshaping

Since the observed ages are different, we need final structure of the data frame to be melted. So, create copies of graunt and halley and extract parts of what we need and give feasible names.

graunt.2 <- graunt
halley.2 <- halley
names(graunt.2) <- c("x", "Graunt")
names(halley.2) <- c("x", "Halley")
graunt.halley.melt <- melt(list(graunt.2, halley.2), id.vars = "x", value.name = "lx", variable.name = "Who")
str(graunt.halley.melt)

## 'data.frame':    94 obs. of  4 variables:
##  $ x  : num  0 6 16 26 36 46 56 66 76 0 ...
##  $ Who: Factor w/ 2 levels "Graunt","Halley": 1 1 1 1 1 1 1 1 1 2 ...
##  $ lx : num  100 64 40 25 16 10 6 3 1 100 ...
##  $ L1 : int  1 1 1 1 1 1 1 1 1 2 ...

graunt.halley.melt <- graunt.halley.melt[-4]
(graunt.halley.melt.g <- subset(graunt.halley.melt, graunt.halley.melt$x %in% graunt$x))

##     x    Who    lx
## 1   0 Graunt 100.0
## 2   6 Graunt  64.0
## 3  16 Graunt  40.0
## 4  26 Graunt  25.0
## 5  36 Graunt  16.0
## 6  46 Graunt  10.0
## 7  56 Graunt   6.0
## 8  66 Graunt   3.0
## 9  76 Graunt   1.0
## 10  0 Halley 100.0
## 16  6 Halley  57.4
## 26 16 Halley  50.2
## 36 26 Halley  45.2
## 46 36 Halley  38.9
## 56 46 Halley  31.3
## 66 56 Halley  22.8
## 76 66 Halley  14.7
## 86 76 Halley   6.3

Survival Function, Step by Step

(gh1 <- ggplot() + 
  geom_line(data = graunt.halley.melt, aes(x = x, y = lx, colour = Who)))

(gh2 <- gh1 + 
  geom_point(data = graunt.halley.melt.g, aes(x = x, y = lx, colour = Who), shape = 21, fill = "white"))

(gh3 <- gh2 + 
  theme_bw() + 
  xlab(x.lab) + 
  ylab(y.lab) + 
  ggtitle(main.title.2))

(gh4 <- gh3 + 
  theme(legend.position = c(0.8, 0.8)) +
  annotate("text", x = c(16, 36), y = c(20, 50), label = c("Graunt", "Halley")) +
  scale_x_continuous(breaks = c(graunt$x, 84)) + 
  scale_y_continuous(breaks = c(graunt$lx.17th, halley$px[halley$age == 6])))

ggsave("../pics/graunt_halley_ggplot.png", gh4)

## Saving 6 x 6 in image

Polygon

Reuse poly.upper data frame and poly.lower data frame.

(ghp4 <- gh4 + 
  geom_polygon(data = poly.upper, aes(x = x, y = y), alpha = 0.3, fill = "red"))

(ghp5 <- ghp4 + 
  geom_polygon(data = poly.lower, aes(x = x, y = y), alpha = 0.3, fill = "green"))

(ghp5 <- ghp5 +
  geom_point(data = data.frame(x = 84, y = halley$px[85]), aes(x = x, y = y),  colour = 3, shape = 21, fill = "white"))

ggsave("../pics/graunt_halley_poly_ggplot.png", ghp5)

## Saving 6 x 6 in image

Graunt, Halley, and US93

Data Reshape

us93.2 <- us.93
names(us93.2) <- c("x", "US93")
ghu.melt <- melt(list(graunt.2, halley.2, us93.2), id.vars = "x", value.name = "lx", variable.name = "Who")
ghu.melt.g <- ghu.melt[ghu.melt$x %in% graunt$x, ]
main.title.3 <- "Survival Function Plots"

Survival Function Plots with ggplot

(ghu1 <- ggplot() + 
  geom_line(data = ghu.melt, aes(x = x, y = lx, colour = Who)))

(ghu2 <- ghu1 + 
  geom_point(data = ghu.melt.g, aes(x = x, y = lx, colour = Who), shape = 21, fill = "white"))

(ghu3 <- ghu2 + 
  theme_bw() + 
  xlab(x.lab) + 
  ylab(y.lab) + 
  ggtitle(main.title.3))

(ghu4 <- ghu3 + 
  theme(legend.position = c(0.2, 0.2)) +
  annotate("text", x = c(36, 36, 70), y = c(25, 50, 90), label = c("Graunt", "Halley", "US93")) +
  scale_x_continuous(breaks = c(graunt$x, 84)) + 
  scale_y_continuous(breaks = c(graunt$lx.17th, halley$px[halley$age == 6])))

ggsave("../pics/graunt_halley_us_ggplot.png", ghu4)

## Saving 6 x 6 in image

Polygon

Coordinates

In order to make the graphs truncated at the age 76, restrict the age of Halley up to 76.

poly.lower.76 <- subset(poly.lower, poly.lower$x <= 76)
poly.3.x <- c(us93.2$x, halley.2$x[85:12], 10.8, graunt.2$x[2:1])
poly.3.y <- c(us93.2$US93, halley.2$Halley[85:12], 52.8, graunt.2$Graunt[2:1])
poly.us <- data.frame(x = poly.3.x, y = poly.3.y)
poly.us.76 <- subset(poly.us, poly.us$x <= 76)

(ghup4 <- ghu4 + 
  geom_polygon(data = poly.upper, aes(x = x, y = y), alpha = 0.3, fill = "red"))

(ghup5 <- ghup4 + 
  geom_polygon(data = poly.lower.76, aes(x = x, y = y), alpha = 0.3, fill = "green"))

(ghup6 <- ghup5 +
  geom_polygon(data = poly.us.76, aes(x = x, y = y), alpha = 0.3, fill = "blue"))

(ghup7 <- ghup6 +
  geom_point(data = data.frame(x = 84, y = halley$px[85]), aes(x = x, y = y),  colour = 3, shape = 21, fill = "white"))

ggsave("../pics/graunt_halley_us_poly_ggplot.png", ghup7)

## Saving 6 x 6 in image

dump() and source()

Check out how to save and retrieve. Use source() and load() for retrieval.

dump("area.R", file = "area.R")
save.image("graunt_halley_160329v2.rda")

Graunt’s and Halley’s Life Table with ggplot

coop711

2016-04-06

Source of Data

Data Input

More data

Data Extraction

Into one single data frame

Life Expectancy

Step by step approach to draw survival function plot

Setting up coordinates for `polygon()` (Clockwise)

Area under the curve

Comparison with US 1993 life table

Setting coordinates for `polygon()`

Life expectancy

Comparison with Halley’s life table

Halley’s life table

R base graphics

Using `subset()`

Polygon

Life expectancy

ggplot

Data Reshape

Plot

Points and Lines

Polygon

Change default annotations

Points and Lines

Polygon

Graunt and Halley

Data Reshaping

Survival Function, Step by Step

Polygon

Graunt, Halley, and US93

Data Reshape

Survival Function Plots with ggplot

Polygon

Coordinates

dump() and source()

Graunt’s and Halley’s Life Table with ggplot

coop711

2016-04-06

Source of Data

Data Input

More data

Data Extraction

Into one single data frame

Life Expectancy

Step by step approach to draw survival function plot

Setting up coordinates for polygon() (Clockwise)

Area under the curve

Comparison with US 1993 life table

Setting coordinates for polygon()

Life expectancy

Comparison with Halley’s life table

Halley’s life table

R base graphics

Using subset()

Polygon

Life expectancy

ggplot

Data Reshape

Plot

Points and Lines

Polygon

Change default annotations

Points and Lines

Polygon

Graunt and Halley

Data Reshaping

Survival Function, Step by Step

Polygon

Graunt, Halley, and US93

Data Reshape

Survival Function Plots with ggplot

Polygon

Coordinates

dump() and source()

Setting up coordinates for `polygon()` (Clockwise)

Setting coordinates for `polygon()`

Using `subset()`