```Econometrics and Free Software by Bruno Rodrigues.
Check out my package that adds logging to R functions, {chronicler}.
Or read my free ebooks, to learn some R and build reproducible analytical pipelines..
You can also watch my youtube channel or find the slides to the talks I've given here.
Buy me a coffee, my kids don't let me sleep.
```

# Why you should(n't) care about Monads if you're an R programmer

R

Update: I also made a video out of this blog post; watch it on youtube.

## Introduction: functions

To understand Monads, I think it’s useful to first think about functions; why do we use functions? Why don’t we simply write scripts with the required operations one after the other? For instance, to compute the average height by species in a data set of individuals from the famous space opera “Star Wars”, we could very well write this code:

``````suppressPackageStartupMessages(library(dplyr))

data(starwars)

sum_humans <- 0
sum_others <- 0
n_humans <- 0
n_others <- 0

for(i in seq_along(1:nrow(starwars))){

if(!is.na(unlist(starwars[i, "species"])) &
unlist(starwars[i, "species"]) == "Human"){
if(!is.na(unlist(starwars[i, "height"]))){
sum_humans <- sum_humans + unlist(starwars[i, "height"])
n_humans <- n_humans + 1
} else {

0

}

} else {
if(!is.na(unlist(starwars[i, "height"]))){
sum_others <- sum_others + unlist(starwars[i, "height"])
n_others <- n_others + 1
} else {
0
}
}
}

mean_height_humans <- sum_humans/n_humans
mean_height_others <- sum_others/n_others``````

Well, we could do it like this, but we definitely shouldn’t:

• what this code does is not immediately obvious. If the code blocks aren’t commented, readers of this code will have to read line by line to understand what is going on;
• this code is not reusable. If now I need the average height by species and sex, I need to copy and paste the code, and modify it, and in some cases modify it substantially;
• this code handles missing values in a cumbersome way, with nested `if...else...`s;
• this code is not easy to test;
• this code cannot be composed (meaning, chained) with other code without substantially altering it (to be precise, chaining and composing are two different things, strictly speaking, but for simplicity’s sake, let’s just assume it is the same. Whenever I’m talking about “composing” something, I mean “chaining” something.)

But it’s not just shortcomings, this imperative code has one advantage; it uses only some very fundamental building blocks: `if...else...`, for loops and that’s almost it (it does use some functions provided by a base installation of R, namely `is.na()`, `!()`, `unlist()` and `[()`, so strictly speaking, the code above is not purely imperative, but maybe closer to being procedural?).

Using functions solves all the issues from imperative programming. Here is a base solution to the problem above, using a declarative, or functional, approach:

``````aggregate(starwars\$height,
by = list(starwars\$species == "Human"),
FUN = \(x)(mean(x, na.rm = TRUE)))``````
``````##   Group.1        x
## 1   FALSE 172.4043
## 2    TRUE 176.6452``````

• what this code does is obvious, but only if you know what `aggregate()` does. But if you read its documentation you’ll know, and you’ll know every time you’ll see `aggregate()` unlike a loop like the loop above where you’ll have to read it each time to understand;
• this code is reusable. Replace the data frame by another, and that’s it;
• Missing values are now ignored easily using the `na.rm` argument of `mean()`;
• this code is easy to test (using unit tests);
• this code can be composed, for instance like this:
``````aggregate(starwars\$height,
by = list(starwars\$species == "Human"),
FUN = \(x)(mean(x, na.rm = TRUE))) |>
setNames(c("is_human", "mean_height"))``````
``````##   is_human mean_height
## 1    FALSE    172.4043
## 2     TRUE    176.6452``````

The issue with the functional approach (at least that’s the issue that many people I spoke to about this raise) is that… in some way people that don’t like this approach feel like they “lose” control over what’s going on. You don’t know what happens inside these functions. I remember, while working my first job, that my boss required that I don’t use any functions nor packages, but instead write all the loops explicitely, because she wanted to understand what was going on (of course, I completely ignored this request and just did as I pleased). As discussed above, the imperative approach requires minimum knowledge of the language, and almost anyone with an ounce of programming experience can understand imperative code. That’s not the case with a functional approach. Readers will have to be familiar with the individual functions like `aggregate()`, but also anonymous functions (I had to use `\(x)(mean(x, na.rm = TRUE))` to set `na.rm = TRUE`, which is `FALSE` by default) and also `|>` for composition/chaining.

It may same more complex, and maybe it is, but the advantages far outweigh the shortcoming.

For completeness, here is a `{dplyr}` version:

``````starwars %>%
group_by(is_human = species == "Human") %>%
summarise(mean_height = mean(height, na.rm = TRUE))``````
``````## # A tibble: 3 × 2
##   is_human mean_height
##   <lgl>          <dbl>
## 1 FALSE           172.
## 2 TRUE            177.
## 3 NA              181.``````

`{dplyr}` code is even more concise than base functional code. Here again, users will have to know about the individual functions and `%>%`. But personally, I think that the only hurdle is understanding what `%>%` does, and once you know this, `{dplyr}` code can be understood quite easily, thanks to very explicit function names.

So functions are great. They’re easy to test, easy to document, easy to package, easy to reuse, and easy to compose. Composition is really important. For example, let’s go back to the imperative code, and put the result in a neat data frame object, like the functional solutions do:

``````sum_humans <- 0
sum_others <- 0
n_humans <- 0
n_others <- 0

for(i in seq_along(1:nrow(starwars))){

if(!is.na(unlist(starwars[i, "species"])) &
unlist(starwars[i, "species"]) == "Human"){
if(!is.na(unlist(starwars[i, "height"]))){
sum_humans <- sum_humans + unlist(starwars[i, "height"])
n_humans <- n_humans + 1
} else {

0

}

} else {
if(!is.na(unlist(starwars[i, "height"]))){
sum_others <- sum_others + unlist(starwars[i, "height"])
n_others <- n_others + 1
} else {
0
}
}
}

mean_height_humans <- sum_humans/n_humans
mean_height_others <- sum_others/n_others

# These two lines are new
data.frame(list("is_human" = c(TRUE, FALSE),
"mean_height" = c(mean_height_others, mean_height_humans)))``````
``````##   is_human mean_height
## 1     TRUE    172.9400
## 2    FALSE    176.6452``````

It’s just two lines (right at the end), but the implications are huge; because imperative code cannot be composed, I had to write separate code to put the result into a data frame. More code that I need to write, more opportunities for mistakes. I actually did a mistake, did you notice? This kind of mistake could go unnoticed for eons. But if you use functions, you don’t have this problem, and can focus on getting (even complex) things done:

``````starwars %>%
filter(skin_color == "light") %>%
select(species, sex, mass) %>%
group_by(sex, species) %>%
summarise(
total_individuals = n(),
min_mass = min(mass, na.rm = TRUE),
mean_mass = mean(mass, na.rm = TRUE),
sd_mass = sd(mass, na.rm = TRUE),
max_mass = max(mass, na.rm = TRUE)
) %>%
select(-species) %>%
tidyr::pivot_longer(-sex, names_to = "statistic", values_to = "value")``````
``````## `summarise()` has grouped output by 'sex'. You can override using the `.groups`
## argument.``````
``````## # A tibble: 10 × 3
## # Groups:   sex [2]
##    sex    statistic         value
##    <chr>  <chr>             <dbl>
##  1 female total_individuals   6
##  2 female min_mass           45
##  3 female mean_mass          56.3
##  4 female sd_mass            16.3
##  5 female max_mass           75
##  6 male   total_individuals   5
##  7 male   min_mass           79
##  8 male   mean_mass          90.5
##  9 male   sd_mass            19.8
## 10 male   max_mass          120``````

Needless to say, trying to write the above code using only for loops and `if...else...` is not something I’d wish to do, especially passing the result of all the `{dplyr}` calls to `pivot_longer()`. Creating that last data frame by hand is error prone, and there would definitely be mistakes in there.

I hope I don’t need to convince you any more that functions are great, and that one of the great things they offer is their ability to be chained, or composed. But strictly speaking, you don’t need them. You could write your code without any function whatsoever, and use the most basic building blocks there are (loops and `if...else...` and little more). However, doing this would result in much messier code. It’s the same with monads. You can live without them. But there will be situations where not using them will result in messier code.

One more thing: as I was writing this blog post, I happened on this tweet:

This is a fine example of all that I’ve been discussing until now. The person who wrote this code was very likely trying to get the diagonal elements of a matrix. That person was likely a beginner in R and used for loops to try to get the answer. We have all been there; what I’m trying to articulate is this: imperative programming can be useful, but it can get messy very quickly…

## When functions are not enough

Functions are awesome, but there are situations which functions simply can’t easily deal with. Situations in which you would like your functions to do a little extra more, and the only way forward you see is to rewrite them to do something totally unrelated. For example, suppose you would like to time your code. Most people would to something such as:

``````tic <- Sys.time()
starwars %>%
filter(skin_color == "light") %>%
select(species, sex, mass) %>%
group_by(sex, species) %>%
summarise(
total_individuals = n(),
min_mass = min(mass, na.rm = TRUE),
mean_mass = mean(mass, na.rm = TRUE),
sd_mass = sd(mass, na.rm = TRUE),
max_mass = max(mass, na.rm = TRUE)
) %>%
select(-species) %>%
tidyr::pivot_longer(-sex, names_to = "statistic", values_to = "value")``````
``````## `summarise()` has grouped output by 'sex'. You can override using the `.groups`
## argument.``````
``````## # A tibble: 10 × 3
## # Groups:   sex [2]
##    sex    statistic         value
##    <chr>  <chr>             <dbl>
##  1 female total_individuals   6
##  2 female min_mass           45
##  3 female mean_mass          56.3
##  4 female sd_mass            16.3
##  5 female max_mass           75
##  6 male   total_individuals   5
##  7 male   min_mass           79
##  8 male   mean_mass          90.5
##  9 male   sd_mass            19.8
## 10 male   max_mass          120``````
``````toc <- Sys.time()

(running_time <- toc - tic)``````
``## Time difference of 0.04228544 secs``

You could totally do that. But now you’re back to square one. You have to deal with this tic-toc nonsense separately, have to keep track it, overburdening you mentally and polluting your code. To keep track of it, you’ll want to add the running times in a separate data frame, in which you could have all the running times of all your operations you need to run:

``````data.frame(list("operations" = seq(1:3),
"running_time" = c(running_time, running_time * 2, running_time * 3)))``````
``````##   operations    running_time
## 1          1 0.04228544 secs
## 2          2 0.08457088 secs
## 3          3 0.12685633 secs``````

This data frame is the consequence of this tic-toc nonsense not being composable and now you have to deal with it, but you don’t want to. So what now? You might be tempted to do something like this:

``````tic_filter <- function(...){

tic <- Sys.time()

result <- filter(...)

toc <- Sys.time()

message("Running time: ", toc - tic)

return(result)

}

starwars %>%
tic_filter(species == "Human")``````
``## Running time: 0.00481176376342773``
``````## # A tibble: 35 × 14
##    name     height  mass hair_color skin_color eye_color birth_year sex   gender
##    <chr>     <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr>
##  1 Luke Sk…    172    77 blond      fair       blue            19   male  mascu…
##  2 Darth V…    202   136 none       white      yellow          41.9 male  mascu…
##  3 Leia Or…    150    49 brown      light      brown           19   fema… femin…
##  4 Owen La…    178   120 brown, gr… light      blue            52   male  mascu…
##  5 Beru Wh…    165    75 brown      light      blue            47   fema… femin…
##  6 Biggs D…    183    84 black      light      brown           24   male  mascu…
##  7 Obi-Wan…    182    77 auburn, w… fair       blue-gray       57   male  mascu…
##  8 Anakin …    188    84 blond      fair       blue            41.9 male  mascu…
##  9 Wilhuff…    180    NA auburn, g… fair       blue            64   male  mascu…
## 10 Han Solo    180    80 brown      fair       brown           29   male  mascu…
## # … with 25 more rows, and 5 more variables: homeworld <chr>, species <chr>,
## #   films <list>, vehicles <list>, starships <list>``````

But that’s actually worse: not only do you have to change all the functions you need, and wrap them around tic-toc, but the running time is only shown as a message, so you can’t reuse it. You could then try to rewrite the function like this:

``````tic_filter <- function(...){

tic <- Sys.time()

result <- filter(...)

toc <- Sys.time()

running_time <- toc - tic

list("result" = result,
"running_time" = running_time)

}

starwars %>%
tic_filter(species == "Human")``````
``````## \$result
## # A tibble: 35 × 14
##    name     height  mass hair_color skin_color eye_color birth_year sex   gender
##    <chr>     <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr>
##  1 Luke Sk…    172    77 blond      fair       blue            19   male  mascu…
##  2 Darth V…    202   136 none       white      yellow          41.9 male  mascu…
##  3 Leia Or…    150    49 brown      light      brown           19   fema… femin…
##  4 Owen La…    178   120 brown, gr… light      blue            52   male  mascu…
##  5 Beru Wh…    165    75 brown      light      blue            47   fema… femin…
##  6 Biggs D…    183    84 black      light      brown           24   male  mascu…
##  7 Obi-Wan…    182    77 auburn, w… fair       blue-gray       57   male  mascu…
##  8 Anakin …    188    84 blond      fair       blue            41.9 male  mascu…
##  9 Wilhuff…    180    NA auburn, g… fair       blue            64   male  mascu…
## 10 Han Solo    180    80 brown      fair       brown           29   male  mascu…
## # … with 25 more rows, and 5 more variables: homeworld <chr>, species <chr>,
## #   films <list>, vehicles <list>, starships <list>
##
## \$running_time
## Time difference of 0.004878759 secs``````

At least now you save the running time along with the object. But the problem of rewriting many functions remains, and these rewritten `{dplyr}` functions now return a list, and not a data frame anymore so something like this:

``````starwars %>%
tic_filter(species == "Human") %>%
tic_select(species, sex)``````

wouldn’t work, because `tic_select()` expects a data frame, not a list where the first element is a data frame and the second a double.

So what else can be done? Perhaps you’d be tempted to use a global variable for this:

``````tic_filter <- function(..., running_time = 0){

tic <- Sys.time()

result <- filter(...)

toc <- Sys.time()

running_time <<- toc - tic + running_time

result

}``````

Functions written like this would save the running time in a global variable called `running_time` and each of them would take turns overwriting it:

``````running_time <- 0

one <- starwars %>%
tic_filter(species == "Human", running_time = running_time)

running_time``````
``## Time difference of 0.00490284 secs``
``````two <- one %>%
tic_select(species, sex, running_time = running_time)

running_time``````
``## Time difference of 0.007258415 secs``

(I defined `tic_select()` but am not showing it here.)

This has the advantage that the wrapped functions now return data frames as well, and can thus be composed/chained. But these functions are not pure functions, because they change something (the global variable `running_time`) outside their scope. Impure functions can be tricky; for instance here, because the code keeps overwriting the same variable, if you run the whole script and then separate chunks to try some things, `running_time` will keep getting incremented. Once again, you have to be extra careful and keep track of it, once again overburdening you mentally.

## The solution

The solution to this problem looks like one of the previous things we tried, namely:

``````tic_filter <- function(...){

tic <- Sys.time()

result <- filter(...)

toc <- Sys.time()

running_time <- toc - tic

list("result" = result,
"running_time" = running_time)

}``````

While it is true that it returns a list, this function has the yuge advantage of being pure. But still, we need to solve two problems:

• how to avoid having to rewrite every function;
• how to compose these functions so that the output of one function can be ingested as the input of the next.

Solving the first problem consists in writing a new function that builds functions, what Hadley Wickham calls function factories. Let’s try:

``````timeit <- function(.f, ..., running_time = 0){

function(..., .running_time = running_time){

tic <- Sys.time()

result <- .f(...)

toc <- Sys.time()

list(result = result,
running_time = toc - tic + .running_time)
}

}``````

`timeit()` is a function that takes a function (and its arguments as an input), and returns a new function. This function returns the result of the original function (`.f`) evaluated on its arguments (`...`) as well as the time it took to run as a list. You’ll notice as well that this function takes another argument, called `running_time` with a default value of 0. This will become useful below, for now, ignore it.

``````t_sqrt <- timeit(sqrt)

t_sqrt(10)``````
``````## \$result
## [1] 3.162278
##
## \$running_time
## Time difference of 8.34465e-06 secs``````

That’s great, but we can’t compose these functions. This fails:

``````t_log <- timeit(log)

10 |>
t_sqrt() |>
t_log()``````
``Error in .f(...) : non-numeric argument to mathematical function``

because `t_log()` expects a number, not a list. The solution? Write another functions to help! Let’s call this function bind:

``````bind <- function(.l, .f, ...){

.f(.l\$result, ..., .running_time = .l\$running_time)

}``````

`bind()` takes a list object returned by a timed function (`.l`, with elements `\$result` and `\$running_time`) and applies another timed function `.f()` to the `\$result` element of `.l` as well as any further arguments `...` and finally sets the `running_time` argument of `.f` equal to `.l\$running_time`. `.l\$running_time` is the running time of the previous timed function call, so now this running time gets added to the running time of `.f` (see the definition of the list of `timeit()`).

An example might help:

``````t_log <- timeit(log)

10 |>
t_sqrt() |>
bind(t_log)``````
``````## \$result
## [1] 1.151293
##
## \$running_time
## Time difference of 8.368492e-05 secs``````

What’s nice with this solution, is that it works with any function:

``````t_filter <- timeit(filter)
t_select <- timeit(select)
t_group_by <- timeit(group_by)
t_summarise <- timeit(summarise)
t_p_longer <- timeit(tidyr::pivot_longer)

starwars %>%
t_filter(skin_color == "light") %>% # no need to use bind here
bind(t_select, species, sex, mass) %>%
bind(t_group_by, sex, species) %>%
bind(t_summarise,
total_individuals = n(),
min_mass = min(mass, na.rm = TRUE),
mean_mass = mean(mass, na.rm = TRUE),
sd_mass = sd(mass, na.rm = TRUE),
max_mass = max(mass, na.rm = TRUE)
) %>%
bind(t_select, -species) %>%
bind(t_p_longer, -sex, names_to = "statistic", values_to = "value")``````
``````## `summarise()` has grouped output by 'sex'. You can override using the `.groups`
## argument.``````
``````## \$result
## # A tibble: 10 × 3
## # Groups:   sex [2]
##    sex    statistic         value
##    <chr>  <chr>             <dbl>
##  1 female total_individuals   6
##  2 female min_mass           45
##  3 female mean_mass          56.3
##  4 female sd_mass            16.3
##  5 female max_mass           75
##  6 male   total_individuals   5
##  7 male   min_mass           79
##  8 male   mean_mass          90.5
##  9 male   sd_mass            19.8
## 10 male   max_mass          120
##
## \$running_time
## Time difference of 0.09293914 secs``````

There is some overhead compared to the solution that simply calls `tic` at the beginning of all the `{dplyr}` calls and then `toc` at the end, but this overhead becomes negligible the longer the base operations run for. And now the advantage is that you don’t have to think about keeping track of running times. Re-running separate chunks will also not interfere with the running time of any other chunk.

So here we are, ready to learn what monads are, or rather, we’re done, because you already know what monads are. The solution described before is a monad:

• a function factory to create functions that return a special, wrapped value (here it simply was a list of elements `\$result` and `\$running_time`). This wrapped value is also called a monadic value.
• a function to compose, or chain, these special functions together.

Some other pieces can be added to the list, and one would need to check so-called monadic laws to make extra sure we’re dealing with a monad, but that’s outside the scope of this blog post.

There are many monads, for instance the so-called `Maybe` monad, available on R thanks to Andrew McNeil who implemented this monad as an R package. I have also developed a monad for logging (which also logs execution time), which I called `{chronicler}`, read more about it here.

If this blog post was not enough to satiate your curiosity, here are some more nice resources:

Hope you enjoyed! If you found this blog post useful, you might want to follow me on twitter for blog post updates and buy me an espresso or paypal.me, or buy my ebook on Leanpub. You can also watch my videos on youtube. So much content for you to consoom!