## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "##",
  warning = FALSE,
  message = FALSE,
  fig.width = 7,
  fig.height = 5
)
library(HistData)
library(ggplot2)

## ----arbuthnot-duplicate------------------------------------------------------
data(Arbuthnot)
Arbuthnot$Christenings <- Arbuthnot$Males + Arbuthnot$Females
Arbuthnot[Arbuthnot$Year %in% c(1674, 1704),
          c("Year", "Males", "Females", "Christenings", "Ratio")]

## ----arbuthnot-plot-----------------------------------------------------------
copy_level <- Arbuthnot$Christenings[Arbuthnot$Year == 1704]

ggplot(Arbuthnot, aes(x = Year, y = Christenings)) +
  geom_line(color = "grey50") +
  geom_point(color = "grey50") +
  geom_hline(yintercept = copy_level, linetype = "dotted", color = "black") +
  annotate(
    "segment", x = 1677, xend = 1701, y = copy_level, yend = copy_level,
    color = "darkred", linewidth = 1.3,
    arrow = arrow(length = unit(0.3, "cm"), type = "closed")
  ) +
  geom_point(
    data = subset(Arbuthnot, Year %in% c(1674, 1704)),
    color = "firebrick", size = 4
  ) +
  labs(
    title = "Arbuthnot's annual christenings in London, 1629-1710",
    subtitle = "1674 and 1704 (red) are identical -- likely a copying error",
    y = "Christenings (Males + Females)"
  ) +
  theme_minimal()

## ----minard-ggplot, fig.height = 3.5------------------------------------------
data(Minard.cities); data(Minard.troops)

ggplot(Minard.troops, aes(x = long, y = lat)) +
  geom_path(aes(linewidth = survivors, colour = direction, group = group),
            lineend = "round") +
  scale_linewidth(range = c(0.5, 15), guide = "none") +
  scale_colour_manual(values = c(A = "grey50", R = "firebrick"), guide = "none") +
  geom_point(data = Minard.cities, size = 1.2) +
  geom_text(data = Minard.cities, aes(label = city), vjust = 1.5, hjust = 0.5, size = 3) +
  labs(title = "Napoleon's Russian campaign of 1812", x = NULL, y = NULL) +
  theme_minimal() +
  theme(axis.text = element_blank(), panel.grid = element_blank())

## ----perozzo-setup------------------------------------------------------------
data(Perozzo)
Pmat <- xtabs(Survivors ~ Year + Age, data = Perozzo)
years <- as.numeric(rownames(Pmat))
ages  <- as.numeric(colnames(Pmat))

## ----perozzo-base, fig.height = 5---------------------------------------------
ages_rev <- -rev(ages)
Pmat_rev <- Pmat[, rev(seq_along(ages))]
persp(years, ages_rev, Pmat_rev,
      xlab = "Year", ylab = "Age", zlab = "Survivors",
      theta = 0, phi = 25, expand = 0.6,
      col = adjustcolor("lightblue", alpha.f = 0.5), shade = 0.5)

## ----perozzo-enhanced, fig.height = 6-----------------------------------------
colnames(Pmat_rev) <- as.character(ages_rev)  # match column names to ages_rev's sign flip

# re-angled to bring the Age = 0 "wall" of large young-survivor counts forward
pmat <- persp(years, ages_rev, Pmat_rev,
      xlab = "Year", ylab = "Age", zlab = "Survivors",
      theta = -30, phi = 20, expand = 0.6, r = 4,
      col = adjustcolor("lightblue", alpha.f = 0.4), shade = 0.5)

# red: highlight both cross-section families -- Survivors-by-Age at each Year,
# and Survivors-by-Year at each Age -- every 25 years/years-of-age
hi_years <- years[years %% 25 == 0]
for (yr in hi_years) {
  z <- Pmat_rev[as.character(yr), ]
  lines(trans3d(rep(yr, length(ages_rev)), ages_rev, z, pmat), col = "red", lwd = 2)
}
hi_ages <- ages[ages %% 25 == 0]
for (ag in hi_ages) {
  z <- Pmat_rev[, as.character(-ag)]
  lines(trans3d(years, rep(-ag, length(years)), z, pmat), col = "red", lwd = 2)
}

# blue: birth-cohort diagonals -- Year = birth year + Age -- every 25 years,
# starting from the Age = 0 edge (cohorts born during the observation window)
# and the Year = 1750 edge (cohorts already alive when it starts)
cohort_starts <- rbind(
  data.frame(year0 = years[years %% 25 == 0], age0 = 0),
  data.frame(year0 = 1750, age0 = ages[ages %% 25 == 0 & ages > 0])
)
for (i in seq_len(nrow(cohort_starts))) {
  b  <- cohort_starts$year0[i]
  a0 <- cohort_starts$age0[i]
  k  <- 0:((100 - a0) / 5)
  Y  <- b + 5 * k
  A  <- a0 + 5 * k
  keep <- Y >= min(years) & Y <= max(years) & A <= 100
  Y <- Y[keep]; A <- A[keep]
  if (length(Y) < 2) next
  z <- mapply(function(y, a) Pmat[as.character(y), as.character(a)], Y, A)
  lines(trans3d(Y, -A, z, pmat), col = "blue", lwd = 2)
}

## ----perozzo-ggplot-contour---------------------------------------------------
ggplot(Perozzo, aes(x = Year, y = Age, z = Survivors)) +
  geom_raster(aes(fill = Survivors)) +
  geom_contour(color = "white", alpha = 0.6) +
  scale_fill_distiller(palette = "Blues", direction = 1) +
  scale_y_reverse() +
  labs(title = "Perozzo's survivorship grid, as a heatmap") +
  theme_minimal()

## ----perozzo-isometric--------------------------------------------------------
shear <- 0.6   # horizontal points per year of age; tune to taste
lift  <- 400   # vertical Survivors-units per year of age

Perozzo_iso <- transform(Perozzo,
  x_iso = Year + shear * Age,
  y_iso = Survivors + lift * Age
)

ggplot(Perozzo_iso, aes(x = x_iso, y = y_iso, group = Age, color = Age)) +
  geom_line() +
  scale_color_distiller(palette = "Blues", direction = 1) +
  labs(
    title = "Perozzo's stereogram as a 2D isometric shear",
    x = "Year (sheared by Age)", y = "Survivors (lifted by Age)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank())

## ----perozzo-plotly-----------------------------------------------------------
library(plotly)
plot_ly(x = ~years, y = ~ages, z = ~Pmat, type = "surface") |>
  layout(scene = list(
    xaxis = list(title = "Year"),
    yaxis = list(title = "Age"),
    zaxis = list(title = "Survivors")
  ))

## ----snow-density-------------------------------------------------------------
data(Snow.deaths); data(Snow.pumps); data(Snow.streets)

ggplot(Snow.deaths, aes(x = x, y = y)) +
  geom_line(data = Snow.streets, aes(x = x, y = y, group = street),
            color = "black", linewidth = 0.5) +
  geom_point(alpha = 0.6) +
  geom_density_2d_filled(bins = 4, show.legend = FALSE) +
  geom_density_2d(bins = 4, color = "black", alpha = 0.6) +
  scale_fill_manual(values = c("#FFFFFF00", "#F0808030", "#FF000030", "#8B000080")) +
  geom_label(data = Snow.pumps, aes(x = x, y = y, label = label),
             size = 3, color = "grey20", fill = "lightblue", alpha = 0.9) +
  theme_void() +
  labs(
    title = "John Snow's Cholera Map",
    subtitle = "With contour lines showing densities of death",
    caption = "Data: HistData::Snow.deaths / Snow.pumps / Snow.streets"
  )