Stat545A Homework #3 Data aggregation

In this homework assignment I will continue to leverage Hadley's impressive collection of R toolboxes. In particular, I will examine the Gapminder data with plyr. Unfortunately, no figures will be used here–it is very painful to try to examine trends without graphics–hopefully this pain is evident in the tables below. Just for fun, I made the tables a bit easier to read, but this should be no solace at all in the face of no figures!

Let's begin:

Basic import and quality control

library(plyr)
library(xtable)
library(MASS)
gDat <- read.delim("gapminderDataFiveYear.txt")
str(gDat)

## 'data.frame':    1704 obs. of  6 variables:
##  $ country  : Factor w/ 142 levels "Afghanistan",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ year     : int  1952 1957 1962 1967 1972 1977 1982 1987 1992 1997 ...
##  $ pop      : num  8425333 9240934 10267083 11537966 13079460 ...
##  $ continent: Factor w/ 5 levels "Africa","Americas",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ lifeExp  : num  28.8 30.3 32 34 36.1 ...
##  $ gdpPercap: num  779 821 853 836 740 ...

tail(gDat)

##       country year      pop continent lifeExp gdpPercap
## 1699 Zimbabwe 1982  7636524    Africa   60.36     788.9
## 1700 Zimbabwe 1987  9216418    Africa   62.35     706.2
## 1701 Zimbabwe 1992 10704340    Africa   60.38     693.4
## 1702 Zimbabwe 1997 11404948    Africa   46.81     792.4
## 1703 Zimbabwe 2002 11926563    Africa   39.99     672.0
## 1704 Zimbabwe 2007 12311143    Africa   43.49     469.7

Since I will be making many tables, here is a utility function to make the process (a little) less painful:

printTable <- function(df) {
    print(xtable(df), type = "html", include.rownames = F)
}

Easy analysis of GDP “trends”

I am going to explore the GDP per capita data extensively in the next few steps. At first we will look to see the min and max GDP per capita on a continent by continent basis.

minMaxGDPwide <- ddply(gDat, ~continent, summarize, minGDP = min(gdpPercap), 
    maxGDP = max(gdpPercap))
minMaxGDPwide <- arrange(minMaxGDPwide, minGDP)  # sort on minGDP
printTable(minMaxGDPwide)

Nothing is too surprising here: as expected Africa and China have the smallest GDP per capita while Europe, Americas, and Oceania lead the group. Interestingly, the trends do not carry onto the maxGDP by continent–here we observe that Asia actually has the largest maxGDP. An obvious question is: for which country and in what year did these min and max values occur? We look into this in the below:

minMaxGDPElements <- ddply(gDat, ~continent, summarize, countryMin = country[which.min(gdpPercap)], 
    timeMin = year[which.min(gdpPercap)], minGDP = min(gdpPercap), countryMax = country[which.max(gdpPercap)], 
    timeMax = year[which.max(gdpPercap)], maxGDP = max(gdpPercap))
minMaxGDPElements <- arrange(minMaxGDPElements, minGDP)
printTable(minMaxGDPElements)

A few disturbing cases pop out here:

• Africa's max GDP occurred in 1977, and in as recently as 2002 Congo achieved the lowest per capita GDP in Africa's history.
• In the Americas, both the highest and lowest GDP per capita were observed in 2007 (!), by America and Haiti, respectively.
• Australia has clearly experienced a higher quality of life over time.

What is the spread of GDP across continents?

spreadGDP <- ddply(gDat, ~continent, summarize, SD = sd(gdpPercap), MAD = mad(gdpPercap), 
    IQR = IQR(gdpPercap))
spreadGDP <- arrange(spreadGDP, SD)
printTable(spreadGDP)

In looking at the spread data, I have examined three different measures of spread:

• Standard deviation
• Median average deviation
• Interquartile range

The second two measures are meant to be less sensitive to outlines, and thus more robust. Interestingly the ordering in GDP spread is not invariant to measure. Most notably, I would want to point out that the dispersion in GDP per capita in the Americas seems much smaller on a MAD basis than on a standard deviation basis.

Life expectancy trends

For some preliminary results on life expectancy, I looked to see the average life expectancy as a function of year. In the table below, we see that average life expectancy has increased through time. To make the measure more robust, I also examined the trimmed mean–that is, the mean of a trimmed dataset. In this instance the data has been trimmed to exclude the 10% of the data at the ends of the ordered life expectancy set.

# Trimmed mean, trim describes the fraction of observations to be trimmed at
# each end--btwn 0 - 0.5
lifeExpTrim <- ddply(gDat, ~year, summarize, StandardMean = mean(lifeExp), TrimmedMean = mean(lifeExp, 
    trim = 0.1), `StandardMean - TrimmedMean` = StandardMean - TrimmedMean)
printTable(lifeExpTrim)

In the earlier years (1952-1962) the standard (untrimmed) mean is greater than the trimmed mean which indicates that there were countries with very high life expediencies as outliers to the whole data set. This trend reversed in 1967 which indicates that the outliers are now pulling the standard mean down.

It may also be of interest to see how life expectancy is changing over time on different continents:

lifeExpCont <- ddply(gDat, .(continent, year), summarize, meanLifeExp = mean(lifeExp), 
    medianLifeExp = median(lifeExp))
# How is life expectancy changing over time on different continents? 'Wide'
# format
lifeExpContWide <- ddply(lifeExpCont, ~year, function(t) setNames(t$meanLifeExp, 
    unique(t$continent)))
printTable(lifeExpCont)

printTable(lifeExpContWide)

The long data is presented first, and then is converted to a wide format through plyr. Note, however, that this is really a reshaping through plyr (which could have been done through reshape2) and not a native plyr command to generate wide data. In the actual data, it is encouraging to see that within each continent the general trend has been increased life expectancy, though clearly at different rates. This is consistent with our earlier findings on the trimmed mean!

How many of these continents, and in turn, countries have low life expectancy? We define low life expectancy to be less than or equal to the 5% quantile.

lowLifeExp <- as.numeric(quantile(gDat$lifeExp, probs = 0.05))
continentLifeExp <- ddply(gDat, .(continent, year), summarize, lowLifeInstances = sum(lifeExp <= 
    lowLifeExp))  # this is tall, aka not fun
continentLifeExp <- ddply(continentLifeExp, ~year, function(t) setNames(t$lowLifeInstances, 
    unique(t$continent)))
printTable(continentLifeExp)

I noticed immediately that the low life expectancy data was predominately in Africa and Asia and stopped in 1997 (though this is really an artifact of how a low life expectancy is defined). It is more informative to have the proportion of countries which have low life expectancy than just the raw counts:

continentLifeExpProp <- ddply(gDat, .(continent, year), summarize, lowLifeInstances = sum(lifeExp <= 
    lowLifeExp)/length(unique(country)))
continentLifeExpProp <- ddply(continentLifeExpProp, ~year, function(t) setNames(round(t$lowLifeInstances, 
    2) * 100, unique(t$continent)))
printTable(continentLifeExpProp)

Wow, in 1952 44% of Africa's countries had a “low life expectancy!”

Interesting story lines

• Following on an idea from class, I have run a linear regression of life expectancy vs. time for each country and examined the minimum coefficients

yearMin <- min(gDat$year)
lifeExpRegression <- function(x) {
    estCoefs <- coef(lm(lifeExp ~ I(year - yearMin), x))
    names(estCoefs) <- c("intercept", "slope")
    return(estCoefs)
}

lifeExpChanges <- arrange(ddply(ddply(gDat, .(continent, country), lifeExpRegression), 
    ~continent, summarize, minInt = min(intercept), minIntCountry = country[which.min(intercept)], 
    minSlope = min(slope), minSlopeCountry = country[which.min(slope)]), minInt)
printTable(lifeExpChanges)

Slightly alarming is that Zimbabwe is actually experiencing deteriorating life expectancy.

We can modify our regression function to provide us with some interesting results based on residuals from the least-squares fit:

lifeExpRegression <- function(x) {
    fit <- lm(lifeExp ~ I(year - yearMin), x)
    res <- max(fit$residuals)
    names(res) <- c("maxResidual")
    return(res)
}

lifeExpRes <- arrange(ddply(ddply(gDat, .(continent, country), lifeExpRegression), 
    ~continent, summarize, maxRes = max(maxResidual), maxResCountry = country[which.max(maxResidual)]), 
    maxRes)
printTable(lifeExpRes)

It would be interesting to do some analysis of these 5 countries to see why the life expectancy trend broke down…But I need figures and that has to wait until next week :)