Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
190 lines (151 sloc) 7.6 KB
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 https://github.com/vanatteveldt/learningr/tree/master/data. If you're accessing this tutorial from the githup project, these files should be in your 'data' sub folder automatically.

Plotting data in R

In this hands-on we continue with the capital variable created in the transforming data howto. You can also download this variable from the course pages:

load("data/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:

library(reshape2)
wide = dcast(capital, Year ~ Country, value.var="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:

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.

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:

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;

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:

library(plyr)
income = read.csv("data/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:

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.

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:

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:

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:

pairs(total[-1], panel=panel.smooth)

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

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:

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.