---
title: "Package rvec"
author:
  - name: John Bryant
output:
  bookdown::html_document2:
    base_format: rmarkdown::html_vignette
    fig_caption: yes
    toc: true
    toc_depth: 2
    number_sections: true
pkgdown:
  as_is: true
vignette: >
  %\VignetteIndexEntry{Package rvec}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: references.bib
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```


# Introduction {#sec:intro}

The `rvec` package provides tools for working with random draws. The draws are held in a structure called an rvec. An rvec can, for many purposes, be treated like an ordinary R vector, and manipulated using standard functions from base R and the [tidyverse](https://tidyverse.org/). `rvec` also provides tools for summarising across draws to, for instance, produce point estimates or credible intervals.

# Examples {#sec:examples}

We start with some examples. Along with the `rvec` package, we will need some tidyverse packages,

```{r, message = FALSE}
library(rvec)
library(dplyr)
library(tidyr)
library(ggplot2)
```

When `rvec` is loaded, R emits messages about functions such as `sd()` and `var()`
being "masked". However, these functions retain their usual behavior for ordinary
vectors, while gaining new behavior for rvecs. No action is needed.

We create an rvec called `theta` containing three draws,

```{r}
l <- list(c(3, 1, 0))
theta <- rvec(l)
theta
```

The header `<rvec_dbl<3>[1]>` describes the structure of `theta`:

- `_dbl` indicates that `theta` is composed of [doubles](https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/double)
- `<3>` indicates that `theta` holds three random draws, and
- `[1]` indicates that each draw has length 1.

We can perform standard mathematical operations on `theta`,

```{r}
theta^2 + 1
```

`theta` recycles to match the length of other vectors,

```{r}
beta <- theta + c(1, -1)
beta
```

The new rvec `beta` consists of three draws of a vector of length 2:

1. the first vector, `c(4, 2)`, is obtained by adding `3` to `c(1, -1)`,
2. the second vector, `c(2, 0)`, is obtained by adding `1` to `c(1, -1)`, and
3. the third vector, `c(1, -1)`, is obtained by adding `0` to `c(1, -1)`.

To summarize across random draws, we use `draws_*` functions, e.g.,

```{r}
draws_mean(beta)
```
or 
```{r}
draws_sd(beta)
```

Our next example is based on draws from the posterior distribution from a Bayesian analysis of divorce in New Zealand. The posterior distribution refers to divorces per thousand people per year, disaggregated by age and sex. In their original form, the draws are not stored as an rvec, but instead in a one-draw-per-row format,

```{r}
divorce
```

We use function `collapse_to_rvec()` to convert from the one-draw-per-row format to an rvec,

```{r}
divorce_rvec <- divorce |>
  collapse_to_rvec(values = rate)
divorce_rvec
```

`collapse_to_rvec()` ensures that, internally, draw 1 for females aged 15-19 lines up with draw 1 for females aged 20-24, which lines up with draw 1 for females aged 25-29, and so on.

In this example, the rvec has many draws, so, rather than showing each draw, the print method shows the median value followed by the 2.5% and 97.5% quantiles.

We calculate the ratio between female and male divorce rates,

```{r}
divorce_ratio <- divorce_rvec |>
  pivot_wider(names_from = sex, values_from = rate) |>
  mutate(ratio = Female / Male)
divorce_ratio
```

Note that although `rate`, `Female`, and `Male` in these calculations are all rvecs, the code is identical to the code that would be needed for an ordinary R vector.

To obtain point estimates and uncertainty measures for rvecs, we use the `draw` functions. `draws_ci()`, for instance, returns medians and 95% credible intervals.

```{r}
divorce_ratio  <- divorce_ratio |>
  mutate(draws_ci(ratio))
divorce_ratio
```

Next we graph the results,

```{r, fig.width = 7, fig.height = 4}
ggplot(divorce_ratio,
       aes(x = age, 
           ymin = ratio.lower, 
           y = ratio.mid,
           ymax = ratio.upper)) +
  geom_pointrange() +
  ylab("Ratio") +
  ggtitle("Ratio between female divorce rate and male divorce rate")
```

# Structure of an rvec {#sec:structure}

The class `"rvec"` has four subclasses:

- `"rvec_dbl"`, which holds doubles, e.g. `3.142`, `-1.01`;
- `"rvec_int"`, which holds integers, e.g. `42`, `-1`;
- `"rvec_lgl"`, which holds `TRUE`, `FALSE`, and `NA`; and
- `"rvec_chr"`, which holds characters, e.g. `"a"`, `"Thomas Bayes"`.

Internally, an rvec is a matrix, with each row representing one unknown quantity, and each column representing one draw from the joint distribution of the unknown quantities,

Table: **The internal structure of an rvec**

|              | Draw 1        | Draw 2        | $\dots$  | Draw $n$      |
|--------------|:-------------:|:-------------:|:--------:|:-------------:|
| Quantity 1   | $\theta_{11}$ | $\theta_{12}$ | $\dots$  | $\theta_{1n}$ |
| Quantity 2   | $\theta_{21}$ | $\theta_{22}$ | $\dots$  | $\theta_{2n}$ |
| $\vdots$     | $\vdots$      | $\vdots$      | $\ddots$ | $\vdots$      |
| Quantity $m$ | $\theta_{m1}$ | $\theta_{m2}$ | $\dots$  | $\theta_{mn}$ |


Ordinary functions are applied independently to each column. For instance, calling `sum()` on an rvec creates a new rvec with structure

Table : **The result of summing along the vector**

|              | Draw 1                    | Draw 2                    | $\dots$  | Draw $n$                  |
|--------------|:-------------------------:|:-------------------------:|:--------:|:-------------------------:|
| Quantity 1   | $\sum_{i=1}^m\theta_{i1}$ | $\sum_{i=1}^m\theta_{i2}$ | $\dots$  | $\sum_{i=1}^m\theta_{in}$ |


Functions with a `draws_` prefix are applied independently to each row. For instance, calling `draws_mean()` on an rvec creates a new numeric vector with structure

Table: **The result of taking means across draws**

|              | Value                                |
|--------------|:------------------------------------:|
| Quantity 1   | $\frac{1}{n}\sum_{j=1}^n\theta_{1j}$ |
| Quantity 2   | $\frac{1}{n}\sum_{j=1}^n\theta_{2j}$ |
| $\vdots$     | $\vdots$                             |
| Quantity $m$ | $\frac{1}{n}\sum_{j=1}^n\theta_{mj}$ |


Each rvec holds a fixed number of draws. Two rvecs can only be used together in a function if

1. both rvecs contain the same number of draws, or
2. one of the rvecs contains a single draw.



# Creating rvecs

An individual rvec can be created from a list of vectors,
```{r}
x <- list(LETTERS, letters)
rvec(x)
```

a matrix, 
```{r}
x <- matrix(rnorm(2000), nrow = 2)
rvec(x)
```

or an atomic vector
```{r}
x <- c(TRUE, FALSE)
rvec(x)
```

Function `rvec()` chooses from classes `"rvec_dbl"`, `"rvec_int"`, `"rvec_lgl"`, and `"rvec_chr"`, based on the input. To enforce a particular choice of class, use function `rvec_dbl()`, `rvec_int()`, `rvec_lgl()`, or `rvec_chr()`,

```{r}
x <- list(1:3)
rvec(x)
rvec_dbl(x)
rvec_chr(x)
```
When the raw data take the form of a data frame with one draw per row, the most efficient way to create rvecs is to use `collapse_to_rvec()`. See Section \@ref(sec:examples) for an example.

Section \@ref(sec:prob) shows how to create an rvec consisting of draws from a standard probability distribution.

# Mathematical and logical operations

Mathematical and logical operations are applied independently to each draw. 
```{r}
x <- rvec(list(c(TRUE, FALSE),
               c(TRUE, TRUE)))
x	       
all(x)
any(x)
```
User-defined functions that consist entirely of standard mathematical and logical operations should work with no modifications.

```{r}
logit <- function(p) log(p / (1-p))
tibble(
  x = rvec(list(c(0.2, 0.4),
                c(0.6, 0.9))),
  y = logit(x)
)
```

Multiplying an rvec by a matrix produces an rvec (though only with R version 4.3.0 and higher)
```{r}
if (getRversion() >= "4.3.0") {
  m <- rbind(c(1, 1),
             c(0, 1))
  x <- rvec(list(1:2,
                 3:4))
  m %*% x
}
```

`rvec` contains a suite of functions for summarising weighted data, such as `weighted_mean()` and `weighted_var()`.

# Probability distributions {#sec:prob}

Standard R probability functions such as `dnorm()` or `rbinom()` do not allow rvec arguments. Package `rvec` provides modified functions that do. For instance,
```{r}
y <- rvec(list(c(-1, 0.2),
               c(3, -7)))
mu <- rvec(list(c(0, 1),
                c(0, -1)))
dnorm_rvec(y, mean = mu, sd = 3)
rbinom_rvec(n = 2, size = round(y+10), prob = 0.8)
```

The return value from an `rvec` probability function is an rvec if and only if at least one argument to the function is an rvec, with one exception. The exception is random variate functions, where a value can be supplied for a special argument called  `n_draw`. When a value for `n_draw` is supplied, the return value is an rvec with `n_draw` draws,
```{r}
rnorm_rvec(n = 3, mean = 100, sd = 10, n_draw = 2)
```
This is a convenient way to create inputs to a simulation. 


# Manipulating rvecs

## Compatibility with base R

Most code for manipulating ordinary R vectors should continue to work when applied to rvecs. There are, however, important base R functions where things go wrong. In each case, however, there are alternatives.

| Base R function(s) | What goes wrong | Alternative |
|:------------|:-------------|:----------------------|
| `rbind()`  | Creates a list matrix | `dplyr::bind_rows` or `vctrs::vec_rbind()`  |
| `cbind()` | Creates a list matrix | `dplyr::bind_cols` or `vctrs::vec_cbind()` |
| `ifelse()` | Doesn't preserve the rvec structure. | `rvec::if_else_rvec()`; `dplyr::if_else()` also works when `condition` is an ordinary vector |
| `sapply()`, `vapply()` | Doesn't combine results into an rvec | `rvec::map_rvec()` or `lapply()`  |

## Subsetting

Standard R ways of selecting elements from vectors work with rvecs. 
```{r}
x <- rvec(list(a = 1:2,
               b = 3:4,
               c = 5:6))
x[1]  ## element number
x[c("a", "c")]  ## element name
x[c(TRUE, FALSE, TRUE)]  ## logical flag
```


## If-Else {#sec:ifelse}

The standard R function `ifelse()` does not preserve the structure of an rvec.

The tidyverse function `if_else()` works when the `true`, `false`, or `missing` arguments are rvecs,
```{r}
x <- rvec(list(1:2,
               3:4))
if_else(condition = c(TRUE, FALSE), 
        true = x,
        false = -x)        
```

However, `if_else()` does not work when the `condition` argument is an rvec. For this we need `rvec` function `if_else_rvec()`,
```{r}
if_else_rvec(x <= 2, x, 2)
```

Function `if_else_rvec()` can be used to independently transform or recode values across different draws,
```{r}
x <- rvec(list(c(1, 3.3),
               c(NA, -2)))
x
x_recode <- if_else_rvec(is.na(x), 99, x)
x_recode
```

## Combining

The standard R concatenation function `c()` works with rvecs,
```{r}
x1 <- rvec(list(c(0.1, 0.2),
                c(0.3, 0.4)))
x2 <- rvec(list(c(0.5, 0.6),
                c(0.7, 0.8)))
c(x1, x2)
```

Unfortunately, `cbind()` and `rbind()` cannot be made to work properly on raw rvecs,
```{r}
rbind(x1, x2)
cbind(x1, x2)
```
though `cbind()` does work if the rvecs are contained in data frames
```{r}
df1 <- data.frame(x1)
df2 <- data.frame(x2)
cbind(df1, df2)
```

With rvec columns in data frames, `dplyr::bind_rows()` and
`vctrs::vec_rbind()` combine rows. To create columns from bare rvecs, use
`dplyr::bind_cols()` or `vctrs::vec_cbind()`,
```{r}
library(vctrs, warn.conflicts = FALSE)
vec_cbind(a = x1, b = x2)
```

When each result is an rvec, base R's `sapply()` does not combine them into
one rvec (`simplify = FALSE` keeps a list). `map_rvec()`, based on map
functions in package [purrr](https://purrr.tidyverse.org), combines
length-one rvec results:
```{r}
l <- list(a = rvec(list(c(1, 4))),
          b = rvec(list(c(9, 16))))
l
map_rvec(l, sqrt)
```
## Coercing

Function `as.matrix()` returns the data underlying an rvec.
```{r}
m <- matrix(1:6, nr = 2)
m
x <- rvec(m)
x
as.matrix(x)
```

Function `as_list_col()` returns a list of vectors
```{r}
as_list_col(x)
```
Functions such as [point_interval](https://mjskay.github.io/ggdist/reference/point_interval.html) in package [ggdist](https://mjskay.github.io/ggdist/) accept lists of vectors. A good way to access the facilities for working with random draws in `ggdist`, or in packages such as [tidybayes](https://mjskay.github.io/tidybayes/) and [bayesplot](https://mc-stan.org/bayesplot/),  is to use `as_list_col()` to convert an rvec to a list column.

Function `expand_from_rvec()` is the inverse of function `collapse_to_rvec()`, introduced in Section \@ref(sec:examples).
```{r}
divorce |>
  head(2)
divorce |>
  collapse_to_rvec(values = rate) |>
  head(2)
divorce |>
  collapse_to_rvec(values = rate) |>
  expand_from_rvec() |>
  head(2)
```

# Summarising distributions

Most functions in `rvec` are concerned with deriving random vectors from other random vectors: that is, with the column-wise calculations described in Section \@ref(sec:structure). But once we have derived the random vectors, we typically want to summarise them, using statistics such as means or quantiles: that is, we want to carry out row-wise calculations.

The functions for carrying out row-wise calculations on rvecs include, for instance, `draws_median()` to calculate the median values across draws, and `draws_var()` to calculate variance across draws.

Most draws functions use [matrixStats](https://CRAN.R-project.org/package=matrixStats) internally and are therefore fast.


# Similar packages

The first R package to provide a specialized object for handling multiple draws was [rv](https://CRAN.R-project.org/package=rv). The specialized object, called an rv, can be manipulated and summarized much like an rvec. However, in software terms, an rv is  not strictly a vector and does not behave like one inside a data frame. It is therefore not well suited to tidyverse-style work flows.

R package [posterior](https://CRAN.R-project.org/package=posterior) provides several data structures for handling multiple draws. One of these -- the rvar -- is similar to an rvec. An rvar is, however, not limited to a single dimension, and has special facilities for dealing with multiple chains (as produced by Markov chain Monte Carlo methods). These features are essential for some analyses, but they can make rvars harder to master, and they are not needed for most tidyverse-style work flows. 

Another important difference between rvecs and rvars is the way that they implement standard functions such as `mean()` and `sum()`. Calling `mean()` on an rvec yields the mean across elements within each draw. Means across draws within each element are obtained by calling `draws_mean()`. In contrast, calling `mean()` on an rvar yields the mean across draws within each element. Means across elements with each draw are obtained by calling `rvar_mean()`.

The differences in implementation imply that code written for an ordinary R vector generalises naturally to rvecs but not to rvars. Code such as `x <- median(y)` that was written originally for a numeric vector will continue to yield medians over elements of `y` if `y` is an rvec. If `y` is an rvar, however, then the code will yield medians over draws, which is a substantial change in meaning.