Homework03

This homeork assignment is an excercise in data aggregation with a focus on writing code that is clean and readable.

# read from internet gDat <-
# read.delim('http://www.stat.ubc.ca/~jenny/notOcto/STAT545A/examples/gapminder/data/gapminderDataFiveYear.txt')

# read locally
gDat <- read.delim(file.path(str_replace(getwd(), "stat545a-2013-hw03", "GapMinder"), 
    "gapminderDataFiveYear.txt"))

I begin with the “harder” sections. I hope this saves the TA some time in marking.

Life Expectancy by Year

Mean versus trimmed mean

We compare the mean life expectancy by year to the trimmed mean (20% from each end).

mean.lifeExp_by_year <- arrange(ddply(gDat, ~year, summarize, mean = mean(lifeExp), 
    trimmed.mean = mean(lifeExp, trim = 0.2)), year)

print(xtable(mean.lifeExp_by_year, digits = 2), type = "html", include.rownames = F)

Curiously, the mean value of earlier years are higher than the trimmed mean, while the reverse is true for the later years…

(Challenge) Life expectancy over time on different continents: the “Wide” format

### performing the task with 1 line of code
mean.lifeExp_by_year_by_cont.wide <- daply(gDat, ~year, daply, ~continent, summarize, 
    mean = mean(lifeExp), median = median(lifeExp))  #note multiple stats can be computed with one call!

This is an odd excercise, that involves essentially “nested” daply() calls. The “inner” call takes data frame for a given year, returns a “array/matrix” of 2 dimensions: 1st dimension the factor levels you split by, and 2nd dimensions the stats you make it calculate. The “outer” call takes the resulting “array/matrix” and inserts another dimension to this array representing the year.

### data structure: array of 3 dimensions
class(mean.lifeExp_by_year_by_cont.wide)

## [1] "array"

dim(mean.lifeExp_by_year_by_cont.wide)

## [1] 12  5  2

To print a 2 dimensional table, one would use the multidimensional array indexing, [,…,]. As an example myMultDimArray[,,“statName”] retrieves the matrix with information on a statistic called “statName”.

### print 2d tables of 'mean'
print(xtable(mean.lifeExp_by_year_by_cont.wide[, , "mean"], digits = 2), type = "html", 
    include.rownames = T)

### print 2d table of 'median', with table caption
print(xtable(mean.lifeExp_by_year_by_cont.wide[, , "median"], digits = 2, caption = "Median life expectancy of 5 continents by year"), 
    type = "html", include.rownames = T, caption.placement = "top")

(Basically the same Challenge) Report the absolute and/or relative abundance of countries with low life expectancy over time by continent. “Wide” format. Extension of another task of unknown difficulty / value.

This task becomes very easy with the use of nested daply()!

(thresh <- quantile(gDat$lifeExp, 0.1))

##   10% 
## 41.51

(low_lifeExp <- daply(gDat, ~year, daply, ~continent, summarize, countLowLifeExp = sum(lifeExp < 
    thresh), propLowLifeExp = mean(lifeExp < thresh)))  #note that all of these stats can be computed in one call!

## , ,  = countLowLifeExp
## 
##       continent
## year   Africa Americas Asia Europe Oceania
##   1952 35     2        11   0      0      
##   1957 29     1        8    0      0      
##   1962 22     0        4    0      0      
##   1967 14     0        3    0      0      
##   1972 10     0        3    0      0      
##   1977 5      0        2    0      0      
##   1982 3      0        1    0      0      
##   1987 3      0        1    0      0      
##   1992 5      0        0    0      0      
##   1997 4      0        0    0      0      
##   2002 4      0        0    0      0      
##   2007 1      0        0    0      0      
## 
## , ,  = propLowLifeExp
## 
##       continent
## year   Africa  Americas Asia    Europe Oceania
##   1952 0.6731  0.08     0.3333  0      0      
##   1957 0.5577  0.04     0.2424  0      0      
##   1962 0.4231  0        0.1212  0      0      
##   1967 0.2692  0        0.09091 0      0      
##   1972 0.1923  0        0.09091 0      0      
##   1977 0.09615 0        0.06061 0      0      
##   1982 0.05769 0        0.0303  0      0      
##   1987 0.05769 0        0.0303  0      0      
##   1992 0.09615 0        0       0      0      
##   1997 0.07692 0        0       0      0      
##   2002 0.07692 0        0       0      0      
##   2007 0.01923 0        0       0      0

Note that printing the object results in this “slicing” of multidimensional arrays by the inner-most dimension- in our case the dimension representing statistics we created.

This task actually embeds my knowledge of how to create counts and proportions using comparisons/masks (e.g. propLowLifeExp= mean(lifeExp < thresh)). It is probably unnecessary to print so many tables for 1 assignment. saves times for the marker.

### print 2d tables of 'mean'
print(xtable(low_lifeExp[, , "countLowLifeExp"], digits = 0, caption = "The counts of countries with low life expectancy in each continent by year"), 
    type = "html", include.rownames = T, caption.placement = "top")

### print 2d table of 'median', with table caption
print(xtable(low_lifeExp[, , "propLowLifeExp"], digits = 2, caption = "The proportion of countries with low life expectancy in each continent by year"), 
    type = "html", include.rownames = T, caption.placement = "top")

GDP per capita by Continent

A wide table

Create a table with columns “Continent”, “Minimum GDP per Capita” and “Maximum GDP per Capita”. The continents are arranged by ascending minimum GDP.

minmax.gdp_by_continent.wide <- arrange(ddply(gDat, ~continent, summarize, minGDP = min(gdpPercap), 
    maxGDP = max(gdpPercap)), minGDP)
print(xtable(minmax.gdp_by_continent.wide, digits = 0), type = "html", include.rownames = F)

A Tall table

An alternative look at minimum and maximum GDP per capita by continent. Here the “type” of statistic occupies its one column followed by the results of applying that summary statistic.

# this is not very elegant
fct <- function(dframe) {
    res <- data.frame(factor_name = "min", GDP = min(dframe$gdpPercap))
    res <- rbind(res, data.frame(factor_name = "max", GDP = max(dframe$gdpPercap)))
}

minmax.gdp_by_continent.tall <- arrange(ddply(gDat, ~continent, fct), factor_name, 
    GDP)
print(xtable(minmax.gdp_by_continent.tall, digits = 0), type = "html", include.rownames = F)

Spread of GDP per capita

Here we investigate the spread of GDP per capita in each continent by 3 different measures, the standard deviation (sd), median absolute deviation (mad), and interquartile range (IQR).

spread.gdp_by_continent <- arrange(ddply(gDat, ~continent, summarize, sd = sd(gdpPercap), 
    mad = mad(gdpPercap), IQR = IQR(gdpPercap)), IQR)

print(xtable(spread.gdp_by_continent, digits = 0), type = "html", include.rownames = F)

The above table is sorted by ascending IQR values. Notice that ordering by standard deviations is very different from this ordering by IQR, while the order by mad follows the order of IQR more closely. This suggests that sd may not always be the most meaningful spread to report and interpret for strangely shaped distributions. One often interprets sd as width containing 60% of the data assuming normal distribution, but for highly skewed data this interpretation is far from correct. An IQR, based on actual proporions from the data may be a more useful measure of spread.