Nightingale.R

Uploaded by:mdevlin

              require(plotrix)
require(zoo)
data('Nightingale', package = 'HistData')

# For some graphs, it is more convenient to reshape death rates to long format
#  keep only Date and death rates

require(reshape2)
Night <- Nightingale[, c(1, 8:10)]
melted <- melt(Night, "Date")
names(melted) <- c("Date", "Cause", "Deaths")
melted$Cause <- sub("\\.rate", "", melted$Cause)
melted$Regime <- ordered(rep(c(rep('Before', 12), rep('After', 12)), 3)
                         , levels = c('Before', 'After'))
Night <- melted

# subsets, to facilitate separate plotting
Night1 <- subset(Night, Date < as.Date("1855-04-01"))
Night2 <- subset(Night, Date >= as.Date("1855-04-01"))

## Not run: 
require(ggplot2)
# Before plot
cxc1 <- ggplot(Night1, aes(x = factor(Date), y=Deaths, fill = Cause)) +
    geom_bar(width = 1, stat="identity", color="black") +   # a stacked bar chart first
    scale_y_sqrt()   # set scale so area ~ Deaths  

# A coxcomb plot = bar chart + polar coordinates
cxc1 + coord_polar(start=3*pi/2) + 
  ggtitle("Causes of Mortality in the Army in the East") + 
  xlab("")

# After plot
cxc2 <- ggplot(Night2, aes(x = factor(Date), y=Deaths, fill = Cause)) +
  geom_bar(width = 1, stat="identity", color="black") +
  scale_y_sqrt()

cxc2 + coord_polar(start=3*pi/2) +
  ggtitle("Causes of Mortality in the Army in the East") + 
  xlab("")

# do both together, with faceting
cxc <- ggplot(Night, aes(x = factor(Date), y=Deaths, fill = Cause)) +
  geom_bar(width = 1, stat="identity", color="black") + 
  scale_y_sqrt() +
  facet_grid(. ~ Regime, scales="free", labeller=label_both)

cxc + coord_polar(start=3*pi/2) +
  ggtitle("Causes of Mortality in the Army in the East") + 
  xlab("")

## End(Not run)

## What if she had made a set of line graphs?

# these plots are best viewed with width ~ 2 * height 
colors <- c("blue", "red", "black")

with(Nightingale, {
  plot(Date, Disease.rate, type="n", cex.lab=1.25, 
       ylab="Annual Death Rate", xlab="Date", xaxt="n",
       main="Causes of Mortality of the British Army in the East");
  
  # background, to separate before, after
  rect(as.Date("1854/4/1"), -10, as.Date("1855/3/1"), 
       1.02*max(Disease.rate), col=gray(.90), border="transparent");
  text( as.Date("1854/4/1"), .98*max(Disease.rate), "Before Sanitary\nCommission", pos=4);
  text( as.Date("1855/4/1"), .98*max(Disease.rate), "After Sanitary\nCommission", pos=4);
  
  # plot the data
  points(Date, Disease.rate, type="b", col=colors[1], lwd=3);
  points(Date, Wounds.rate, type="b", col=colors[2], lwd=2);
  points(Date, Other.rate, type="b", col=colors[3], lwd=2)
}
)

# add custom Date axis and legend

axis.Date(1, at=seq(as.Date("1854/4/1"), as.Date("1856/3/1"), "3 months"), format="%b %Y")
legend(x = 'right', c("Preventable disease", "Wounds and injuries", "Other"),
       col=colors, fill=colors, title="Cause", cex=0.75, inset = .01)

# Alternatively, show each cause of death as percent of total
Nightingale <- within(Nightingale, {
  Total <- Disease + Wounds + Other
  Disease.pct <- 100*Disease/Total
  Wounds.pct <- 100*Wounds/Total
  Other.pct <- 100*Other/Total
})

colors <- c("blue", "red", "black")

with(Nightingale, {
  plot(Date, Disease.pct, type="n",  ylim=c(0,100), cex.lab=1.25,
       ylab="Percent deaths", xlab="Date", xaxt="n",
       main="Percentage of Deaths by Cause");
  # background, to separate before, after
  rect(as.Date("1854/4/1"), -10, as.Date("1855/3/1"), 
       1.02*max(Disease.rate), col=gray(.90), border="transparent");
  text( as.Date("1854/4/1"), .98*max(Disease.pct), "Before Sanitary\nCommission", pos=4);
  text( as.Date("1855/4/1"), .98*max(Disease.pct), "After Sanitary\nCommission", pos=4);
  # plot the data
  points(Date, Disease.pct, type="b", col=colors[1], lwd=3);
  points(Date, Wounds.pct, type="b", col=colors[2], lwd=2);
  points(Date, Other.pct, type="b", col=colors[3], lwd=2)
}
)

# add custom Date axis and legend
axis.Date(1, at=seq(as.Date("1854/4/1"), as.Date("1856/3/1"), "3 months"), format="%b %Y")
legend(x = 'center', c("Preventable disease", "Wounds and injuries", "Other"),
       col=colors, fill=colors, title="Cause", cex=.75)

###

fn <- read.table("../fn.data",header=TRUE,sep="\t")

dates    <- as.Date(as.yearmon(fn$Date, "%b %Y"))
MthlyZym <- (fn$DeathZym/fn$ArmySize) * 1000 
MthlyWnd <- (fn$DeathWnds/fn$ArmySize) * 1000
MthlyAll <- (fn$DeathAll/fn$ArmySize) * 1000

op <- par(mfrow=c(1,2), pty="s")
par(cex.axis=0.5) # shrink outer labels for mfrow mode
par(cex.lab=0.5)
permudates12bef <- c(as.character(fn$Date[7:12]), as.character(fn$Date[1:6]))
fnstart12 <- 2*pi*6/12 # start at pi radians
radial.pie(sqrt(MthlyZym[1:12]), clockwise=TRUE, show.grid.labels=TRUE, 
           start=fnstart12, labels=permudates12bef, label.prop=1.1, sector.colors=rep("cornflowerblue",12))
radialtext("X", center=c(-10.5,0),col="red")
title(main="BEFORE: Apr 1854 - Mar 1855") # main arg doesn't work

permudates12aft <- c(as.character(fn$Date[19:24]), as.character(fn$Date[13:18]))
fnstart12 <- 2*pi*6/12 # start at pi radians
radial.pie(sqrt(MthlyZym[13:24]), clockwise=TRUE, show.grid.labels=TRUE, start=fnstart12, 
           labels=permudates12aft, label.prop=1.1, sector.colors=rep("blue3",12))
radialtext("X", center=c(-5.5,0),col="red")
title(main="AFTER: Apr 1855 - Mar 1856")
par(op)