6_visualization.Rmd

```{r, echo=FALSE}
cat(paste("(C) (cc by-sa) Wouter van Atteveldt, file generated", format(Sys.Date(), format="%B %d %Y")))
```

> Note on the data used in this howto: 
> This data can be downloaded from http://piketty.pse.ens.fr/files/capital21c/en/xls/, 
> but the excel format is a bit difficult to parse at it is meant to be human readable, with multiple header rows etc. 
> For that reason, I've extracted csv files for some interesting tables that I've uploaded to 
> http://vanatteveldt.com/uploads/rcourse/data

Plotting data in R
===

In this hands-on we continue with the `capital` variable created in the [transforming data howto](4_transforming.md).
You can also download this variable from the course pages:

```{r}
download.file("http://vanatteveldt.com/uploads/rcourse/data/capital.rdata", destfile="capital.rdata")
load("capital.rdata")
head(capital)
```

We also make a 'wide' vesion of this data frame using melt, and we turn the first column into the row names so all columns contain actual data:

```{r}
library(reshape)
wide = cast(capital, Year ~ Country, value="Total")
head(wide)
```

Plotting multiple lines with a loop
----

In many cases, just calling ~plot~ on an R object gives a good quick visualizations.
Sometimes, however, we want to have more control over the exact placing and look of a graph, e.g. to prepare a graph for publication.
In such cases, it is usually best to start with an empty plot, and then put in all the plot elements such as the lines but also the axes, legend, etc.
We first determine the needed x and y limits to fit all the data, and then plot an empty plot with those limits:


```{r}
ylim = c(min(wide[-1], na.rm=T), max(wide[-1], na.rm=T))
xlim = c(min(wide$Year), max(wide$Year))
plot(0,  type="n", xlim=xlim, ylim=ylim, frame.plot=F, xlab='Year', ylab='Capital', axes=F,
     main='Capital accumulation as proportion of national income')
```

This gives us an empty 'canvas' to start drawing plot elements in. Let's add lines, axes, a legend, and a dashed horizontal line at y=0.  For plotting multiple lines, it is often best to use a for loop over the data frame columns, and draw from a predefined color set such as the ~rainbow~ colors. 


```{r}
plot(0,  type="n", xlim=xlim, ylim=ylim, frame.plot=F, xlab='Year', ylab='Capital', axes=F,
     main='Capital accumulation as proportion of national income')
axis(2) # normal vertical axis
axis(1, at=seq(xlim[1], xlim[2], 5), las=2) # specify vertical ticks every 5 years
colors = rainbow(ncol(wide)) 
for (i in 2:ncol(wide)) {
  lines(x=wide$Year, y=wide[[i]], col=colors[i])  
}
legend("topleft", legend=colnames(wide)[-1], col=colors[-1], lty=1, cex=.65, ncol=2)
# regression line for US case
abline(lm(U.S. ~ Year, wide), col=colors[2], lty=2)
```

*Note*: In RStudio, the exact positioning of plot elements is determined by the size of the plot window.
If you are preparing a plot for publication, you normally want to produce an external plot.
It is usually best to create the plot as a file directly and open it in a second window in a photo viewer, 
so the positioning is correct for the file, rather than for the RStudio window. 
You create a plot file by opening it using the `png` function, running the plot commands, and then closing the file with `dev.off()`.

```
png("fig1.png", width=1600, height=1200)
... plot commands ...
dev.off()
```

Bar plots
----

The visualizations so far were mainly line graphs.
We could also create a bar plot of e.g. the average accumulation of capital per country:

```{r}
total.capital = aggregate(capital$Total, capital["Country"], FUN=mean)
barplot(total.capital$x, names.arg=total.capital$Country, las=2)
```

We can also put e.g. the public and private wealth side by side.
Annoyingly, the barplot function needs a matrix with the countries in the columns,
so we use `as.matrix` to convert the aggregated data to a matrix form, and then use `t()` to transpose it;

```{r}
total.capital = aggregate(capital[c("Public","Private")], capital["Country"], FUN=mean)
m = t(as.matrix(total.capital[-1]))
m
barplot(m, names.arg=total.capital$Country, las=2, beside=T, col=rainbow(2), ylim=c(0,6))
legend("top", c("Public", "Private"), fill=rainbow(2), ncol=2)
```

Combined bar/line plots
---

Sometimes it is more intuitive to combine bars and lines into a single plot. 
For example, we could combine the capital data with the inequality data used earlier and plot the latter as a line.
First, we download the income data again, select the years from 1970, and melt the data into the same form as `capital`:
```{r}
library(plyr)
download.file("http://vanatteveldt.com/wp-content/uploads/rcourse/data/income_toppercentile.csv", 
              destfile="income_toppercentile.csv")
income = read.csv("income_toppercentile.csv")
income = rename(income, c("US" = "U.S."))
income = income[income$Year >= 1970, ]
income.long = melt(income, id.vars="Year", na.rm=T)
colnames(income.long) = c("Year", "Country", "income")
head(income.long)
```

Now, we can compute the total income inequality and combine that with the capital accumulation:

```{r}
total.income = aggregate(income.long["income"], income.long["Country"], FUN=mean)
total = merge(total.capital, total.income)
total
```

Note that by default, merge only keeps rows for with both data frames have values.
By specifying `all=T` you can keep all rows, getting NA for missing values. 
Note that we need to specify a secondary axis for the inequality to account for the different scale. 

```{r}
par(mar=c(5,5,5,5))
x = barplot(total$Private, names.arg=total$Country, las=2, ylab="Capital accumulation")
lines(x=x, y=total$income*20)
axis(side=4,at=0:5, labels=0:5 * 20)
mtext("Top percentile share of income", 4, line=3)
```

Scatter plots
---

We can also create a scatter plot of inequality versus capital accumulation:

```{r}
plot(x=total$Private, y=total$income, xlim=c(3,6), frame.plot=F)
text(x=total$Private, y=total$income, labels=total$Country, pos=4)
```

Note that calling plot on a data frame with all interval columns automatically creates pairwise scatter plots:

```{r}
plot(total[-1])
```

This can also be done using the `pairs` function, which moreover can add a 'panel' line to give a visual indication of the seeming linear or curved relation between two variables:

```{r}
pairs(total[-1], panel=panel.smooth)
```

Finally, let's create a scatter plot for all yearly values rather than the totals per country. 

```{r}
d = merge(capital, income.long)
au = d[d$Country == "Australia",]
ca = d[d$Country == "Canada",]
us = d[d$Country == "U.S.",]
fr = d[d$Country == "France",]
plot(0, ylim=c(2, 6), xlim=c(0.05, 0.25), type="n")
points(y=au$Private, x=au$income, col="red")
points(y=ca$Private, x=ca$income, col="blue")
points(y=us$Private, x=us$income, col="darkgreen")
points(y=fr$Private, x=fr$income, col="purple")
```

It seems that for australia and canada there is a fairly linear relation between capital and income inequality. 
Using the `abline` command based on a linear model fit, we can plot the regression lines for all points as well:

```{r}
plot(0, ylim=c(2, 6), xlim=c(0.05, 0.25), type="n", xlab="Income inequality", ylab="Capital Accumulation", frame.plot=F)
points(y=au$Private, x=au$income, col="red")
abline(lm(au$Private ~ au$income), col="red")
points(y=ca$Private, x=ca$income, col="blue")
abline(lm(ca$Private ~ ca$income), col="blue")
points(y=us$Private, x=us$income, col="darkgreen")
abline(lm(us$Private ~ us$income), col="darkgreen")
points(y=fr$Private, x=fr$income, col="purple")
abline(lm(fr$Private ~ fr$income), col="purple")
legend(.17, 3, legend=c("Australia", "Canada", "US", "France"), ncol=2, fill=c("red", "blue", "darkgreen", "purple"))
```

(Of course, this code can be made more generic by using a for loop rather than copying the commands for each country,
which is left as an exercise for the reader.