
<style>
body {
  width: 600px;
}
</style>

```{r echo=FALSE, message=FALSE, warnings=FALSE}
opts_chunk$set(echo = FALSE, messages = FALSE, 
               fig.height = 6, fig.width = 6)
```

# Disasters worldwide from 1900-2008

```{r results='hide', message=FALSE, warning=FALSE, error=FALSE}
library("ggplot2")
library("reshape")
library("rworldmap")
library("colorspace")
library("maptools")
library("plotrix")

data("wrld_simpl")

source("mapCountryData.R")
source("helpers.R")
```

This document reproduces the example "Disasters worldwide from 1900-2008"
presented in the paper "[Probabilistic archetypal analysis](http://arxiv.org/abs/1312.7604)" by Seth and Eugster (2014). The [Rmd file](Disaster-report.Rmd) lists the 
corresponding [R](http://www.r-project.org) code.


```{r results='asis'}
cat("* Read tourist data from file `emdata.csv`\n")
emdata <- read.table("data_Disaster/emdata.csv", stringsAsFactors = FALSE, sep = ";", header = TRUE)
emdata_mat <- as.matrix(emdata[, -1])

cat("* Load precomputed PAA solution from file `Disaster_paa.mat`\n")
paa <- load_paa("data_Disaster/Disaster_paa.mat", emdata_mat)
```


#### Cost plot:

```{r fig.height=4}
op <- par(mar = c(4, 4, 1, 1))
plot(paa$obj, type = "l", xlab = "Iterations", ylab = "Cost")
par(op)
```


#### Disaster profiles:

```{r results='asis'}
colnames(paa$Z) <- sprintf("A%s", 1:7)

print(xtable(paa$Z), type = "html")
```


#### Barplot of profiles (% of maximum value):

```{r}
aa_perc <- t(perc(paa$Z, emdata_mat))

colnames(aa_perc) <- colnames(emdata_mat)
rownames(aa_perc) <- sprintf("A%s", seq(ncol(paa$Z)))
colnames(aa_perc) <- c("CD", "DR", "EQ", "EP", "ET", "FL",
                       "IA", "II", "MD", "MW", "MA", "ST",
                       "TA", "VO", "WF")

aa_perc[aa_perc > 1] <- 1.0
aa_perc <- aa_perc * 100

aa_perc <- melt(aa_perc)
aa_perc <- rename(aa_perc, c(X2 = "Disaster", X1 = "Archetype"))

ot <- theme_set(theme_bw(base_size = 9))

ggplot(aa_perc, aes(Disaster, value)) +
    geom_bar(stat = "identity") +
    facet_grid(Archetype ~ ., scales = "fixed") +
    theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust = 0.5),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()) +
    scale_y_continuous(name = "Archetype profiles (%)")
```


#### World maps colored based on factor $H$ for each archetype:

```{r }
plot_map <- function(H, i, borderLwd = 0.01, borderCol = "white") {
  dat <- data.frame(H3 = H[i, ], Country = emdata$Country, stringsAsFactors = FALSE)
  map <- joinCountryData2Map(dat, joinCode = "NAME", nameJoinColumn = "Country")

  cols <- rev(heat_hcl(51, c = c(80, 30), l = c(30, 90), power = c(1/5, 2)))
  cols[1] <- "lightgray"

  op <- par(mai = c(0,0,0,0), xaxs = "i", yaxs = "i")
  mp <- mapCountryData(map, nameColumnToPlot = "H3",
                       catMethod = "pretty", numCat = 50,
                       colourPalette = cols, addLegend = FALSE,
                       mapTitle = "", borderCol = borderCol,
                       borderLwd = borderLwd,
                       missingCountryCol = "lightgray")
  text(-152, 40, sprintf("A%s", i))

  invisible(mp)
}
```

```{r message=FALSE, warning=FALSE, results='hide', fig.height=3, fig.width=6}
plot_map(paa$H, 1)
plot_map(paa$H, 2)
plot_map(paa$H, 3)
plot_map(paa$H, 4)
plot_map(paa$H, 5)
plot_map(paa$H, 6)
plot_map(paa$H, 7)
```


#### Simplex visualization showing the model deviance:

```{r}
paa_perc <- paa
paa_perc$Z <- perc(paa$Z, emdata_mat)
```

```{r}
dev <- poission_deviance(paa_perc, emdata_mat)
val <- (dev / max(dev))
val <- (1 - (exp(-val * 10))) / (1 - exp(-10))
ci <- classIntervals(val, style = "pretty", n = 10)
point_cols <- findColours(ci, sequential_hcl(11))
```

```{r}
op <- par(mar = c(1.5, 1.5, 1.5, 1.5), xpd = TRUE)
set.seed(1239)
simplexplot_paa(paa_perc, show_circle = FALSE, 
             show_points = TRUE, points_col = point_cols,
             projection = tspsimplex_projection, show_edges = TRUE, labels_cex = 1)

t <- names(attr(point_cols, "table"))

legend("bottomleft", fill = attr(point_cols, "palette"), 
       border = attr(point_cols, "palette"),
       legend = t, bg = "white",
       horiz = FALSE, cex = 0.8, box.lwd = NA)

par(op)
```


#### Simplex visualization showing insect infestations variable:

```{r}
## Variable of interest:
i <- 8

## Color coding:
point_cols <- ifelse(emdata_mat[, i] > 0, "black", "darkgray")
dir_cols <- ifelse(emdata_mat[, i] > 0, "black", NA)
iicntry <- emdata_mat[, i]/max(emdata_mat[, i])
point_cexs <- ifelse(iicntry == 0, 0.4, iicntry)


## Simplex plot:
op <- par(mar = c(1.5, 1.5, 1.5, 1.5), xpd = TRUE)
set.seed(1239)
s <- simplexplot_paa(paa_perc, show_circle = FALSE, points_col = point_cols, show_direction=TRUE,
                  show_points = TRUE, directions_col = dir_cols, points_cex = point_cexs * 2,
                  projection = tspsimplex_projection, show_edges = TRUE, labels_cex = 1)

iso2 <- match(emdata$Country, wrld_simpl$NAME)
iso2 <- wrld_simpl$ISO2[iso2]
 
w <- which(iicntry > 0.5)
w <- w[order(iicntry[w], decreasing = TRUE)]

pc <- function(i, xshift = 0, yshift = 0) {
  boxed.labels(s$proj_h[i, 1, drop = FALSE] + xshift, s$proj_h[i, 2, drop = FALSE] + yshift,
               labels = iso2[i], bg = "white", col = "black", border = FALSE, cex = 1)
}

pc(w[1], yshift = 0.8)
pc(w[2], yshift = 0.8)
pc(w[3], yshift = -0.8)
pc(w[4], yshift = -0.8)
pc(w[5], yshift = 0.8)
pc(w[6], xshift = -0.8)
pc(w[7], yshift = -0.8)
pc(w[8], yshift = -0.8)
pc(w[9], yshift = 0.8)

p <- seq(0, 1, length.out = 5)

legend("bottomleft", pch = 19, col = "black", 
       title = "Insect\ninfestation",
       pt.cex = p * 2,
       legend = paste(p*100, "%", sep = ""), bg = "white",
       horiz = FALSE, cex = 1, box.lwd = NA)

par(op)
```

