I keep my eye out for blog posts illustrating data analysis tasks in R that I can use to test the functionality of my chez-stats and dataframe libraries for Scheme (R6RS). A post comparing pandas (Python) and dplyr (R) in a basic analysis of the gapminder dataset provides a nice little test case. In this post, I will also include base R code used to accomplish the same tasks as a contrast to both the Scheme code and the dplyr code from the other post.

Data Loading

The R package, gapminder, provides an excerpt of the data available at Gapminder.org. I've written the data from that package to a CSV file (available here).

Base R

I've added a little function, head10, to simplify subsequent code. It is not necessary. By default, head provides six rows. I've changed it to 10 to match the default for dataframe-display. The main thing to note here is that because we are using base R, we don't need to load any packages.

> gapminder <- read.csv("gapminder.csv")
> head10 <- function(data) head(data, n = 10)
> head10(gapminder)

   continent     country year lifeExp      pop gdpPercap
1       Asia Afghanistan 1952  28.801  8425333  779.4453
2       Asia Afghanistan 1957  30.332  9240934  820.8530
3       Asia Afghanistan 1962  31.997 10267083  853.1007
4       Asia Afghanistan 1967  34.020 11537966  836.1971
5       Asia Afghanistan 1972  36.088 13079460  739.9811
6       Asia Afghanistan 1977  38.438 14880372  786.1134
7       Asia Afghanistan 1982  39.854 12881816  978.0114
8       Asia Afghanistan 1987  40.822 13867957  852.3959
9       Asia Afghanistan 1992  41.674 16317921  649.3414
10      Asia Afghanistan 1997  41.763 22227415  635.3414

Scheme

The code here is more verbose than in base R. We need to import a couple of libraries. Also, read-delim reads data row-wise, which I'm calling a rowtable, and needs to be converted to a column-wise dataframe. rowtable->dataframe does minimal type conversion; only converting strings to numbers.

(import (chez-stats)
        (dataframe))

(define gapminder
  (-> (read-delim "gapminder.csv")
      (rowtable->dataframe)))

> (dataframe-display gapminder)

 dim: 1704 rows x 6 cols
  continent      country   year  lifeExp       pop  gdpPercap 
       Asia  Afghanistan  1952.  28.8010  8.425E+6   779.4453 
       Asia  Afghanistan  1957.  30.3320  9.241E+6   820.8530 
       Asia  Afghanistan  1962.  31.9970  1.027E+7   853.1007 
       Asia  Afghanistan  1967.  34.0200  1.154E+7   836.1971 
       Asia  Afghanistan  1972.  36.0880  1.308E+7   739.9811 
       Asia  Afghanistan  1977.  38.4380  1.488E+7   786.1134 
       Asia  Afghanistan  1982.  39.8540  1.288E+7   978.0114 
       Asia  Afghanistan  1987.  40.8220  1.387E+7   852.3959 
       Asia  Afghanistan  1992.  41.6740  1.632E+7   649.3414 
       Asia  Afghanistan  1997.  41.7630  2.223E+7   635.3414 

Filtering

Problem 1

Filter the dataset to retain only rows where year is 2007.

Base R

> head10(gapminder[gapminder$year == 2007, ])

    continent     country year lifeExp       pop  gdpPercap
12       Asia Afghanistan 2007  43.828  31889923   974.5803
24     Europe     Albania 2007  76.423   3600523  5937.0295
36     Africa     Algeria 2007  72.301  33333216  6223.3675
48     Africa      Angola 2007  42.731  12420476  4797.2313
60   Americas   Argentina 2007  75.320  40301927 12779.3796
72    Oceania   Australia 2007  81.235  20434176 34435.3674
84     Europe     Austria 2007  79.829   8199783 36126.4927
96       Asia     Bahrain 2007  75.635    708573 29796.0483
108      Asia  Bangladesh 2007  64.062 150448339  1391.2538
120    Europe     Belgium 2007  79.441  10392226 33692.6051

Scheme

This example introduces the thread-first operator (->), which takes the result of the previous procedure and passes it to the first argument of the next procedure. For data analysis, I strongly prefer the threading (or piping) approach over writing the code inside-out or creating lots of intermediate datasets.

> (-> gapminder
      (dataframe-filter* (year) (= year 2007))
      dataframe-display)

 dim: 142 rows x 6 cols
  continent      country   year  lifeExp       pop  gdpPercap 
       Asia  Afghanistan  2007.  43.8280  3.189E+7   974.5803 
     Europe      Albania  2007.  76.4230  3.601E+6  5937.0295 
     Africa      Algeria  2007.  72.3010  3.333E+7  6223.3675 
     Africa       Angola  2007.  42.7310  1.242E+7  4797.2313 
   Americas    Argentina  2007.  75.3200  4.030E+7 12779.3796 
    Oceania    Australia  2007.  81.2350  2.043E+7 34435.3674 
     Europe      Austria  2007.  79.8290  8.200E+6 36126.4927 
       Asia      Bahrain  2007.  75.6350  7.086E+5 29796.0483 
       Asia   Bangladesh  2007.  64.0620  1.504E+8  1391.2538 
     Europe      Belgium  2007.  79.4410  1.039E+7 33692.6051 

Problem 2

Filter the dataset to retain only rows where year is 2007 and continent is Americas.

Base R

In this case, I've used subset rather than the more conventional [ subsetting used in Problem 1. The advantage of subset is that I only need to type gapminder once instead of three times.

> head10(subset(gapminder, year == 2007 & continent == "Americas"))

    continent            country year lifeExp       pop gdpPercap
60   Americas          Argentina 2007  75.320  40301927 12779.380
144  Americas            Bolivia 2007  65.554   9119152  3822.137
180  Americas             Brazil 2007  72.390 190010647  9065.801
252  Americas             Canada 2007  80.653  33390141 36319.235
288  Americas              Chile 2007  78.553  16284741 13171.639
312  Americas           Colombia 2007  72.889  44227550  7006.580
360  Americas         Costa Rica 2007  78.782   4133884  9645.061
396  Americas               Cuba 2007  78.273  11416987  8948.103
444  Americas Dominican Republic 2007  72.235   9319622  6025.375
456  Americas            Ecuador 2007  74.994  13755680  6873.262

Scheme

dataframe-filter* was designed to avoid having to repeat the dataframe name in each sub-expression, i.e., (= year 2007) and (string=? continent "Americas"), but requires typing year and continent two times each in this example.

> (-> gapminder
      (dataframe-filter*
       (year continent)
       (and (= year 2007)
            (string=? continent "Americas")))
      dataframe-display)

 dim: 25 rows x 6 cols
  continent             country   year  lifeExp       pop  gdpPercap 
   Americas           Argentina  2007.  75.3200  4.030E+7 12779.3796 
   Americas             Bolivia  2007.  65.5540  9.119E+6  3822.1371 
   Americas              Brazil  2007.  72.3900  1.900E+8  9065.8008 
   Americas              Canada  2007.  80.6530  3.339E+7 36319.2350 
   Americas               Chile  2007.  78.5530  1.628E+7 13171.6389 
   Americas            Colombia  2007.  72.8890  4.423E+7  7006.5804 
   Americas          Costa Rica  2007.  78.7820  4.134E+6  9645.0614 
   Americas                Cuba  2007.  78.2730  1.142E+7  8948.1029 
   Americas  Dominican Republic  2007.  72.2350  9.320E+6  6025.3748 
   Americas             Ecuador  2007.  74.9940  1.376E+7  6873.2623 

Problem 3

Filter the dataset to retain only rows where year is 2007 and continent is Americas and country is United States. This last filter is a little silly because the country condition makes the continent condition unnecessary, but the point is to show how code complexity increases with additional conditions.

Base R

Thanks to subset, adding additional conditions is straightforward with zero redundancy.

> subset(gapminder, 
         year == 2007 & 
           continent == "Americas" &
           country == "United States")

     continent       country year lifeExp       pop gdpPercap
1620  Americas United States 2007  78.242 301139947  42951.65

Scheme

> (-> gapminder
      (dataframe-filter*
       (year continent country)
       (and (= year 2007)
            (string=? continent "Americas")
            (string=? country "United States")))
      dataframe-display)

 dim: 1 rows x 6 cols
  continent        country   year  lifeExp       pop  gdpPercap 
   Americas  United States  2007.  78.2420  3.011E+8   42951.65 

Summary Statistics

Problem 1

Calculate the average life expectancy worldwide in 2007.

Base R

> mean(subset(gapminder, year == 2007)$lifeExp)
[1] 67.00742

Scheme

This example introduces the $ operator, which was inspired by R, to extract the values of a column (e.g., lifeExp). The difference in verbosity between base R and Scheme in this example is related to the filtering step. Extracting a column and calculating the mean involve nearly the same number of characters. mean is provided by chez-stats, but is simple procedure. (mean x) is equivalent to (/ (apply + x) (length x)).

> (-> gapminder
      (dataframe-filter* (year) (= year 2007))
      ($ 'lifeExp)
      (mean))

67.00742253521126    

Problem 2

Calculate the average life expectancy for every continent in 2007.

Base R

aggregate allows for use of the formula syntax (e.g., lifeExp ~ continent) to concisely describe the summarized value (lifeExp) and grouping variable(s).

> aggregate(lifeExp ~ continent, 
           data = subset(gapminder, year == 2007),
           FUN = mean)

  continent  lifeExp
1    Africa 54.80604
2  Americas 73.60812
3      Asia 70.72848
4    Europe 77.64860
5   Oceania 80.71950

Scheme

One difference to note here is that dataframe-aggregate doesn't automatically sort by the grouping variable. We would have to explicitly add that sorting step.

> (-> gapminder
      (dataframe-filter* (year) (= year 2007))
      (dataframe-aggregate*
       (continent)
       (mean-lifeExp (lifeExp) (mean lifeExp)))
      (dataframe-display))

 dim: 5 rows x 2 cols
  continent  mean-lifeExp 
    Oceania       80.7195 
     Europe       77.6486 
   Americas       73.6081 
       Asia       70.7285 
     Africa       54.8060 

Problem 3

Calculate the total population per continent in 2007 and sort the results in descending order of total population.

Base R

> total_pop = aggregate(pop ~ continent, 
                        data = subset(gapminder, year == 2007),
                        FUN = sum)
> total_pop[order(-total_pop$pop),]

  continent        pop
3      Asia 3811953827
1    Africa  929539692
2  Americas  898871184
4    Europe  586098529
5   Oceania   24549947

Scheme

> (-> gapminder
      (dataframe-filter* (year) (= year 2007))
      (dataframe-aggregate*
       (continent)
       (total-pop (pop) (sum pop)))
      (dataframe-sort* (> total-pop))
      (dataframe-display))

 dim: 5 rows x 2 cols
  continent  total-pop 
       Asia   3.812E+9 
     Africa   9.295E+8 
   Americas   8.989E+8 
     Europe   5.861E+8 
    Oceania   2.455E+7 

Creating Derived Columns

Problem 1

Calculate the total GDP by multiplying pop and gdpPercap.

Base R

The help page for transform advises that it is only intended for interactive use. Alternatively, you could use: gapminder$GPD = gapminder$pop * gapminder$gdpPercap.

> head10(transform(gapminder, GDP = pop * gdpPercap))

   continent     country year lifeExp      pop gdpPercap         GDP
1       Asia Afghanistan 1952  28.801  8425333  779.4453  6567086330
2       Asia Afghanistan 1957  30.332  9240934  820.8530  7585448670
3       Asia Afghanistan 1962  31.997 10267083  853.1007  8758855797
4       Asia Afghanistan 1967  34.020 11537966  836.1971  9648014150
5       Asia Afghanistan 1972  36.088 13079460  739.9811  9678553274
6       Asia Afghanistan 1977  38.438 14880372  786.1134 11697659231
7       Asia Afghanistan 1982  39.854 12881816  978.0114 12598563401
8       Asia Afghanistan 1987  40.822 13867957  852.3959 11820990309
9       Asia Afghanistan 1992  41.674 16317921  649.3414 10595901589
10      Asia Afghanistan 1997  41.763 22227415  635.3414 14121995875

Scheme

> (-> gapminder
      (dataframe-modify*
       (GDP (pop gdpPercap) (* pop gdpPercap)))
      (dataframe-display))

 dim: 1704 rows x 7 cols
  continent      country   year  lifeExp       pop  gdpPercap        GDP 
       Asia  Afghanistan  1952.  28.8010  8.425E+6   779.4453  6.567E+09 
       Asia  Afghanistan  1957.  30.3320  9.241E+6   820.8530  7.585E+09 
       Asia  Afghanistan  1962.  31.9970  1.027E+7   853.1007  8.759E+09 
       Asia  Afghanistan  1967.  34.0200  1.154E+7   836.1971  9.648E+09 
       Asia  Afghanistan  1972.  36.0880  1.308E+7   739.9811  9.679E+09 
       Asia  Afghanistan  1977.  38.4380  1.488E+7   786.1134  1.170E+10 
       Asia  Afghanistan  1982.  39.8540  1.288E+7   978.0114  1.260E+10 
       Asia  Afghanistan  1987.  40.8220  1.387E+7   852.3959  1.182E+10 
       Asia  Afghanistan  1992.  41.6740  1.632E+7   649.3414  1.060E+10 
       Asia  Afghanistan  1997.  41.7630  2.223E+7   635.3414  1.412E+10 

Problem 2

Find the top 10 countries in percentile of gdpPercap.

Base R

> percentile <- function(x){
    rank_x = rank(x)
    rank_x/max(rank_x)
  }

> gapminder2007 = transform(subset(gapminder, year == 2007),
                            percentile = percentile(gdpPercap))
> head10(gapminder2007[order(-gapminder2007$percentile),])

     continent          country year lifeExp       pop gdpPercap percentile
1152    Europe           Norway 2007  80.196   4627926  49357.19  1.0000000
864       Asia           Kuwait 2007  77.588   2505559  47306.99  0.9929577
1368      Asia        Singapore 2007  79.972   4553009  47143.18  0.9859155
1620  Americas    United States 2007  78.242 301139947  42951.65  0.9788732
756     Europe          Ireland 2007  78.885   4109086  40676.00  0.9718310
672       Asia Hong Kong, China 2007  82.208   6980412  39724.98  0.9647887
1488    Europe      Switzerland 2007  81.701   7554661  37506.42  0.9577465
1092    Europe      Netherlands 2007  79.762  16570613  36797.93  0.9507042
252   Americas           Canada 2007  80.653  33390141  36319.24  0.9436620
696     Europe          Iceland 2007  81.757    301931  36180.79  0.9366197

Scheme

This example reveals a weak spot in the dataframe API. The issue is that rank% operates on a list whereas dataframe-modify* takes only scalars as inputs because the expr is mapped over all rows in a dataframe. The workaround to avoid the mapping is to specify that no columns (i.e., ()) from the dataframe are used, then the list created in the subsequent expression (e.g., (rank% ($ df 'gdpPercap))) will be used as the column values.

Because the dataframe library only has thread-first (->) and thread-last (->>) operators, we have to create the awkward lambda procedure to get the output of dataframe-filter* to the right place in dataframe-modify*. SRFI 197 provides a pipeline operator that requires that the location of the output of the previous procedure is explicitly specified with _, which would likely simplify this code. An earlier, simpler version of SRFI 197 was the source for -> and ->> in the dataframe library.

(define (rank% lst)
  (let* ([rank-lst (rank lst 'mean)]
         [max-rank (apply max rank-lst)])
    (map (lambda (x) (exact->inexact (/ x max-rank))) rank-lst)))

> (-> gapminder
      (dataframe-filter* (year) (= year 2007))
      (->> ((lambda (df)
              (dataframe-modify*
               df
               (percentile () (rank% ($ df 'gdpPercap)))))))
      (dataframe-sort* (> gdpPercap))
      (dataframe-display))

 dim: 142 rows x 7 cols
  continent           country   year  lifeExp       pop  gdpPercap  percentile 
     Europe            Norway  2007.  80.1960  4.628E+6   49357.19      1.0000 
       Asia            Kuwait  2007.  77.5880  2.506E+6   47306.99      0.9930 
       Asia         Singapore  2007.  79.9720  4.553E+6   47143.18      0.9859 
   Americas     United States  2007.  78.2420  3.011E+8   42951.65      0.9789 
     Europe           Ireland  2007.  78.8850  4.109E+6   40676.00      0.9718 
       Asia  Hong Kong, China  2007.  82.2080  6.980E+6   39724.98      0.9648 
     Europe       Switzerland  2007.  81.7010  7.555E+6   37506.42      0.9577 
     Europe       Netherlands  2007.  79.7620  1.657E+7   36797.93      0.9507 
   Americas            Canada  2007.  80.6530  3.339E+7   36319.24      0.9437 
     Europe           Iceland  2007.  81.7570  3.019E+5   36180.79      0.9366 

If we create an intermediate dataframe, we get simpler syntax.

(define gapminder2007
  (dataframe-filter* gapminder (year) (= year 2007)))

(-> (dataframe-modify*
     gapminder2007
     (percentile () (rank% ($ gapminder2007 'gdpPercap))))
    (dataframe-sort* (> gdpPercap))
    (dataframe-display))

Conclusion

The main outcome of writing this post was that it led me to completely rewrite dataframe-display because it previously didn't handle large numbers well, which was made abundantly clear in working with the gapminder data. The new version of dataframe-display is much improved and I gained a greater appreciation for all of the decisions that go into how to print a representation of a data structure to the screen.

If you are not familiar with Scheme code, you might find it to be unacceptably verbose in the examples above and, especially, when compared to dplyr code. Because of Scheme's macro system, dataframe could be written in a more terse style. But I made the decision early on to stick to relatively simple macro usage and write dataframe in a way that I thought would be familiar to Scheme programmers. And that often involves more verbose code. For example, nearly all of the dataframe procedures have dataframe in the name, e.g., dataframe-filter**, dataframe-modify*, etc. This is following the example for hashtables, e.g., hashtable-ref, hashtable-values, etc. I also hope that experienced Scheme programmers can see that the dataframe macros mostly exist to reduce the number of times that lambda is written but the shape of the code should still feel familiar. While Scheme code might be more verbose in the small, I find it extremely expressive in the large because the core ideas compose so well.